MODULE ddgbundle_mod  ! DDG-Bundle method

  USE r_precision, ONLY : prec  ! Precision for reals.
  IMPLICIT NONE

! MODULE ddgbundle_mod includes the following subroutines (S) and functions (F).

  PUBLIC :: &
       init_ddgb, &  ! S Initialization for DDG-bundle subroutine.
       ddgbundle     ! S DDG-bundle subroutine for derivative free, nonsmooth,
                     !    large-scale optimization. Contains:
                     !     S restar  Initialization and reinitialization.
                     !     F sclpar  Selection of scaling parameter.
  PRIVATE :: &
       lls_armijo, & ! S Armijo line search. 
       lls_cubic, &  ! S Wolfe line search. Contains:
                     !     F qint    Quadratic interpolation.
                     !     F cint    Cubic interpolation.
       aggre1, &     ! S Simplified discrete gradient aggregation.
       aggre2, &     ! S Discrete gradient aggregation.
       dgrad, &      ! S Computation of discrete gradient
       wprint, &     ! S Printout the error and warning messages.
       rprint        ! S Printout the results.


CONTAINS

!***********************************************************************
!*                                                                     *
!*     * SUBROUTINE init_ddgb *                                        *
!*                                                                     *
!*     Initialization for diagonal discrete gradient bundle subroutine *
!*     for large-scale unconstrained nonsmooth optimization.           *
!*                                                                     *
!***********************************************************************
      
  SUBROUTINE init_ddgb(iterm)
    USE param, ONLY : zero,half,one,small,large
    USE initializat, ONLY : &
         n, &         ! Number of variables.
         mcu, &       ! Maximum number of stored corrections, mcu >= 1.
         iprint, &    ! Printout specification (see initializat for more details).
         iscale, &    ! Selection of the diagonal matrix updating 
                      ! (see initializat for more details).
         tolf, &      ! Tolerance for change of function values.
         tolf2, &     ! Second tolerance for change of function values.
         tolb, &      ! Tolerance for the function value.
         tolg, &      ! Tolerance for the first termination criterion. 
         told, &      ! Tolerance for the diagonal elements. 
         xmax, &      ! Maximum stepsize, 1 < XMAX.
         eta, &       ! Distance measure parameter, eta >= 0. 
         epsl, &      ! Line search parameter, 0 < epsl < 0.25.
         mtesf, &     ! Maximum number of iterations with changes of
                      ! function values smaller than tolf.
         mit, &       ! Maximun number of iterations.
         mfe          ! Maximun number of function evaluations.

    IMPLICIT NONE

! Scalar Arguments
    INTEGER, INTENT(OUT) :: iterm    ! Cause of termination:
                                     !   0  - Everything is ok.
                                     !  -5  - Invalid input parameters.
 
! Initialization and error checking

    IF (iprint < -1) iprint  = 1     ! Printout specification. 
    iterm = 0
      
    IF (n <= 0) THEN
       iterm = -5
       CALL wprint(iterm,iprint,1)
       RETURN 
    END IF

    IF (epsl >= 0.25_prec) THEN
       iterm = -5
       CALL wprint(iterm,iprint,4)
       RETURN
    END IF

    IF (mcu <= 0) THEN
       iterm = -5
       CALL wprint(iterm,iprint,2)
       RETURN
    END IF

! Default values

    IF (mit   <= 0) mit      = 10000           ! Maximum number of iterations.
    IF (mfe   <= 0) mfe      = n*mit           ! Maximum number of function evaluations.
    IF (tolf  <= zero) tolf  = 1.0E-08_prec    ! Tolerance for change of function values.
    IF (tolf2 == zero) tolf2 = 1.0E+04_prec    ! Second tolerance for change of function values.
    IF (tolb  == zero) tolb  = -large + small  ! Tolerance for the function value.
    IF (tolg  <= zero) tolg  = 1.0E-06_prec    ! Tolerance for the termination criterion.
    IF (told  <= zero) told  = 0.01_prec       ! Tolerance for the diagonal elements.
    IF (xmax  <= zero) xmax  = 1.5_prec        ! Maximum stepsize.
    IF (eta   <  zero) eta   = half            ! Distance measure parameter
    IF (epsl  <= zero) epsl  = 1.0E-04_prec    ! Line search parameter,
    IF (mtesf <= 0) mtesf    = 10              ! Maximum number of iterations with changes 
                                               ! of function values smaller than tolf.
    IF (iscale > 6 .OR. iscale < 0) iscale = 4 ! Selection of the diagonal matrix updating method.

     
  END SUBROUTINE init_ddgb


!***********************************************************************
!*                                                                     *
!*     * SUBROUTINE ddgbundle *                                        *
!*                                                                     *
!*     Diagonal discrete gradient bundle subroutine.                   *
!*                                                                     *
!***********************************************************************
      
  SUBROUTINE ddgbundle(f,nit,nfe,noutit,iterm,strtim)
    USE param, ONLY : small,large,zero,half,one  ! Parameters.
    USE exe_time, ONLY : getime  ! Execution time.
    USE initializat, ONLY : &
         n, &         ! Number of variables.
         x, &         ! Vector of variables
         mcu, &       ! Maximum number of stored corrections, mcu >= 1.
         iprint, &    ! Printout specification (see initializat for more details).
         iscale, &    ! Selection of the diagonal matrix updating method
                      ! (see initializat for more details).
         tolf, &      ! Tolerance for change of function values.
         tolf2, &     ! Second tolerance for change of function values.
         tolb, &      ! Tolerance for the function value.
         tolg, &      ! Tolerance for the first termination criterion. 
         told, &      ! Tolerance for the diagonal elements of diag.
         xmax, &      ! Maximum stepsize, 1 < XMAX.
         eta, &       ! Distance measure parameter, eta >= 0. 
         epsl, &      ! Line search parameter, 0 < epsl < 0.25.
         mtesf, &     ! Maximum number of iterations with changes of
                      ! function values smaller than tolf.
         mit, &       ! Maximun number of iterations.
         mfe, &       ! Maximun number of function evaluations.
         time         ! Maximum time
    USE obj_fun, ONLY : &
         myf          ! Computation of the value of the objective function. 
    USE subpro, ONLY : &
         vdot, &      ! Dot product.
         vxdiag, &    ! Vector is multiplied by a diagonal matrix y:=diag*x.
         xdiffy, &    ! Difference of two vectors.
         copy, &      ! Copying of a vector.
         copy2        ! Copying of two vectors.


    IMPLICIT NONE

! Scalar Arguments
    REAL(KIND=prec), INTENT(OUT) :: &
         f            ! Value of the objective function.
    INTEGER, INTENT(OUT) :: & 
         noutit, &    ! Number of outer iterations.
         nit,nfe      ! Number of iterations and function evaluations.
    INTEGER, INTENT(OUT) :: &
         iterm        ! Cause of termination:
                      !   1  - The problem has been solved with desired accuracy.
                      !   2  - Changes in function values < tolf in mtesf
                      !        subsequent iterations.
                      !   3  - Changes in function value < tolf*small*MAX(|f_k|,|f_(k-1)|,1),
                      !        where small is the smallest positive number such that 
                      !        1.0 + small > 1.0.
                      !   4  - Number of function calls > mfe.
                      !   5  - Number of iterations > mit.
                      !   6  - Time limit exceeded. 
                      !   7  - f < tolb.
                      !   8  - Failure in attaining the demanded accuracy.
                      !  -1  - Two consecutive restarts.
                      !  -2  - Number of restarts > maximum number of restarts.
                      !  -3  - Failure in function calculations 
                      !        (assigned by the user).
                      !  -4  
                      !  -5  - Invalid input parameters.
                      !  -6  - Unspecified error.


