MODULE DefinedFunctions
  IMPLICIT NONE

CONTAINS
  !  Function That compute J(X) for the input matrix MatxC the vector VecB
  ! ----------------------------------------------------------------------
  REAL Function JJ(MatxC, VecB, X)
    REAL, DIMENSION(:)  , INTENT(IN)  :: X, VecB
    REAL, DIMENSION(:,:), INTENT(IN)  :: MatxC
    INTEGER   :: N
    
    N = SIZE(X)
    jj       = - 0.5*(N+1)*SUM( X*MATMUL(MatxC,X)) - SUM(VecB*X)
  END Function JJ

  !  Function That compute F(rho)=J(MyVecU + rho*MyDk) for the input 
  !  matrix MatxC the vector VecB and for given  vector MyDk and a real Rho
  ! -----------------------------------------------------------------------
  REAL Function MyF(MatxC, VecB, MyVecU,  MyDk, rho)
    REAL, DIMENSION(:)   , INTENT(IN)  :: MyVecU, VecB, MyDk
    REAL, DIMENSION(:,:) , INTENT(IN)  :: MatxC
    REAL                 , INTENT(IN)  :: Rho

    REAL, DIMENSION(SIZE(MyVecU)) :: MyEk

    MyEk      = rho* MyDk

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

  END Function MyF

  !  Function that compute the value  that realized the minimum of 
  !  F(rho)=J(MyVecU + rho*MyDk)  in the interval [ MyA, MyB]  by a 
  !  Dichotomous seach, with a maximum of idmax iterations and
  !  a cut of parameter Myeps
  ! -----------------------------------------------------------------------
  REAL Function Dichotomous(MatxC, VecB,  Myeps, MyDk, MyVecU, MyA, MyB, idmax)
    REAL, DIMENSION(:)   , INTENT(IN)  :: MyVecU, VecB, MyDk
    REAL, DIMENSION(:,:) , INTENT(IN)  :: MatxC
    REAL                 , INTENT(IN)  :: Myeps, MyA, MyB
    INTEGER              , INTENT(IN)  ::  idmax
    
    REAL    ::  Rho1, Rho2
    INTEGER :: id

    DO id = 1, idmax
       Rho1  = ( Mya + Myb - Myeps)/2.0
       Rho2  = ( Mya + Myb + Myeps)/2.0
       IF(       MyF(MatxC, VecB, MyVecU, MyDk, rho1)  &
            & <= MyF(MatxC, VecB, MyVecU, MyDk, 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




  !  Function that compute the value  that realized the minimum of 
  !  F(rho)=J(MyVecU + rho*MyDk)  in the interval [ MyA, MyB]  by a 
  !  Gold section seach, with a maximum of idmax iterations and
  !  a cut of parameter Myeps
  ! -----------------------------------------------------------------------
  REAL Function GoldSection(MatxC, VecB,  Myeps, MyDk, MyVecU, MyA, MyB, idmax)
   REAL, DIMENSION(:)   , INTENT(IN)  :: MyVecU, VecB, MyDk
    REAL, DIMENSION(:,:) , INTENT(IN)  :: MatxC
    REAL                 , INTENT(IN)  :: Myeps, MyA, MyB
    INTEGER              , INTENT(IN)  ::  idmax
    
    REAL    ::  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, MyDk, rho1) &
            & <= MyF(MatxC, VecB, MyVecU, MyDk, 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


  !  Function that compute the gradient of J at the point  MyVecU
  ! -----------------------------------------------------------------------
 FUNCTION GradJJ(MatxC, VecB, MyVecU) RESULT(GradJ)
    REAL, DIMENSION(:)    , INTENT(IN)     ::  VecB,  MyVecU
    REAL, DIMENSION(:,:)   , INTENT(IN)    :: MatxC
    REAL, DIMENSION(SIZE(VecB))            ::  GradJ ! Return vector of the function

   INTEGER              :: N
    
    N = SIZE(VecB)

    GradJ = - (N+1)*MATMUL(MatxC, MyVecU) - VecB

  END FUNCTION GradJJ

END MODULE DefinedFunctions
