MODULE TheFunctions
  IMPLICIT NONE

CONTAINS
  ! First Internal Function
  REAL Function JJ(MatxC, VecB, X)
    REAL, DIMENSION(:)   :: X, VecB
    REAL, DIMENSION(:,:) :: MatxC
    INTEGER              :: N
    
    N = SIZE(X)
    jj       = - 0.5*(N+1)*SUM( X*MATMUL(MatxC,X)) - SUM(VecB*X)
  END Function JJ

  ! third Internal Function, use the JJ function.
  REAL Function MyF(MatxC, VecB, MyVecU,  Myk, rho)
    REAL, DIMENSION(:)            :: MyVecU, VecB
    REAL, DIMENSION(:,:)          :: MatxC
    REAL, DIMENSION(SIZE(MyVecU)) :: MyEk
    REAL    :: Rho
    INTEGER :: Myk

    MyEk      = 0.0
    MyEk(Myk) = rho

    MyF = JJ(MatxC, VecB, MyVecU + MyEk)

  END Function MyF

  ! computing the minimum in a given direction
  REAL Function Dichotomous(MatxC, VecB,  Myeps, Myk, MyVecU, MyA, MyB, idmax)
    REAL, DIMENSION(:)            ::  VecB
    REAL, DIMENSION(:,:)          :: MatxC
    REAL    :: Myeps
    INTEGER :: Myk, idmax
    REAL    :: MyA, MyB, Rho1, Rho2
    INTEGER :: id
    REAL, DIMENSION(:) :: MyVecU

    DO id = 1, idmax
       Rho1  = ( Mya + Myb - Myeps)/2.0
       Rho2  = ( Mya + Myb + Myeps)/2.0
       IF(       MyF(MatxC, VecB, MyVecU, Myk, rho1)  &
            & <= MyF(MatxC, VecB, MyVecU, Myk, rho2) ) THEN
          Myb = rho2
       ELSE
          Mya = Rho1
       END IF

       IF( ABS(Myb-Mya) < 2*Myeps ) Exit
    END DO
    Dichotomous = 0.5*( Mya + Myb)  
  END Function Dichotomous



  REAL Function GoldSection(MatxC, VecB,  Myeps, Myk, MyVecU, MyA, MyB, idmax)
    REAL, DIMENSION(:)            ::  VecB
    REAL, DIMENSION(:,:)          :: MatxC
    REAL    :: Myeps
    INTEGER :: Myk, idmax
    REAL, DIMENSION(:) :: MyVecU
    REAL    :: MyA, MyB, Rho1, Rho2, Gold
    INTEGER :: id

    Gold = (1.0+SQRT(5.0))/2.0

    DO id = 1, idmax
       Rho1  = ( Gold*Mya + Myb       )/(1.0+Gold)
       Rho2  = ( Mya      + Gold*Myb  )/(1.0+Gold)
       IF(       MyF(MatxC, VecB, MyVecU, Myk, rho1) &
            & <= MyF(MatxC, VecB, MyVecU, Myk, rho2) ) THEN
          Myb = rho2
       ELSE
          Mya = Rho1
       END IF

       IF( ABS(Myb-Mya) < 2*Myeps ) Exit

    END DO

    GoldSection = 0.5*( Mya + Myb)

  END Function GoldSection

END MODULE TheFunctions

PROGRAM MoreElaborated
  ! ========================================================================
  ! Opening the program CONTEXT
  ! ========================================================================

  ! This is to use the objects defined in the module Thefunctions
  USE TheFunctions

  IMPLICIT NONE   !  disable the implicit declarations

  ! ========================================================================
  ! SECTION : Declarations and specifications of Variables
  ! ========================================================================

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

  ! The end of the Declarations and specifications Section
  ! is the first executable instruction. In this case it is 
  ! The line with OPEN(UNIT=10, FILE="DataFile.data")

  ! ========================================================================
  ! SECTION :  Executable instructions  
  ! ========================================================================
  ! 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   

  ! READ(5,*)  N ! to read directly on the keyboard
  ! In this case OPEN ans CLOSE sections are remouved.

  ! 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))

  VecB  = VecB0
  VecU  = 10
  eps   = 1.0E-9
  a     = -0.5 
  b     = 1.5
  k     = 1
  write(6,*) Dichotomous(MatxC, VecB0, eps, k, VecU, a, b, 200*N)
  k     = 2
  write(6,*) Dichotomous(MatxC, VecB0, eps, k, VecU, a, b, 200*N)

  VecU  = 10 ! This is the starting point
  ! Save the evolution of minimizing point in the file "Nuage.gnu"
  OPEN(UNIT=11, FILE="Nuage.gnu")
  DO i = 1, 100*N
     DO k = 1, N
        a= -0.5 ; b = 1.5
        VecU(k) = VecU(k) + Dichotomous(MatxC, VecB0, eps, k, VecU, a, b, 200*N)
        IF(k<=2) WRITE(11, *) VecU(1:2), JJ(MatxC, VecB0, VecU)
     END DO

     IF(MOD(i,10) == 0)  Write(6,*) i, " Xs(1:2) = " , VecU(1:2)
  
  END DO

  VecB  = VecB0
  VecU  = -10! This is the starting point
  ! Save the evolution of minimizing point in the file "Nuage2.gnu"
  OPEN(UNIT=11, FILE="Nuage2.gnu")
  DO i = 1, 100*N
     DO k = 1, N
        a= -0.5 ; b = 1. ! fifferent values for a and b.
        VecU(k) = VecU(k) + GoldSection(MatxC, VecB0, eps, k, VecU, a, b, 200*N)
        IF(k<=2) WRITE(11, *) VecU(1:2), JJ(MatxC, VecB0,VecU)
     END DO
     IF(MOD(i,10) == 0)  Write(6,*) i, " Xs(1:2) = " , VecU(1:2)
  
  END DO
CLOSE(11)
  Write(6,*) " This is the actual end of the program"

END PROGRAM MoreElaborated