! Local Arrays
    REAL(KIND=prec), DIMENSION(n) :: &
         xo, &        ! Previous vector of variables.
         dg, &        ! Discrete gradient of the objective function.
         dgp, &       ! Previous  discrete gradient of the ohjective function.
         dga, &       ! Aggregate discrete gradient.
         g, &         ! Direction for discrete gradient.
         s, &         ! Difference of current and previous variables.
         u, &         ! Difference of current and previous discrete gradients.
         d, &         ! Direction vector.
         diag         ! Diagonal matrix.
    REAL(KIND=prec), DIMENSION(n*mcu) :: &
         ssmatrix, &  ! Matrix whose columns are stored differences of 
                      ! variables^2
         sumatrix     ! Matrix whose columns are stored differences of 
                      ! variables and discrete gradients.

! Local Scalars
    REAL(KIND=prec) :: &
         tmpx, &
         delta, &     ! Outer iteration parameter.
         sl, &        ! Size parameter for discrete gradient.
         slout, &     ! Bound for the size parameter sl. 
         alfn, &      ! Locality measure.
         alfv, &      ! Aggregate locality measure.
         epsr, &      ! Line search parameter.
         dnorm, &     ! Euclidean norm of the direction vector.
         dgnorm, &     ! Euclidean norm of the aggregate discrete gradient.
         xnorm, &     ! Stopping criterion.
         pxnorm, &    ! Previous stopping criterion.
         p, &         ! Directional derivative.
         tmax, &      ! Maximum stepsize.
         t, &         ! Stepsize.
         theta, &     ! Correction parameter for stepsize.
         fo, &        ! Previous value of the objective.
         stu, &       ! stu = trans(s)*u. 
         utu, &       ! utu = trans(u)*u. 
         sts, &       ! sts = trans(s)*s. 
         gamma_old, & ! Previous scaling parameter.
         gamma, &     ! Scaling parameter.
         tmp          ! Auxiliary value. !hihu uusi
    INTEGER :: i,j, &
         ngrad, &     ! Number of discrete gradient calls.
         mcc, &       ! Current number of stored corrections.
         inew, &      ! Index for the circular arrays.
         iters, &     ! Null step indicator.
                      !   0  - Null step.
                      !   1  - Serious step.
         nnk, &       ! Consecutive null steps counter.
         ibun, &      ! Index for the circular arrays in bundle updating.
         mal, &       ! Current size of the bundle.
         neps, &      ! Number of consecutive equal stopping criterions.
         ntesf, &     ! Number of tests on function decrease.
         ncres, &     ! Number of restarts.
         nres, &      ! Number of consecutive restarts.
         nout         ! Auxilary printout specification.

! Intrinsic Functions
    INTRINSIC ABS,MAX,SQRT

! Computational Time
    REAL strtim,elapstim
       
! Parameters
    REAL(KIND=prec), PARAMETER :: &
         fmin    = -large, &        ! Smallest acceptable value of the function.
         lengthd = 1.0E+20_prec, &  ! Direction vector length (no need if you use Armijo line-search). 
         tmin    = 1.0E-12_prec, &  ! Minimum stepsize.
         slmin   = 1.0E-5_prec, &   ! Minimum size parameter for the discrete gradient.
         pwt     = 1.0E-8_prec, &   ! Parameter for discrete gradient, pwt = z(sl)*alpha^j*e_j.
         sigma   = 1.0E-1_prec      ! Acceleration multiplier for outer iteration parameter.

    INTEGER, PARAMETER :: &
         maxeps = 20, &             ! Maximum number of consecutive equal stopping criterions.
         maxnrs = 2000              ! Maximum number of restarts.
      
      
    IF (iprint > 3) WRITE (6,FMT='(1X,''Entry to DDG-Bundle:'')')
      
     
! Initialization

    nout   = 0
    nit    = 0
    nfe    = 0
    ntesf  = 0
    nres   = 1
    ncres  = -1
    neps   = 0
    iterm  = 0
    iters  = 1
    nnk    = 0
    alfn   = zero
    alfv   = zero
    stu    = one
    sts    = one
    utu    = one
      
    tmax   = xmax
    xnorm  = large

    epsr   = 0.25_prec+small
    IF (epsl+epsl >= epsr) THEN
       epsr = epsl+epsl + small
       IF (epsr >= half) THEN
          CALL wprint(iterm,iprint,-2)
       END IF
    END IF
            
! Outer iteration

    tmp = one/SQRT(REAL(n, kind=prec)) !
    ngrad = 0
    noutit = 0

    slout = 1.25E+0_prec
    delta = (REAL(n*n, kind=prec)) 
     
    outer_iteration: DO WHILE (delta > tolg*sigma) ! iterm is the latest cause of stopping
                                                   ! if it is positive everything is ok.

       sl = max(2.0E-1_prec * slout,slmin) ! 
       slout = sl
       noutit = noutit + 1
       iterm  = 0
       iters  = 1
       ntesf  = 0
       neps   = 0
       nres   = 1

! First direction for the discrete gradient
      
       DO i=1,n
          g(i)=tmp
       END DO

! Computation of the value and the discrete gradient of the objective
! function and the search direction for the first iteration

       CALL myf(n,x,f,iterm)
       nfe = nfe + 1

       CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm)
       IF (iterm /= 0) EXIT outer_iteration


       CALL restar(n,mcc,inew,ibun,iters,dgp,dg,nnk,alfv,alfn,gamma,d,diag,mal,ncres)
 
     
! Start of the inner iteration
            
       inner_iteration: DO


! Computational time

          IF (time > 0.0E+00) THEN
             CALL getime(elapstim)
             IF (elapstim-strtim > time) THEN
                iterm = 6
                EXIT outer_iteration  ! Time limit exeeded, real exit
             END IF
          END IF


! Computation of norms

          IF (iters > 0) THEN
             dgnorm = vdot(n,dg,dg)
             dnorm = SQRT(vdot(n,d,d))

             p = vdot(n,dg,d)
             
          ELSE
             dgnorm = vdot(n,dga,dga)
             dnorm = SQRT(vdot(n,d,d))

             p = vdot(n,dga,d)
          END IF

   
! Test on descent direction

          IF (p+small*SQRT(dgnorm)*dnorm <= zero) THEN
             nres = 0
          
          ELSE
             nres = nres + 1
             IF (nres == 1) THEN
                CALL wprint(iterm,iprint,-3)
             
                CALL restar(n,mcc,inew,ibun,iters,dgp,dg,nnk,alfv,alfn,gamma,d,diag,mal,ncres)

                IF (ncres > maxnrs) THEN
                   nout = maxnrs
                   iterm = -2
                   EXIT outer_iteration  ! real exit
                END IF
 
                CYCLE inner_iteration
             END IF
             nout = -1
             iterm = -1
             EXIT outer_iteration ! real exit 
          END IF
       
       
! Stopping criterion

          nit = nit + 1
          pxnorm = xnorm
          xnorm = -p + 2.0_prec*alfv
       
     
! Tests for termination

          IF(xnorm <= delta) THEN 
             
             iterm = 1
             EXIT inner_iteration
          END IF

          IF (nfe >= mfe) THEN  ! too many function calls
             nout = mfe
             iterm = 4
             EXIT outer_iteration  ! real exit
          END IF

      
          IF (nit >= mit) THEN  ! too many iterations
             nout = mit
             iterm = 5
             EXIT outer_iteration  ! real exit
          END IF

      
          IF (f <= tolb) THEN  ! unbounded minimum
             iterm = 7
             EXIT outer_iteration  ! real exit
          END IF
    
      
          IF (iters == 0) THEN
             IF (ABS(xnorm - pxnorm) <= small) THEN
                neps = neps + 1
          
                IF (neps > maxeps) THEN
                   iterm = 8
                   EXIT inner_iteration  ! try again with different delta
                END IF

             ELSE
                neps = 0
             END IF

          ELSE
             neps = 0
          END IF


          IF (pxnorm < xnorm .AND. nnk > 2) THEN
             CALL wprint(iterm,iprint,-4)
          END IF
      

          CALL rprint(n,noutit,nit,nfe,ngrad,x,f,xnorm,half*dgnorm+alfv,iterm,iprint)
      
     
