PROGRAM Simple
  IMPLICIT NONE 
  INTEGER :: N, i, k, m, kk
  REAL    :: eps,  Pi, Ju, JuE, DJ
  REAL, DIMENSION(:), POINTER   :: VecX, VecEk, VecB0, VecU, VecB, GradJ
  REAL, DIMENSION(:,:), POINTER :: MatxC, MatxA, MatxL


  ! Opening and rewind the file "DataFile.data" and
  ! associated it to the unit number 10
  ! Units 5 (keeboard) and 6 (sceen) are reserved
  ! -------------------------------------------------
  OPEN(UNIT=10, FILE="DataFile.data")

  ! By this declaration data in "DataFile.data"  are assumed
  ! to be in formatted form (readable).

  ! read the first line of the file and go to the next line
  ! The values readed are of the type of Nx and Lx (INTEGER and REAL) 
  ! -------------------------------------------------------
  READ(10,*) N

  ! Close the Unit 1O and the associated file
  ! -----------------------------------------
  CLOSE(10)

  ! set the value of Pi to the result of the  operation
  !     at the right of equality symbol.
  ! ----------------------------------------------------
  Pi = 4.0*ATAN(1.0)

  ! ALLOCATE the vectors Coor and Var to the range 1 to  Nx
  ! =================================================-======
  ALLOCATE( VecX(1:N), VecEk(1:N), VecU(1:N) )
  ALLOCATE( VecB(1:N), VecB0(1:N), GradJ(N)  )
  ALLOCATE(MatxC(N,N), MatxA(N,N), MatxL(N,N) ) 


  ! The memory to store the components of these variables
  ! is now up to date : we can make operations on it.
  ! =====================================================


  ! Initilisation of the matrix C (MatxC)
  MatxC =0
  MatxC(1,1)   = -2.0  ; MatxC(1,2) = 1.0 
  MatxC(N,N-1) =  1.0  ; MatxC(N,N) = -2.0 
  DO i = 2, N -1
     MatxC(i,i-1) =  1.0  ; MatxC(i,i) = -2.0  ; MatxC(i,i+1) = 1.0       
  END DO

  ! Initilisation of the vector B0 (VecB0)
  VecB0    = 0.0
  VecB0(N) = REAL( (N+1))
  VecU     = 1.0/(N+1)



  ! Question 1.1
  ! -----------------------------------------------
  write(6,'(90("-"))')
  Write(6,*) " This is the solution of the question 1.1  ==> U=1.0/(N+1)"

  VecB  = VecB0
  

  VecX  = VecU 
  Ju    = JJNL(VecX)
  Eps   = 1.0E-02

  DO k = 1, N
     VecEk = 0.0  ; VecEk(k) = eps
     VecX  = VecU + VecEk
     JuE   =  JJNL(VecX)    
     DJ    = JuE - Ju 
     write(6,*) "k=",k, " J(U-eps*Ek) - J(U) = ", DJ
  END DO


  ! Question 2
  ! -----------------------------------------------
  write(6,'(90("-"))')
  Write(6,*) " This is the solution of questions 1.2-1.4 ==> U=1.0/(N+1)"
  VecB  = VecB0
  Eps   = 1.0E-01
  DO m = 1, 4
     VecX  = VecU 
     Ju    = 0.5*(N+1)*SUM( VecX*MATMUL(MatxC,VecX)) + SUM(VecB*VecX)
     write(6,'(70("-"))')
     Write(6,*) " Eps = ", Eps, Ju
     write(6,'(1("-"))')

     DO k = 1, N
        VecEk = 0.0  ; VecEk(k) = eps
        VecX  = VecU + VecEk
        JuE   = JJNL(VecX)
        GradJ(k)   = (JuE - Ju)/eps
        IF(k<3 .OR. k>N-2) write(6,*) "k=",k, " ( J(U-eps*Ek) - J(U) )/eps = ", GradJ(k)
     END DO

     write(6,*) " Norm of approximated GradJ(U) :", SQRT(SUM(GradJ*GradJ))
     write(6,*) "                      J(U)     = ",Ju

     Eps = Eps**2
  END DO


  write(6,'(90("-"))')
  Write(6,*) " This is the solution of the question 4 ==> U(i)= i/(N+1)"
  VecB  = VecB0
  Eps   = 1.0E-01
  
  DO i = 1, N
     VecU(i) = REAL(i)/(N+1)
  END DO


  DO m = 1, 10
     Eps = 10.0**(- m)
     VecX  = VecU 
     Ju    = JJNL(VecX)
     write(6,'(70("-"))')
     Write(6,*) " Eps = ", Eps
     write(6,'(1("-"))')

     DO k = 1, N
        VecEk = 0.0  ; VecEk(k) = eps
        VecX  = VecU + VecEk
        JuE   = JJNL(VecX)
        GradJ(k)   = (JuE - Ju)/eps
        IF(k<3 .OR. k>N-2) write(6,*) "k=",k, " ( J(U-eps*Ek) - J(U) )/eps = ", GradJ(k)
     END DO

     write(6,*) " Norm of approximated GradJ(U) = ", SQRT(SUM(GradJ*GradJ))
     write(6,*) "                      J(U)     = ",Ju

  END Do
     write(6,'(70("-"))')
     write(6,'(70("-"))')


 ! Open a file to an other unit 11
  ! (but it an be 10 because this unit is now free after close(10))
  ! To save the initial solution  in a formatted form.
  ! -----------------------------------------------------------
  OPEN(UNIT=11, FILE="InitVar.gnu")

  DO i = 1, N
     WRITE(11, *) 
  END DO

  Write(6,*) " This is the actual end of the program"
CONTAINS

  REAL Function JJ(X)
    REAL, DIMENSION(:) :: X
    jj = 0.5*(N+1)*SUM( X*MATMUL(MatxC,X)) + SUM(VecB*X)
  END Function JJ

  REAL Function JJNL(X)
    REAL :: Alpha =0.1, Beta=0.2
    REAL, DIMENSION(:) :: X
    jjNL= 0.5*(N+1)*SUM( X*MATMUL(MatxC,X) ) &
         & +  0.5*alpha*SUM( X*X*X) + SUM(VecB*X) + Beta*SUM( X*X)
  END Function JJNL

END PROGRAM Simple