!     Preparation of line search

          fo = f
       
          IF (iters > 0) THEN
             CALL copy2(n,x,xo,dg,dgp)
          END IF

          IF (dnorm > zero) tmax = MAX(xmax/dnorm,sl)
       

! Initial step size

          IF (iters > 0) THEN
             t=min(2.0_prec,tmax)
          ELSE
             t=min(1.0_prec,tmax)
          END IF
     
     
! Armijo line search
          CALL lls_armijo(n,x,dg,d,xo,t,fo,f,p,alfn,sl,pwt,tmin,dnorm,xnorm, &
               epsl,epsr,eta,iters,nfe,ngrad,nnk,iterm)

! Line search with cubic interpolation
!          theta = one
!          IF (dnorm > lengthd) THEN
!             theta=lengthd/dnorm
!          END IF

!          CALL lls_cubic(n,x,dg,d,xo,t,fo,f,p,alfn,sl,pwt,tmin,dnorm,xnorm, &
!               theta,epsl,epsr,eta,iters,nfe,ngrad,nnk,iterm)
    

          IF (iterm /= 0) EXIT outer_iteration  ! real exit
       
          IF (tolf2 >= 0) THEN
             IF (ABS(fo-f) <= tolf2*small*MAX(ABS(f),ABS(fo),one) &
                  .AND. iters == 1) THEN
             
                iterm = 3
                EXIT inner_iteration
             END IF
          END IF

          IF (ABS(fo-f) <= tolf) THEN
             ntesf = ntesf + 1
             
             IF (ntesf >= mtesf .AND. iters == 1) THEN
                iterm = 2
                EXIT inner_iteration
             END IF
          
          ELSE
             ntesf = 0
          END IF
      

! Computation of variables difference 

          CALL xdiffy(n,x,xo,s)
  

! Computation of aggregate values and gradients difference

          IF (iters == 0) THEN
             nnk = nnk + 1

             IF (nnk == 1) THEN
                CALL xdiffy(n,dg,dgp,u)
             
                CALL aggre1(n,dg,dgp,dga,u,d,diag,alfn,alfv,gamma)
             
             ELSE
                CALL xdiffy(n,dg,dgp,u)
                CALL aggre2(n,dg,dgp,dga,d,diag,alfn,alfv,gamma)
             END IF
          
             CALL copy(n,xo,x)
             f = fo
          
          ELSE
             nnk = 0
             CALL xdiffy(n,dg,dgp,u)
          END IF

     
! Serious step initialization

          IF (iters > 0) THEN
             alfn = zero
             alfv = zero
          END IF

      
! Direction finding

          gamma_old = gamma

          IF (iters > 0) THEN  ! Serious step update and direction determination
       
             SELECT CASE(iscale)

             CASE(4)  ! Herskovits diagonal update 

                CALL update(n,mcu,mcc,inew,s,u,ssmatrix,sumatrix)  ! original with sts and stu

                DO i = 1, n
                   stu = zero
                   sts = zero
                   DO j = 1, mcc
                      stu = stu + sumatrix((j-1)*n+i)
                      sts = sts + ssmatrix((j-1)*n+i)
                   END DO
                   IF (stu > small) THEN  
                      diag(i) = sts/stu  
!                   IF (diag(i) < 1.0E-2_prec) diag(i) = 1.0E-2_prec
                   IF (diag(i) < 1.0E-2_prec) diag(i) = 1.0E-2_prec
                   IF (diag(i) > told) diag(i) = told
!                     IF (diag(i) < told) THEN
!                         diag(i) = told
!                      END IF
                   ELSE

                      diag(i) = one

                   END IF
                END DO

                CALL vxdiag(n,diag,dg,d)  ! 
                d = -one*d
                gamma = one


             CASE(5)  ! Complement of Herskovits's diagonal update 

                CALL update(n,mcu,mcc,inew,u,s,ssmatrix,sumatrix)  ! complement with utu and stu

                DO i = 1, n
                   stu = zero
                   sts = zero
                   DO j = 1, mcc
                      stu = stu + sumatrix((j-1)*n+i)
                      sts = sts + ssmatrix((j-1)*n+i)
                   END DO
                   IF (sts > small) THEN  ! utu > 0, complement with stu/utu
                      diag(i) = stu/sts 
                   IF (diag(i) < 1.0E-6_prec) diag(i) = 1.0E-6_prec
                   IF (diag(i) > told) diag(i) = told
!                      IF (diag(i) < told) THEN
!                         diag(i) = told
!                      END IF
                   ELSE

                      diag(i) = one

                   END IF
                END DO

                CALL vxdiag(n,diag,dg,d)  ! 
                d = -one*d
                gamma = one


             CASE(0,2)  ! scaling with stu/utu
                stu = vdot(n,s,u)

                IF (stu > zero) THEN  ! pos.def. matrix
                   utu = vdot(n,u,u)

                   gamma = sclpar(iscale,sts,stu,utu)
                   diag = gamma  

                ELSE
!                   diag = gamma_old   ! skipping
!                   gamma = gamma_old  !
                   diag = one         ! setting diag = I
                   gamma = one        !
                END IF

                d = -gamma*dg  ! diag = gamma*I

             CASE(1,3)  ! scaling with sts/stu
                stu = vdot(n,s,u)

                IF (stu > zero) THEN  ! pos.def. matrix
                   sts = vdot(n,s,s)

                   gamma = sclpar(iscale,sts,stu,utu)
                   diag = gamma  

                ELSE
!                   diag = gamma_old    !skipping
!                   gamma = gamma_old   !
                   diag = one          ! setting diag = I
                   gamma = one         !
                END IF

                d = -gamma*dg  ! diag = gamma*I


             CASE DEFAULT ! no scaling: diag = I
                          ! no need to update diag or gamma
                d = -one*dg ! 


             END SELECT

          ELSE  ! Update and direction determination for null steps.
     

             SELECT CASE(iscale)

             CASE(4,5)  ! Herskovits's diagonal update

                CALL vxdiag(n,diag,dga,d)  ! diag = diag at null steps.
                d = -one*d

!                diag = one  ! diag = 1 at null steps.
!                d = -one*dga


!                IF (nnk == 1) THEN  ! the first null step
!                   CALL update(n,mcu,mcc,inew,s,u,ssmatrix,sumatrix)  ! original with sts and stu

!                   DO i = 1, n
!                      stu = zero
!                      sts = zero
!                      DO j = 1, mcc
!                         stu = stu + sumatrix((j-1)*n+i)
!                         sts = sts + ssmatrix((j-1)*n+i)
!                      END DO
!                      IF (stu > small) THEN  ! stu > 0, original with sts/stu
!                         diag(i) = sts/stu  
!                         IF (diag(i) < told) diag(i) = told
!                      ELSE
!                         diag(i) = one  ! you could also skip this

!                      END IF
!                   END DO

!                END IF

!                CALL vxdiag(n,diag,dga,d)
!                d = -one*d


!                stu = vdot(n,s,u)  ! diag = gamma*I at null steps if gamma < gamma_old

!                IF (stu > zero) THEN  ! pos.def. matrix
!                   sts = vdot(n,s,s)

!                   gamma = sclpar(iscale,sts,stu,utu)  
!                   CALL update(n,mcu,mcc,inew,s,u,ssmatrix,sumatrix)  

!                   IF (gamma < gamma_old) THEN
!                      d = -gamma*dga
!                      diag = gamma
!                   ELSE 
!                      CALL vxdiag(n,diag,dga,d)  ! diag = diag
!                      d = -one*d
!                      gamma = gamma_old
!                   END IF

!                ELSE  ! stu <= 0
!                   CALL vxdiag(n,diag,dga,d)  ! diag = diag
!                   d = -one*d

!                END IF


             CASE(0,2)
! diag = 1 at null steps
!                diag = one
!                d = -one*dga

! diag = diag at null steps (skipping the update)
                d = -gamma*dga  ! no need to update diag or gamma

! diag = gamma*I at null steps if gamma < gamma_old
!                stu = vdot(n,s,u)

!                IF (stu > zero) THEN  ! pos.def. matrix
!                   utu = vdot(n,u,u)

!                   gamma = sclpar(iscale,sts,stu,utu)

!                ELSE
!                   gamma = gamma_old
!!                   gamma = one 

!                END IF

!                IF (gamma > gamma_old) gamma = gamma_old
!!                IF (gamma > gamma_old) gamma = one
!                d = -gamma*dga
!                diag = gamma


             CASE(1,3)
! diag = 1 at null steps
!                diag = one
!                d = -one*dga

! diag = diag at null steps (skipping the update)
                d = -gamma*dga  ! no need to update diag or gamma

! diag = gamma*I at null steps if gamma < gamma_old
!                stu = vdot(n,s,u)

!                IF (stu > zero) THEN  ! pos.def. matrix
!                   sts = vdot(n,s,s)

!                   gamma = sclpar(iscale,sts,stu,utu)

!                ELSE
!                   gamma = gamma_old
!!                  gamma = one 

!                END IF

!                IF (gamma > gamma_old) gamma = gamma_old
!!                IF (gamma > gamma_old) gamma = one
!                d = -gamma*dga
!                diag = gamma


             CASE DEFAULT ! no scaling: diag = I
           
                d = -one*dga  ! no need to update diag or gamma
             
             END SELECT
          END IF


       END DO inner_iteration
      
       delta  = MIN(delta*sigma,xnorm)

    END DO outer_iteration

     

! Printout the final results

    IF (iprint > 3) WRITE (6,FMT='(1X,''Exit from DDG-Bundle:'')')

    CALL wprint(iterm,iprint,nout)
    CALL rprint(n,noutit,nit,nfe,ngrad,x,f,xnorm,half*dgnorm+alfv,iterm,iprint)

  CONTAINS


    SUBROUTINE restar(n,mcc,inew,ibun,iters,dgp,dg,nnk,alfv,alfn,gamma,d,diag,mal,ncres)  ! Initialization
      
!      USE param, ONLY : zero,one  ! given in host
!      USE subpro, ONLY : copy  ! given in host
      IMPLICIT NONE

! Array Arguments
      REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
           dgp        ! Basic discrete gradient of the objective function.
      REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
           dg         ! Current (auxiliary) discrete gradient of the objective function.
      REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: &
           d, &      ! Search direction.
           diag      ! Diagonal varaible metric matrix.

! Scalar Arguments
      INTEGER, INTENT(IN) :: & 
           n         ! Number of variables
      INTEGER, INTENT(OUT) :: &
           mcc, &
           inew, &
           ibun, &   ! Index for the circular arrays in bundle updating.
           nnk, &    ! Consecutive null steps counter.
           mal       ! Current size of the bundle.
      INTEGER, INTENT(INOUT) :: &
           iters, &  ! Null step indicator.
                     !   0  - Null step.
                     !   1  - Serious step.
           ncres     ! Number of restarts.
      REAL(KIND=prec), INTENT(OUT) :: & 
           alfn, &   ! Locality measure.
           alfv, &   ! Aggregate locality measure.
           gamma     ! Scaling parameter.


! Restart
      mcc   = 0
      inew  = 1
      ibun  = 1
      mal   = 0
      ncres = ncres + 1
        
      IF (iters == 0) THEN
         CALL copy(n,dgp,dg)
         iters = 1
         nnk = 0
         alfv=zero
         alfn=zero
      END IF

!      gamma = large
      gamma = one
      d = -one*dg
      diag = one
      
    END SUBROUTINE restar
    

    FUNCTION sclpar(iscale,sts,stu,utu) RESULT(spar) ! Calculation of the scaling parameter for DDG-Bundle.
      
!    USE param, ONLY : small,one,half ! given in host
!    USE initializat, ONLY : told     ! tolerance for the diagonal elements. 
      IMPLICIT NONE

! Scalar Arguments
      REAL(KIND=prec), INTENT(IN) :: & 
           sts, &    ! sts = trans(s)*s. 
           stu, &    ! stu = trans(s)*u. 
           utu       ! utu = trans(u)*u. 
      REAL(KIND=prec) :: & 
           spar      ! Scaling parameter.
      INTEGER, INTENT(IN) :: & 
           iscale    ! Selection of the scaling:
                     !   0  - Scaling at every iteration with STU/UTU.
                     !   1  - Scaling at every iteration with STS/STU.
                     !   2  - Interval scaling with STU/UTU.
                     !   3  - Interval scaling with STS/STU.
                     !   6  - No scaling.

! Intrinsic Functions
      INTRINSIC SQRT


! Computation of scaling parameter.

      SELECT CASE(iscale)
     
! Scaling parameter = STU/UTU

      CASE(0,2)
         IF (utu < SQRT(small)) THEN
            spar = one
            RETURN
         ELSE
            spar = stu/utu
         END IF
    
! Scaling parameter = STS/STU

      CASE(1,3,4,5)
         IF (stu < SQRT(small)) THEN
            spar = one
            RETURN
         ELSE
            spar = sts/stu
         END IF

! No scaling

      CASE DEFAULT
         spar = one
         RETURN
      END SELECT

               
! Scaling
            
      SELECT CASE(iscale)

! Lower limit      
      CASE(0,1)

         IF (spar < one) spar = one ! 

         
! Interval scaling
      CASE(2)
!         IF (spar < 0.1_prec .OR. spar > 10.0_prec) spar = one !
         IF (spar < half .OR. spar > 5.0_prec) spar = one !

      CASE(3)
!         IF (spar < 0.1 .OR. spar > 10.0_prec) spar = one !
         IF (spar < half .OR. spar > 5.0_prec) spar = one !

! Heskovits diagonal (null step)
      CASE(4)
         IF (spar < told ) spar = told 

! Scaling at every iteration
      CASE DEFAULT
         CONTINUE

      END SELECT

        
    END FUNCTION sclpar
      
  END SUBROUTINE ddgbundle

    
!************************************************************************
!*                                                                      *
!*     * SUBROUTINE lls_armijo *                                        *
!*                                                                      *
!*     Special line search.                                             *
!*                                                                      *
!************************************************************************

  SUBROUTINE lls_armijo(n,x,dg,d,xo,t,fo,f,p,alfn,sl,pwt,tmin,dnorm,&
       wk,epsl,epsr,eta,iters,nfe,ngrad,nnk,iterm)

    USE param, ONLY : zero,half,one
    USE obj_fun, ONLY : &
         myf       ! Computation of the value of the objective function. 
    USE subpro, ONLY : &
         scsum, &  ! Sum of a vector and the scaled vector.
         vdot      ! Dot product of vectors.
    IMPLICIT NONE

! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
         d, &      ! Direction vector.
         xo        ! Previous vector of variables.
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
         x         ! Vector of variables.
    REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: &
         dg        ! Discrete gradient of the objective function.

! Scalar Arguments
    REAL(KIND=prec), INTENT(INOUT) :: & 
         epsl,&    ! Linesearch parameter.
         epsr,&    ! Linesearch parameter. 
         p, &      ! Directional derivative.
         t         ! Stepsize
    REAL(KIND=prec), INTENT(IN) :: & 
         sl, &
         pwt, &
         eta, &    ! Distance measure parameter.
         fo, &     ! Previous value of the objective function.
         wk, &     ! Stopping parameter.
         dnorm, &  ! Euclidean norm of the direction vector.
         tmin      ! Lower limit for stepsize parameter.
    REAL(KIND=prec), INTENT(OUT) :: & 
         f, &      ! Value of the objective function.
         alfn      ! Locality measure.
    INTEGER, INTENT(IN) :: & 
         n, &      ! Number of variables
         nnk       ! Number of consequtive null steps.
    INTEGER, INTENT(INOUT) :: & 
         nfe, &    ! Number of function evaluations.
         ngrad     ! Number of discrete gradient calls.
    INTEGER, INTENT(OUT) :: &
         iters, &  ! Null step indicator.
                   !   0  - Null step.
                   !   1  - Serious step.
         iterm     ! Cause of termination:
                   !   0  - Everything is ok.
                   !  -3  - Failure in function or subgradient calculations.

! Local arrays
    REAL(KIND=prec), DIMENSION(n) :: &
         g, &      ! Initial direction.
         xu        ! trial point.

! Local Scalars
    REAL(KIND=prec) :: &
         fu, &          ! Auxiliary values of the objective function.
         alpha, &       ! Increase parameter.
         a_orig, &      ! Original value of increase parameter. 
         epslwk         ! Auxiliary scalar.
    INTEGER :: i,nin, & ! Number of interpolations.
         maxnin

! Intrinsic Functions
    INTRINSIC ABS,MAX

      

! Initialization

    nin = 0
    IF (nnk >= 1) THEN
       maxnin = 40
!       maxnin = 50
    ELSE
       maxnin = 40
!       maxnin = 50
    END IF
    alpha  = 2.0_prec
    a_orig = alpha
    epslwk = epsl*wk


! Function evaluation at a new point

    DO i=1,n
       g(i)=d(i)/dnorm
    END DO

      
    CALL scsum(n,t,d,xo,x)
    CALL myf(n,x,f,iterm)
    nfe = nfe + 1
    IF (iterm /= 0) RETURN

    CALL scsum(n,alpha*t,d,xo,xu)
    CALL myf(n,xu,fu,iterm)
    nfe = nfe + 1
    IF (iterm /= 0) RETURN
       
    IF (f <= fo - t*epslwk) THEN ! descent condition

       IF(fu > fo - alpha*t*epslwk) THEN


! Serious step
          iters=1
          CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm)
          RETURN

       ELSE

! Increase t
          inct: DO
             nin=nin+1 
             f=fu
             x=xu
             alpha = alpha * a_orig
             CALL scsum(n,alpha*t,d,xo,xu)
             CALL myf(n,xu,fu,iterm)
             nfe = nfe + 1

             IF( fu > fo - alpha*t*epslwk .OR. nin > maxnin) THEN

! Serious step
                iters=1
                CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm)
                RETURN
             END IF
          END DO inct
       END IF
    ELSE

       alpha = one 
! Decrease t
       dect: DO
          nin=nin+1
          IF (nin > maxnin) EXIT dect
          alpha = alpha / a_orig ! 
          IF (alpha*t < tmin) THEN
             alpha = alpha * a_orig
             EXIT dect
          END IF
          CALL scsum(n,alpha*t,d,xo,x)
          CALL myf(n,x,f,iterm)
          IF (f <= fo - alpha*t*epslwk) THEN ! descent condition

! Serious step
             iters=1
             CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm)
             RETURN
          END IF
       END DO dect
    END IF

! Null step
    iters = 0


    CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm)
    IF (iterm /= 0) RETURN

    p = vdot(n,dg,d) 
!   alfn = MAX(ABS(fo-f+p*t),eta*(t*dnorm)**2) 
    alfn = MAX(ABS(fo-f+alpha*t*p),eta*(alpha*t*dnorm)**2) 


  END SUBROUTINE lls_armijo


!************************************************************************
!*                                                                      *
!*     * SUBROUTINE lls_cubic *                                         *
!*                                                                      *
!*     Special line search.                                             *
!*                                                                      *
!************************************************************************

  SUBROUTINE lls_cubic(n,x,dg,d,xo,t,fo,f,p,alfn,sl,pwt,tmin,dnorm,&
       wk,theta,epsl,epsr,eta,iters,nfe,ngrad,nnk,iterm)

    USE param, ONLY : zero,half,one
    USE obj_fun, ONLY : &
         myf       ! Computation of the value of the objective function. 
    USE subpro, ONLY : &
         scsum, &  ! Sum of a vector and the scaled vector.
         vdot      ! Dot product of vectors.
    IMPLICIT NONE

! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
         d, &      ! Direction vector.
         xo        ! Previous vector of variables.
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
         x         ! Vector of variables.
    REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: &
         dg        ! Discrete gradient of the objective function.

! Scalar Arguments
    REAL(KIND=prec), INTENT(INOUT) :: & 
         epsl,&    ! Linesearch parameter.
         epsr,&    ! Linesearch parameter.
         p, &      ! Directional derivative.
         t         ! Stepsize
    REAL(KIND=prec), INTENT(IN) :: & 
         sl, &
         pwt, &
         theta, &  ! Scaling parameter.
         eta, &    ! Distance measure parameter.
         fo, &     ! Previous value of the objective function.
         wk, &     ! Stopping parameter.
         dnorm, &  ! Euclidean norm of the direction vector.
         tmin      ! Lower limit for stepsize parameter.
    REAL(KIND=prec), INTENT(OUT) :: & 
         f, &      ! Value of the objective function.
         alfn      ! Locality measure.
    INTEGER, INTENT(IN) :: & 
         n, &      ! Number of variables
         nnk       ! Number of consequtive null steps.
    INTEGER, INTENT(INOUT) :: & 
         nfe, &    ! Number of function evaluations.
         ngrad     ! Number of discrete gradient calls.
    INTEGER, INTENT(OUT) :: &
         iters, &  ! Null step indicator.
                   !   0  - Null step.
                   !   1  - Serious step.
         iterm     ! Cause of termination:
                   !   0  - Everything is ok.
                   !  -3  - Failure in function or subgradient calculations.

! Local arrays
    REAL(KIND=prec), DIMENSION(n) :: &
         g         ! Initial direction.

! Local Scalars
    REAL(KIND=prec) :: &
         tl,tu, &  ! Lower and upper limits for t used in interpolation.
         fl,fu, &  ! Values of the objective function for t=tl and t=tu.
         pl, &     ! Previous directional derivative.
         epsa,epst,epslk,epsrk, & ! Line search parameters.
         thdnorm,epsawk,epstwk,epslwk,epsrwk ! Auxiliary scalars.
    INTEGER :: i,nin ! Number of interpolations.

! Parameters
    INTEGER, PARAMETER :: maxint = 200  ! Maximum number of interpolations.

! Intrinsic Functions
    INTRINSIC ABS,MAX

      

! Initialization

    nin = 0

    epst = epsl+epsl
    epsa = half*(epsr - epst)
    thdnorm = theta*dnorm

    pl = -p
    tl = zero
    tu = t
    fl = fo
    fu = fo

    IF (theta < one) THEN
       epst  = theta*epst
       epsa  = theta*epsa
       epslk = epsl
       epsl  = theta*epsl
       epsrk = epsr
       epsr  = theta*epsr
    END IF
 
    epsawk   = epsa*wk
    epslwk   = epsl*wk
    epsrwk   = epsr*wk
    epstwk   = epst*wk


! Function and gradient evaluation at a new point

    DO i=1,n
       g(i)=d(i)/dnorm
    END DO

    lls_iteration: DO
       
       CALL scsum(n,theta*t,d,xo,x)

       CALL myf(n,x,f,iterm)
       nfe = nfe + 1

       IF (iterm /= 0) RETURN

    
! Null/descent step test (ITERS=0/1)

       iters = 1
       IF (f <= fo - t*epstwk) THEN
          tl = t
          fl = f
       ELSE
          tu = t
          fu = f
       END IF
      

! Additional interpolation

       IF (f > fo - t*epslwk .AND. t > one) THEN
          nin=nin+1
          IF (tl == zero .AND. wk > zero) THEN
             t = qint(tu,fl,fu,pl,one-half/(one-epst)) ! If you are using cubic interpolation
          ELSE
             t = half*(tu+tl)
          END IF
          CYCLE lls_iteration
       END IF


! original additional interpolation
       IF (f > fo - t*epslwk .AND. tu-tl >= tmin*0.1_prec &
            .AND. nnk >= 1 .AND. nin < maxint) THEN

          nin=nin+1
          IF (tl == zero .AND. wk > zero) THEN
             t = qint(tu,fl,fu,pl,one-half/(one-epst)) ! If you are using cubic interpolation
          ELSE
             t = half*(tu+tl)
          END IF
          CYCLE lls_iteration
       END IF

       CALL dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm)
       IF (iterm /= 0) RETURN

       p = theta*vdot(n,dg,d)
       alfn = MAX(ABS(fo-f+p*t),eta*(t*thdnorm)**2)


! Serious step

       IF (f <= fo - t*epslwk .AND. (t >= tmin .OR. alfn > epsawk)) EXIT lls_iteration


! Null step

       IF (p-alfn >= -epsrwk .OR. tu-tl < tmin*0.1_prec .OR. &
            nin >= maxint) THEN
          iters = 0
          EXIT lls_iteration
       END IF


! Interpolation

       nin=nin+1
       IF (tl == zero .AND. wk > zero) THEN
          t = cint(tu,fl,fu,pl,p,one-half/(one-epst))
       ELSE
          t = half*(tu+tl)
       END IF

    END DO lls_iteration

    IF (theta /= one) THEN
       epsl = epslk
       epsr = epsrk
    END IF

  CONTAINS

    FUNCTION qint(tu,fl,fu,wk,kappa) RESULT(t) ! Quadratic interpolation

!      USE param, ONLY : half,one  ! given in host
      IMPLICIT NONE
      
! Scalar Arguments
      REAL(KIND=prec), INTENT(IN) :: & 
           fl, &  ! Value of the objective function.
           fu, &  ! Value of the objective function for t=tu.
           wk, &  ! Directional derivative.
           tu, &  ! Upper value of the stepsize parameter.
           kappa  ! Interpolation parameter.
      REAL(KIND=prec) :: & 
           t      ! Stepsize.

! Local Scalars
      REAL(KIND=prec) :: tmp1,tmp2

! Intrinsic Functions
      INTRINSIC MAX


      tmp1 = (fu-fl)/ (-wk*tu)

     
! Quadratic interpolation with one directional derivative

      tmp2 = 2.0_prec * (one - tmp1)

      IF (tmp2 > one) THEN
   
! Interpolation accepted
       
         t = MAX(kappa*tu,tu/tmp2)
         RETURN
      END IF
      
     
! Bisection

      t = half*tu
      
    END FUNCTION qint

    FUNCTION cint(tu,fl,fu,wk,pu,kappa) RESULT(t) ! Cubic interpolation

!      USE param, ONLY : half,one  ! given in host
      IMPLICIT NONE
      
! Scalar Arguments
      REAL(KIND=prec), INTENT(IN) :: & 
           fl, &  ! Value of the objective function.
           fu, &  ! Value of the objective function for t=tu.
           wk, &  ! Directional derivative.
           pu, &  ! Directional derivative at t=tu.
           tu, &  ! Upper value of the stepsize parameter.
           kappa  ! Interpolation parameter.
      REAL(KIND=prec) :: & 
           t      ! Stepsize.

! Local Scalars
      REAL(KIND=prec) :: tmp1,tmp2,tmp3,tmp4,tmp5,tmp6

! Intrinsic Functions
      INTRINSIC MAX,SQRT


! Cubic interpolation
      tmp1 = (fu-fl)/ (-wk*tu)
      tmp3 = pu/ (-wk)  
      tmp4 = tmp3 - 2.0_prec*tmp1 + one
      tmp5 = tmp3 - 3.0_prec*tmp1 + 2.0_prec
      tmp6 = tmp5*tmp5 - 3.0_prec*tmp4

      IF (tmp6 < zero) THEN 
         t = half*tu  ! bisection
         RETURN
      END IF

      tmp2 = tmp5 + SQRT(tmp6)

      IF (tmp2 > one) THEN
   
! Interpolation accepted
       
         t = MAX(kappa*tu,tu/tmp2)
         RETURN
      END IF
      
! Bisection

      t = half*tu
      
    END FUNCTION cint

  END SUBROUTINE lls_cubic



!************************************************************************
!*                                                                      *
!*     * SUBROUTINE update *                                            *
!*                                                                      *
!*     Matrix update.                                                   *
!*                                                                      *
!************************************************************************
    
  SUBROUTINE update(n,mcu,mcc,inew,s,u,ssmatrix,sumatrix)
      
    USE subpro, ONLY : &
         copy2     ! Copying of two vectors.
    IMPLICIT NONE

! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
         s, &         ! Difference of current and previous variables.
         u            ! Difference of current and previous subgradients.
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
         ssmatrix, &  ! Matrix whose columns are stored corrections.
         sumatrix     ! Matrix whose columns are stored subgradient differences.

! Scalar Arguments
    INTEGER, INTENT(IN) :: & 
         n, &         ! Number of variables
         mcu          ! Declared number of stored corrections.
    INTEGER, INTENT(INOUT) :: & 
         mcc, &       ! Current number of stored corrections.
         inew         ! Index for circular arrays.
      
! Local Arrays
    REAL(KIND=prec), DIMENSION(n) :: &
         ss,su     ! Componentwise multiplication of s*s and s*u.

! Local Scalars
    INTEGER :: i

     
! Compute new values

    DO i = 1, n
       ss(i)=s(i)*s(i)
       su(i)=s(i)*u(i)
    END DO
         
     
! Update ssmatrix and sumatrix

    CALL copy2(n,ss,ssmatrix((inew-1)*n+1:),su,sumatrix((inew-1)*n+1:))
           
    inew = inew + 1
    IF (inew > mcu) inew = 1
    IF (mcc < mcu) mcc = mcc + 1
          
  END SUBROUTINE update


!************************************************************************
!*                                                                      *
!*     * SUBROUTINE aggre1 *                                            *
!*                                                                      * 
!*     Computation of aggregate values after first null step.           *
!*                                                                      *
!************************************************************************

  SUBROUTINE aggre1(n,dg,dgp,dga,u,d,diag,alfn,alfv,gamma)

    USE param, ONLY : zero,half,one
    USE subpro, ONLY : &
         vxdiag, & ! Vector multiplied by a diagonal matrix.
         vdot      ! Dot product.

    IMPLICIT NONE

! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
         diag, &   ! Diagonal variable metric update.
         d, &      ! Direction vector.
         dg, &      ! Current (auxiliary) subgradient of the objective function.
         dgp, &     ! Previous subgradient of the objective function.
         u         ! Difference of trial and aggregate gradients.
    REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: &
         dga        ! Next aggregate subgradient of the objective function.

! Scalar Arguments
    REAL(KIND=prec), INTENT(OUT) :: & 
         alfv      ! Aggregate locality measure.
    REAL(KIND=prec), INTENT(IN) :: & 
         gamma, &  ! Scaling parameter.
         alfn      ! Locality measure.
    INTEGER, INTENT(IN) :: & 
         n         ! Number of variables

! Local arrays
    REAL(KIND=prec), DIMENSION(n) :: tmpn

! Local Scalars
    REAL(KIND=prec) :: &
         p, &      ! p = trans(d)*u - alfn.
         q, &      ! q = trans(u)*dm*u, where dm is the inverse approximation of 
                   ! the Hessian calculated by using the L-BFGS formula.
         lam, &    ! Multiplier used to calculate aggregate values.
         w         ! Correction.
    INTEGER :: i

! Intrinsic Functions
    INTRINSIC MAX,MIN,SIGN

      
! Computation of p=trans(d)*u - alfn and the product q=trans(u)*dm*u

    p = vdot(n,d,u) - alfn
    CALL vxdiag(n,diag,u,tmpn)
    q = vdot(n,u,tmpn)

    
    lam = half + SIGN(half,p)

    IF (q > zero) lam = MIN(one,MAX(zero,p/q))
      

! Computation of the aggregate values

    p = one - lam
    DO i=1,n
       dga(i)=lam*dg(i) + p*dgp(i)
    END DO
      
    alfv = lam*alfn
      
  END SUBROUTINE aggre1


!************************************************************************
!*                                                                      *
!*     * SUBROUTINE aggre2 *                                            *
!*                                                                      *
!*     Computation of aggregate values after consequtive null steps.    *
!*                                                                      *
!************************************************************************
      
  SUBROUTINE aggre2(n,dg,dgp,dga,d,diag,alfn,alfv,gamma)
      
    USE param, ONLY : zero,one,small
    USE subpro, ONLY : &
         vxdiag, & ! Vector multiplied by a diagonal matrix.
         vdot, &   ! Dot product.
         xdiffy    ! Difference of two vectors.
    IMPLICIT NONE

! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
         diag, &   ! Diagonal variable metric update.
         d, &      ! Direction vector.
         dg, &     ! Current (auxiliary) subgradient of the objective function.
         dgp       ! Previous subgradient of the objective function.
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
         dga       ! Aggregate subgradient of the objective function.

! Scalar Arguments
    REAL(KIND=prec), INTENT(INOUT) :: & 
         alfv      ! Aggregate locality measure.
    REAL(KIND=prec), INTENT(IN) :: & 
         gamma, &  ! Scaling parameter.
         alfn      ! Locality measure.
    INTEGER, INTENT(IN) :: & 
         n         ! Number of variables
      
! Local arrays
    REAL(KIND=prec), DIMENSION(n) :: tmpn2,tmpn3,tmpn4

! Local Scalars
    REAL(KIND=prec) :: &
         pr, &     ! pr = trans(dgp-dga) dm (dgp-dga), where dm
                   ! presents the L-SR1- approximation of Hessian.
         rrp, &    ! rrp = trans(dgp-dga) dm dga - alfv.
         prqr, &   ! prqr = trans(dgp-dga) dm (dg-dga).
         rrq, &    ! rrq = trans(dg-dga) dm dga - alfv + alfn.
         qr, &     ! qr = trans(dg-dga) dm (dg-dga).
         pq, &     ! pq = trans(dg-dgp) dm (dg-dgp).
         qqp, &    ! qqp = trans(dg-dgp) dm dg + alfn.
         lam1, &   ! Multiplier used to calculate aggregate values.
         lam2, &   ! Multiplier used to calculate aggregate values.
         tmp1, &   ! Auxiliary scalar.
         tmp2, &   ! Auxiliary scalar.
         tmp3      ! Auxiliary scalar.
    INTEGER :: i 

      
! Intrinsic Functions
    INTRINSIC MIN,MAX


    CALL xdiffy(n,dgp,dga,tmpn2)
     
      
! Calculation of tmpn3 = trans(dgp - dga)dm

    CALL vxdiag(n,diag,tmpn2,tmpn3)

         
    pr = vdot(n,tmpn3,tmpn2)
    rrp = vdot(n,tmpn3,dga) 
    CALL xdiffy(n,dg,dga,tmpn4)
    prqr = vdot(n,tmpn3,tmpn4)
    rrq = -vdot(n,tmpn4,d)


! calculation of qr = trans(dg - dga) dm (dg - dga)
       
    CALL vxdiag(n,diag,tmpn4,tmpn3)
    qr = vdot(n,tmpn4,tmpn3)
      
    pq = qr - prqr - prqr + pr
    qqp = pq + prqr + rrq - pr - rrp + alfn
    rrp = rrp - alfv
    rrq = rrq + alfn - alfv


! computation of multipliers lam1 and lam2

    IF (pr > zero .AND. qr > zero) THEN
       tmp1 = rrq/qr
       tmp2 = prqr/qr
       tmp3 = pr - prqr*tmp2
       IF (tmp3 /= zero) THEN
          lam1 = (tmp1*prqr - rrp)/tmp3
          lam2 = -tmp1 - lam1*tmp2
          IF (lam1*(lam1 - one) < zero .AND. &
               lam2*(lam1 + lam2 - one) < zero) GO TO 200
       END IF
    END IF


! Minimum on the boundary

100 continue
    lam1 = zero
    lam2 = zero
    IF (alfn <= alfv) lam2 = one
    IF (qr > zero) lam2 = MIN(one,MAX(zero,-rrq/qr))
    tmp3 = (lam2*qr + rrq+rrq)*lam2
    tmp1 = zero
    IF (alfv >= zero) tmp1 = one
    IF (pr > zero) tmp1 = MIN(one,MAX(zero,-rrp/pr))
    tmp2 = (tmp1*pr + rrp+rrp)*tmp1
    IF (tmp2 < tmp3) THEN
       tmp3 = tmp2
       lam1 = tmp1
       lam2 = zero
    END IF
    
    IF (qqp*(qqp - pq) < zero) THEN
       IF (qr + rrq + rrq - qqp*qqp/pq < tmp3) THEN
          lam1 = qqp/pq
          lam2 = one - lam1
       END IF
    END IF
    
200 CONTINUE
    IF (lam1 == zero .AND. lam2*(lam2 - one) < zero &
         .AND. -rrp - lam2*prqr > zero .AND. pr > zero) &
         lam1 = MIN(one - lam2, (-rrp-lam2*prqr)/pr)


! Computation of the aggregate values
      
    tmp1 = one - lam1 - lam2
    DO i=1,n
       dga(i)=lam1*dgp(i)+lam2*dg(i)+tmp1*dga(i)
    END DO
    
    alfv = lam2*alfn + tmp1*alfv
    
  END SUBROUTINE aggre2
      

!************************************************************************
!*                                                                      *
!*     * SUBROUTINE dgrad *                                             *
!*                                                                      *
!*     Calculation of discrete gradients.                               *
!*                                                                      *
!************************************************************************

  SUBROUTINE dgrad(n,x,sl,g,dg,f,ngrad,nfe,pwt,iterm)  

    USE param, ONLY : zero,half,one,small
    USE obj_fun, ONLY : myf        ! Computation of the value of the objective function.
    USE subpro, ONLY : vdot        ! Dot product
    IMPLICIT NONE

! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
         x, &      ! Current iteration point.
         g         ! Direction of discrete gradient.
    REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: &
         dg        ! Discrete gradient.

! Scalar Arguments
    REAL(KIND=prec), INTENT(IN) :: & 
         f, &      ! Value of the objective at x.
         sl, &     ! Parameter for discrete gradient. (lambda)
         pwt       ! Parameter: pwt = z(sl)*alpha^j*e_j.
    INTEGER, INTENT(IN) :: & 
         n         ! Number of variables
    INTEGER, INTENT(INOUT) :: & 
         ngrad, &  ! Number of discrete gradient calls.
         nfe       ! Number of function evaluations.
    INTEGER, INTENT(OUT) :: & 
         iterm     ! Cause of termination:
                   !   0  - Everything is ok.
                   !  -3  - Failure in function calculations.

! Local arrays
    REAL(KIND=prec), DIMENSION(n) :: x1

! Local Scalars
    REAL(KIND=prec) :: &
         a1, &     ! a1 = argmax{|g_i|}
         r2, &     ! Function value at (x + sl*g)
         r3, &     ! Function value at (x^(j-1))
         r4, &     ! Function value at (x^j)
         dsum, &   ! trans(dg)*g - dg(imax)*g(imax) 
         tmp1      ! Auxiliary scalar.
    INTEGER :: k, &
         imax      ! Index of the largest component of g.
    
! Intrinsic Functions 
    INTRINSIC ABS

! Compute a1 = argmax{|g_i|} and set imax 

    imax = 1
    a1 = zero
    DO k=1,n
       tmp1 = ABS(g(k))
       IF (a1 < tmp1) THEN
          a1 = tmp1
          imax = k
       END IF
    END DO

    DO k=1,n
       x1(k) = x(k)+sl*g(k)
    END DO

    CALL myf(n,x1,r2,iterm) 
    IF (iterm /= 0) RETURN
    nfe = nfe+1

    dsum = zero

    r4 = r2
    DO k=1,n
       IF (k /= imax) THEN
          r3 = r4
          x1(k) = x1(k) + pwt 
          CALL myf(n,x1,r4,iterm) 
          nfe = nfe+1
          IF (iterm /= 0) RETURN
          dg(k) = (r4 - r3)/pwt
          dsum = dsum + dg(k)*g(k)
       END IF
    END DO

    dg(imax) = (r2 - f)/(sl*g(imax)) - dsum/g(imax)

    ngrad = ngrad+1

  END SUBROUTINE dgrad

!************************************************************************
!*                                                                      *
!*     * SUBROUTINE wprint *                                            *
!*                                                                      *
!*     Printout the warning and error messages.                         *
!*                                                                      *
!************************************************************************
      
  SUBROUTINE wprint(iterm,iprint,nout)
    IMPLICIT NONE

! Scalar Arguments
    INTEGER, INTENT(IN) :: &
         iprint, &    ! Printout specification:
                      !  -1  - No printout.
                      !   0  - Only the error messages.
                      !   1  - The final values of the objective
                      !        function.
                      !   2  - The final values of the objective
                      !        function and the most serious
                      !        warning messages.
                      !   3  - The whole final solution. 
                      !   4  - At each iteration values of the
                      !        objective function.
                      !   5  - At each iteration the whole
                      !        solution
         nout, &      ! Auxilary printout specification.
         iterm        ! Cause of termination:
                      !   1  - The problem has been solved with desired accuracy.
                      !   2  - Changes in function values < tolf in mtesf
                      !        subsequent iterations.
                      !   3  - Changes in function value < tolf*small*MAX(|f_k|,|f_(k-1)|,1),
                      !        where small is the smallest positive number such that 
                      !        1.0 + small > 1.0.
                      !   4  - Number of function calls > mfe.
                      !   5  - Number of iterations > mit.
                      !   6  - Time limit exceeded. 
                      !   7  - f < tolb.
                      !   8  - Failure in attaining the demanded accuracy.
                      !  -1  - Two consecutive restarts.
                      !  -2  - Number of restarts > maximum number of restarts.
                      !  -3  - Failure in function calculations (assigned by the user).
                      !  -4  
                      !  -5  - Invalid input parameters.


    IF (iprint >= 0) THEN

! Initial error messages

       IF (iterm == -5) THEN
          IF (nout == 1) WRITE (6,FMT='(1X,''Error: '' &
               ''Number of variables (n) is too small, iterm='',I3)') iterm
          IF (nout == 2) WRITE (6,FMT='(1X,''Error: '' &
               ''The maximum number of stored corrections (mcu) '' &
               ''is too small, iterm='',I3)') iterm
          IF (nout == 4) WRITE (6,FMT='(1X,''Error: '' &
               ''Line search parameter epsl >= 0.25, iterm='',I3)') iterm
          RETURN
       END IF

        
! Warning messages

       IF (iprint >= 2) THEN
          IF (iterm == 0) THEN
             IF (nout == -1) WRITE (6,FMT='(1X,''Warning: '' &
                  ''mc > mcu. Assigned mc = mcu.'')')
             IF (nout == -2) WRITE (6,FMT='(1X,''Warning: '' &
                  ''A line search parameter epsr >= 0.5.'')')
             IF (nout == -3) WRITE (6,FMT='(1X,''Warning: '' &
                  ''A nondescent search direction occured. Restart.'')')
             IF (nout == -4) WRITE (6,FMT='(1X,''Warning: '' &
                  ''Does not converge.'')')
             IF (nout == -5) WRITE (6,FMT='(1X,''Warning: '' &
                  ''tmax < tmin. Restart.'')')
             RETURN
          END IF
         

! Printout the final results
            
          IF (iterm == 6) WRITE (6,FMT='(1X,''Abnormal exit: Time is up.'')')
          IF (iterm == 7) WRITE (6,FMT='(1X,''Abnormal exit: f < tolb.'')')
          IF (iterm == 2) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''Too many steps without significant progress.'')')
          IF (iterm == 3) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''The value of the function does not change.'')')
          IF (iterm == 5) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''Number of iterations > '',I5)') nout
          IF (iterm == 4) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''Number of function evaluations > '',I5)') nout
          IF (iterm == -1) THEN
             IF (nout == -1) THEN
                WRITE (6,FMT='(1X,''Abnormal exit: Two consecutive restarts.'')')
             ELSE
                WRITE (6,FMT='(1X,''Abnormal exit: tmax < tmin in two subsequent iterations.'')')
             END IF
          END IF
          IF (iterm == -2) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''Number of restarts > '',I5''.'')') nout
          IF (iterm == -3) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''Failure in function or subgradient calculations.'')')
          IF (iterm == 8) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''Failure in attaining the demanded accuracy.'')')
       END IF
    END IF
      
  END SUBROUTINE wprint

            
!************************************************************************
!*                                                                      *
!*     * SUBROUTINE rprint *                                            *
!*                                                                      *
!*     Printout the (final) results.                                    *
!*                                                                      *
!************************************************************************
      
  SUBROUTINE rprint(n,noutit,nit,nfe,ngrad,x,f,wk,qk,iterm,iprint)
    IMPLICIT NONE

! Scalar Arguments
    INTEGER, INTENT(IN) :: & 
         n, &         ! Number of variables 
         noutit, &    ! Number of used outer iterations.
         nit, &       ! Number of used iterations.
         nfe, &       ! Number of used function evaluations.
         ngrad, &     ! Number of used discrete gradient calls.
         iprint, &    ! Printout specification:
                      !  -1  - No printout.
                      !   0  - Only the error messages.
                      !   1  - The final values of the objective
                      !        function.
                      !   2  - The final values of the objective
                      !        function and the most serious
                      !        warning messages.
                      !   3  - The whole final solution. 
                      !   4  - At each iteration values of the
                      !        objective function.
                      !   5  - At each iteration the whole
                      !        solution
         iterm        ! Cause of termination:
                      !   1  - The problem has been solved with desired accuracy.
                      !   2  - Changes in function values < tolf in mtesf
                      !        subsequent iterations.
                      !   3  - Changes in function value < tolf*small*MAX(|f_k|,|f_(k-1)|,1),
                      !        where small is the smallest positive number such that 
                      !        1.0 + small > 1.0.
                      !   4  - Number of function calls > mfe.
                      !   5  - Number of iterations > mit.
                      !   6  - Time limit exceeded. 
                      !   7  - f < tolb.
                      !   8  - Failure in attaining the demanded accuracy.
                      !  -1  - Two consecutive restarts.
                      !  -2  - Number of restarts > maximum number of restarts.
                      !  -3  - Failure in function or subgradient calculations 
                      !        (assigned by the user).
                      !  -4  
                      !  -5  - Invalid input parameters.


    REAL(KIND=prec), INTENT(IN) :: &
         f, &         ! Value of the objective function.
         wk, &        ! Value of the first stopping criterion.
         qk           ! Value of the second stopping criterion.

! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
         x            ! Vector of variables
         
! Local Scalars
    INTEGER :: i

! Intermediate results
    
    IF (iterm == 0) THEN
       IF (iprint > 3) WRITE (6,FMT='(1X,''nout='',I3,2X,''nit='',I5,2X, &
            ''nfe='',I7,2X,''ndg='',I7,2X,''f='',D15.8,2X,''wk='',D11.4,2X, &
            ''qk='',D11.4,2X)') noutit,nit,nfe,ngrad,f,wk,qk
       IF (iprint == 5) WRITE (6,FMT='(1X,''x='', &
            5D15.7:/(4X,5D15.7))')(x(i),i=1,n)
       RETURN
    END IF
         

! Final results

    IF (iprint > 0) WRITE (6,FMT='(1X,''nout='',I3,2X,''nit='',I5,2X, &
         ''nfe='',I7,2X,''ndg='',I7,2X,''f='',D15.8,2X,''wk='',D11.4,2X, &
         ''qk='',D11.4,2X,''iterm='',I3)') noutit,nit,nfe,ngrad,f,wk,qk,iterm
    IF (IPRINT .EQ. 3 .OR. IPRINT .EQ. 5) &
         WRITE (6,FMT='(1X,''x='',5D15.7:/(4X,5D15.7))')(x(i),i=1,n)
      
  END SUBROUTINE rprint
      
END MODULE ddgbundle_mod
