MODULE dbundle_mod  ! D-Bundle method

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

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

  PUBLIC :: &
       init_db, &    ! S Initialization for D-bundle subroutine.
       dbundle       ! S D-bundle subroutine for nonsmooth large-scale
                     !   optimization. Contains:
                     !     S restar  Initialization and reinitialization.
                     !     S dobun   Bundle construction.
                     !     F sclpar  Selection of scaling parameter.
  PRIVATE :: &
       tinit, &      ! S Calculation of initial step size. Contains:
                     !     S destep  Stepsize selection using polyhedral 
                     !                 approximation for descent steps.
                     !     S nulstep Stepsize selection using polyhedral 
                     !                 approximation for null steps.
       lls, &        ! S Line search. Contains:
                     !     F qint    Quadratic interpolation.
       aggre1, &     ! S Simplified subgradient aggregation.
       aggre2, &     ! S Subgradient aggregation.
       wprint, &     ! S Printout the error and warning messages.
       rprint        ! S Printout the results.

CONTAINS

!***********************************************************************
!*                                                                     *
!*     * SUBROUTINE init_db *                                          *
!*                                                                     *
!*     Initialization for diagonal bundle subroutine for               *
!*     large-scale unconstrained nonsmooth optimization.               *
!*                                                                     *
!***********************************************************************
      
  SUBROUTINE init_db(iterm)
    USE param, ONLY : zero,half,one,small,large
    USE initializat, ONLY : &
         n, &         ! Number of variables.
         na, &        ! Maximum bundle dimension, na >= 2.
         mcu, &       ! Maximum number of stored corrections, mcu >= 1.
         iprint, &    ! Printout specification (see initializat for more details).
         iscale, &    ! Selection of the scaling (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 (na < 2) THEN
       iterm = -5
       CALL wprint(iterm,iprint,3)
       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 ldgbm_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 first termination criterion.
    IF (told  <= zero) told  = one             ! 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 scaling.

     
  END SUBROUTINE init_db


!***********************************************************************
!*                                                                     *
!*     * SUBROUTINE dbundle *                                          *
!*                                                                     *
!*     Diagonal bundle subroutine for nonsmooth optimization.          *
!*                                                                     *
!***********************************************************************
      
  SUBROUTINE dbundle(f,nit,nfe,nge,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
         na, &        ! Maximum bundle dimension, na >= 2.
         mcu, &       ! Maximum number of stored corrections, mcu >= 1.
         iprint, &    ! Printout specification (see initializat for more details).
         iscale, &    ! Selection of the scaling (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. 
         myg          ! Computation of the subgradient 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) :: & 
         nit,nfe,nge  ! Number of iterations, and function and subgradient 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 or subgradient 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.
         g, &         ! Subgradient of the objective function.
         gp, &        ! Previous subgradient of the ohjective function.
         ga, &        ! Aggregate subgradient.
         s, &         ! Difference of current and previous variables.
         u, &         ! Difference of current and previous subgradients.
         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 subgradients.
    REAL(KIND=prec), DIMENSION(n*na) :: &
         ax,ag        ! Matrix whose columns are bundle points and subgradients.
    REAL(KIND=prec), DIMENSION(na) :: &
         af           ! Vector of bundle values.

! Local Scalars
    REAL(KIND=prec) :: &
         alfn, &      ! Locality measure.
         alfv, &      ! Aggregate locality measure.
         epsr, &      ! Line search parameter.
         dnorm, &     ! Euclidean norm of the direction vector.
         gnorm, &     ! Euclidean norm of the aggregate subgradient.
         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.
    INTEGER :: i,j, &
         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.
         nress, &     ! Number of consecutive restarts in case of tmax < tmin.
         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.
         tmin    = 1.0E-12_prec, &  ! Minimum stepsize.
         lengthd = 1.0E+20_prec     ! Direction vector length.
    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 D-Bundle:'')')
      
     
! Initialization

    nout   = 0
    nit    = 0
    nfe    = 0
    nge    = 0
    ntesf  = 0
    nres   = 1
    ncres  = -1
    nress  = 0
    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
            
     
! Computation of the value and the subgradient of the objective
! function and the search direction for the first iteration

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

    IF (iterm /= 0) GO TO 900
 
    CALL restar(n,mcc,inew,ibun,iters,gp,g,nnk,alfv,alfn,gamma,d,diag,mal,ncres)
          
    CALL dobun(n,na,mal,x,g,f,ax,ag,af,iters,ibun)

     
! Start of the iteration
            
    iteration: DO


! Computational time

       IF (time > 0.0E+00) THEN
          CALL getime(elapstim)
          IF (elapstim-strtim > time) THEN
             iterm = 6
             EXIT iteration
          END IF
       END IF


! Computation of norms

       IF (iters > 0) THEN
          gnorm = vdot(n,g,g)
          dnorm = SQRT(vdot(n,d,d))

          p = vdot(n,g,d)

       ELSE
          gnorm = vdot(n,ga,ga)
          dnorm = SQRT(vdot(n,d,d))

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

   
! Test on descent direction

       IF (p+small*SQRT(gnorm)*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,gp,g,nnk,alfv,alfn,gamma,d,diag,mal,ncres)

             IF (ncres > maxnrs) THEN
                nout = maxnrs
                iterm = -2
                EXIT iteration
             END IF
             
             CALL dobun(n,na,mal,x,g,f,ax,ag,af,iters,ibun)
             
             CYCLE iteration
          END IF
          nout = -1
          iterm = -1
          EXIT iteration
       END IF
       
       
! Stopping criterion

       nit = nit + 1
       pxnorm = xnorm
       xnorm = -p + 2.0_prec*alfv

     
! Tests for termination

       IF(xnorm <= tolg) THEN   ! desired accuracy
          iterm = 1
          EXIT iteration
       END IF


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

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

      
       IF (f <= tolb) THEN  ! too small f
          iterm = 7
          EXIT iteration
       END IF
    
      
       IF (iters == 0) THEN
          IF (ABS(xnorm - pxnorm) <= small) THEN
             neps = neps + 1
          
             IF (neps > maxeps) THEN
                iterm = 8
                EXIT iteration
             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,nit,nfe,nge,x,f,xnorm,half*gnorm+alfv,iterm,iprint)
      
     
!     Preparation of line search

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

       IF (dnorm > zero) tmax = xmax/dnorm
       
       IF (tmax > tmin) THEN
          nress = 0

       ELSE
          nress = nress + 1
          IF (nress == 1) THEN
             CALL wprint(iterm,iprint,-5)
             
             CALL restar(n,mcc,inew,ibun,iters,gp,g,nnk,alfv,alfn,gamma,d,diag,mal,ncres)
             
             IF (ncres > maxnrs) THEN
                nout = maxnrs
                iterm = -2
                EXIT iteration
             END IF
             
             CALL dobun(n,na,mal,x,g,f,ax,ag,af,iters,ibun)
             CYCLE iteration
          END IF
          iterm = -1
          EXIT iteration
       END IF
      

! Initial step size

       CALL tinit(n,na,mal,x,af,ag,ax,ibun,d,f,p,t,tmax,tmin, &
            eta,iters,iterm)

       IF (iterm /= 0) EXIT iteration
     

! Line search with directional derivatives which allows null steps

       theta = one
       IF (dnorm > lengthd) THEN
          theta=lengthd/dnorm
       END IF
      
       CALL lls(n,x,g,d,xo,t,fo,f,p,alfn,tmin,dnorm,xnorm, &
            theta,epsl,epsr,eta,iters,nfe,nge,nnk,iterm)
    
       IF (iterm /= 0) EXIT iteration
       
       IF (tolf2 >= 0) THEN
          IF (ABS(fo-f) <= tolf2*small*MAX(ABS(f),ABS(fo),one) &
               .AND. iters == 1) THEN
             
             iterm = 3
             EXIT iteration
          END IF
       END IF

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

! Bundle updating

       CALL dobun(n,na,mal,x,g,f,ax,ag,af,iters,ibun)
  

! 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,g,gp,u)
             
             CALL aggre1(n,g,gp,ga,u,d,diag,alfn,alfv,gamma)
             
          ELSE
                CALL xdiffy(n,g,gp,u)
             CALL aggre2(n,g,gp,ga,d,diag,alfn,alfv,gamma)
          END IF
          
          CALL copy(n,xo,x)
          f = fo
          
       ELSE
          nnk = 0
          CALL xdiffy(n,g,gp,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(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*g  ! 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*g  ! diag = gamma*I


          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-6_prec) diag(i) = 1.0E-6_prec
                   IF (diag(i) > told) diag(i) = told
                ELSE

                   diag(i) = told

                END IF
             END DO

             CALL vxdiag(n,diag,g,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
                ELSE

                   diag(i) = told

                END IF
             END DO

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


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


          END SELECT

       ELSE  ! Update and direction determination for null steps.
     

          SELECT CASE(iscale)

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

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

! 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*ga
!             diag = gamma


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

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

! 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*ga
!             diag = gamma

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

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

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


!             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,ga,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*ga
!                   diag = gamma
!                ELSE 
!                   CALL vxdiag(n,diag,ga,d)  ! diag = diag
!                   d = -one*d
!                   gamma = gamma_old
!                END IF

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

!             END IF


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


    END DO iteration
      

900 CONTINUE
      

! Printout the final results

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

    CALL wprint(iterm,iprint,nout)
    CALL rprint(n,nit,nfe,nge,x,f,xnorm,half*gnorm+alfv,iterm,iprint)

  CONTAINS


    SUBROUTINE restar(n,mcc,inew,ibun,iters,gp,g,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) :: &
           gp        ! Basic subgradient of the objective function.
      REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
           g         ! Current (auxiliary) subgradient 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,gp,g)
         iters = 1
         nnk = 0
         alfv=zero
         alfn=zero
      END IF

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

    FUNCTION sclpar(iscale,sts,stu,utu) RESULT(spar) ! Calculation of the scaling parameter for D-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.
                     !   4  - Preliminary scaling with STU/UTU. ! tämmöisiä ei varmaan tarvi
                     !   5  - Preliminary scaling with STS/STU. ! tämmöisiä ei varmaan tarvi
                     !   6  - No scaling.      !ok for Dbundle

! 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 < 1.0E-11_prec) spar = 1.0E-11_prec ! no need for correction
!      IF (spar < 1.0E-5_prec) spar = 1.0E-5_prec ! tolg
      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
      

    SUBROUTINE dobun(n,ma,mal,x,g,f,ay,ag,af,iters,ibun)  ! Bundle construction

!      USE subpro, ONLY : copy2 ! given in host
      IMPLICIT NONE

! Array Arguments
      REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
           x, &      ! Vector of variables
           g         ! Subgradient of the objective function.
      REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
           ay, &     ! Matrix whose columns are bundle points 
                     ! (stored in one-dimensional n*ma -array).
           ag, &     ! Matrix whose columns are bundle subgradients 
                     ! (stored in one-dimensional n*ma -array).
           af        ! Vector of values of bundle functions 
                     ! (stored in one-dimensional ma -array).

! Scalar Arguments
      INTEGER, INTENT(IN) :: & 
           n, &      ! Number of variables
           iters, &  ! Null step indicator.
                     !   0  - Null step.
                     !   1  - Serious step.
           ma        ! Maximum size of the bundle.
      INTEGER, INTENT(INOUT) :: &
           ibun, &   ! Index for the circular arrays in bundle updating.
           mal       ! Current size of the bundle.
      REAL(KIND=prec), INTENT(IN) :: & 
           f         ! Value of the objective function.

! Local Scalars
      INTEGER :: i,j
      
      IF (iters == 1) THEN
     
! Serious step

         af(ibun) = f
         i = (ibun-1)*n + 1
         CALL copy2(n,g,ag(i:),x,ay(i:))
         
      ELSE

! Null step

         IF (mal < ma) THEN
            
            af(ibun) = af(mal)
            af(mal) = f
            
            i = mal*n + 1
            CALL copy2(n,ag(i-n:),ag(i:),ay(i-n:),ay(i:))
            CALL copy2(n,g,ag(i-n:),x,ay(i-n:))
            
         ELSE
            i = ibun-1
            IF (i < 1) i = mal
            af(ibun) = af(i)
            af(i) = f
            
            i = (ibun-2)*n + 1
            IF (i < 1) i = (mal-1)*n + 1
            j = (ibun-1)*n + 1
            CALL copy2(n,ag(i:),ag(j:),ay(i:),ay(j:))
            CALL copy2(n,g,ag(i:),x,ay(i:))
         END IF
         
      END IF
      
      mal = mal + 1
      IF (mal > ma) mal = ma
      
      ibun = ibun + 1
      IF (ibun > ma) ibun = 1
      
    END SUBROUTINE dobun
    
  END SUBROUTINE dbundle


!************************************************************************
!*                                                                      *
!*     * SUBROUTINE tinit *                                             *
!*                                                                      *
!*     Initial stepsize selection for limited memory bundle method      *
!*                                                                      *
!************************************************************************

  SUBROUTINE tinit(n,na,mal,x,af,ag,ay,ibun,d,f,p,t,tmax,tmin,eta,iters,iterm)

    USE param, ONLY : zero,half,one,large
!    USE param, ONLY : zero,one  ! these are the one needed in tinit itself 
                                 ! half and large are used in destep and nulstep
    IMPLICIT NONE

! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
         x, &      ! Vector of variables (n array).
         d, &      ! Direction vector (n array).
         ay, &     ! Matrix whose columns are bundle points (stored in one-dimensional n*ma -array).
         ag        ! Matrix whose columns are bundle subgradients (stored in one-dimensional n*ma -array).
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
         af        ! Vector of values of bundle functions (stored in one-dimensional ma -array).

! Scalar Arguments
    REAL(KIND=prec), INTENT(OUT) :: & 
         t         ! Initial stepsize
    REAL(KIND=prec), INTENT(IN) :: & 
         p, &      ! Directional derivative.
         eta, &    ! Distance measure parameter.
         f, &      ! Value of the objective function.
         tmax, &   ! Upper limit for stepsize parameter.
         tmin      ! Lower limit for stepsize parameter.
    INTEGER, INTENT(IN) :: & 
         n, &      ! Number of variables
         na, &     ! Maximum size of the bundle.
         mal, &    ! Current size of the bundle.
         ibun, &   ! Index for the circular arrays in bundle updating.
         iters     ! Null step indicator.
                   !   0  - Null step.
                   !   1  - Serious step.
    INTEGER, INTENT(OUT) :: &
         iterm     ! Cause of termination:
                   !   0  - Everything is ok.
                   !  -6  - Error.


! Intrinsic Functions
    INTRINSIC MAX,MIN

    t = MIN(one,tmax)

    IF (p == zero) RETURN

    IF (iters == 1) THEN
       CALL destep(n,na,mal,x,af,ag,ay,ibun,d,f,p,t,eta,iterm)
    ELSE
       CALL nulstep(n,na,mal,x,af,ag,ay,ibun,d,f,p,t,eta,iterm)
    END IF

    t = MIN(MAX(t,tmin),tmax)

  CONTAINS

    SUBROUTINE destep(n,ma,mal,x,af,ag,ay,ibun,d,f,df,t,eta,iterm)  ! Stepsize selection

!      USE param, ONLY : zero,half,one,large  ! given in host
      IMPLICIT NONE

! Scalar Arguments
      REAL(KIND=prec), INTENT(INOUT) :: & 
           t         ! Initial stepsize
      REAL(KIND=prec), INTENT(IN) :: & 
           df, &     ! Directional derivative.
           eta, &    ! Distance measure parameter.
           f         ! Value of the objective function.
      INTEGER, INTENT(IN) :: & 
           n, &      ! Number of variables
           ma, &     ! Maximum size of the bundle.
           mal, &    ! Current size of the bundle.
           ibun      ! Index for the circular arrays in bundle updating.
      INTEGER, INTENT(OUT) :: &
           iterm     ! Cause of termination:
                     !   0  - Everything is ok.
                     !  -6  - Error.

! Array Arguments
      REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
           x, &      ! Vector of variables (n array).
           d, &      ! Direction vector (n array).
           ay, &     ! Matrix whose columns are bundle points (stored in one-dimensional n*ma -array).
           ag, &     ! Matrix whose columns are bundle subgradients (stored in one-dimensional n*ma -array).
           af        ! Vector of values of bundle functions (stored in one-dimensional ma -array).

! Local Arrays
      REAL(KIND=prec), DIMENSION(2*ma) :: &
           tmparray  ! Auxiliary array.

! Local Scalars
      REAL(KIND=prec) :: alf,alfl,alfr,bet,betl,betr,dx,q,r,w
      INTEGER :: i,j,jn,k,l,lq,ib

! Intrinsic Functions
      INTRINSIC ABS,REAL,MAX,MIN,SQRT

      iterm = 0
      alfl = zero
      betl = zero
      
      w = df*t* (one - t*half)
      
     
! Initial choice of possibly active lines
      
      k = 0
      l = -1
      jn = (ibun-1)*n
      betr = - large
      DO j=1,mal-1
         ib = ibun - 1 + j
         IF (ib > mal) ib = ib - mal
         IF (jn >= mal*n) jn = jn - mal*n
         r = zero
         bet = zero
         alfl = af(ib) - f
         DO i=1,n
            dx = x(i) - ay(jn+i)
            q = ag(jn+i)
            r = r + dx*dx
            alfl = alfl + dx*q
            bet = bet + d(i)*q
         END DO
         alf = MAX(ABS(alfl),eta*r)
         r = one - bet/df
         IF (r*r + (alf+alf)/df > 1.0E-6_prec) THEN
            k = k + 1
            tmparray(k) = alf
            tmparray(ma+k) = bet
            r = t*bet - alf
            IF (r > w) THEN
               w = r
               l = k
            END IF
         END IF
         
         betr = MAX(betr,bet-alf)
         jn = jn + n
      END DO
      lq = -1
      IF (betr <= df*half) RETURN
      lq = 1
      IF (l <= 0) RETURN
      betr = tmparray(ma+l)
      IF (betr <= zero) THEN
         IF (t < one .OR. betr == zero) RETURN
         lq = 2
      END IF
      
      alfr = tmparray(l)


! Iteration loop

      ds_iteration: DO 
         
         IF (lq >= 1) THEN
            q = one - betr/df
            r = q + SQRT(q*q + (alfr+alfr)/df)
            IF (betr >= zero) r = - (alfr+alfr)/ (df*r)
            r = MIN(1.95_prec,MAX(zero,r))
         ELSE
            IF (ABS(betr-betl)+ABS(alfr-alfl)< -1.0E-4_prec*df) RETURN
            IF (betr-betl  == zero) THEN
               iterm = -6
               RETURN
            END IF
            r = (alfr-alfl)/ (betr-betl)
         END IF
       
         IF (ABS(t-r) < 1.0E-4_prec) RETURN
         t = r
         tmparray(l) = - one
         w = t*betr - alfr
         l = -1
         DO j = 1,k
            alf = tmparray(j)
            IF (alf < zero) EXIT
            bet = tmparray(ma+j)
            r = t*bet - alf
            IF (r > w) THEN
               w = r
               l = j
            END IF
         END DO
         IF (l < 0) RETURN
       
         bet = tmparray(ma+l)
         IF (bet == zero) RETURN
     

!     New interval selection

         alf = tmparray(l)
         IF (bet < zero) THEN
            IF (lq == 2) THEN
               alfr = alf
               betr = bet
               
            ELSE
               alfl = alf
               betl = bet
               lq = 0
            END IF
       
         ELSE
            IF (lq == 2) THEN
               alfl = alfr
               betl = betr
               lq = 0
            END IF
       
            alfr = alf
            betr = bet
         END IF
       
      END DO ds_iteration
      
    END SUBROUTINE destep

     
    SUBROUTINE nulstep(n,ma,mal,x,af,ag,ay,ibun,d,f,df,t,eta,iterm)  ! Stepsize selection

!      USE param, ONLY : zero,half,one,large  ! given in host
      IMPLICIT NONE

! Scalar Arguments
      REAL(KIND=prec), INTENT(INOUT) :: & 
           t         ! Initial stepsize
      REAL(KIND=prec), INTENT(IN) :: & 
           df, &     ! Directional derivative.
           eta, &    ! Distance measure parameter.
           f         ! Value of the objective function.
      INTEGER, INTENT(IN) :: & 
           n, &      ! Number of variables
           ma, &     ! Maximum size of the bundle.
           mal, &    ! Current size of the bundle.
           ibun      ! Index for the circular arrays in bundle updating.
      INTEGER, INTENT(OUT) :: &
           iterm     ! Cause of termination:
                     !   0  - Everything is ok.
                     !  -6  - Error.

! Array Arguments
      REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
           x, &      ! Vector of variables (n array).
           d, &      ! Direction vector (n array).
           ay, &     ! Matrix whose columns are bundle points (stored in one-dimensional n*ma -array).
           ag, &     ! Matrix whose columns are bundle subgradients (stored in one-dimensional n*ma -array).
           af        ! Vector of values of bundle functions (stored in one-dimensional 4*ma -array).

! Local Arrays
      REAL(KIND=prec), DIMENSION(2*ma) :: &
           tmparray  ! Auxiliary array.

! Local Scalars
      REAL(KIND=prec) :: alf,alfl,alfr,bet,betl,betr,dx,q,r,w
      INTEGER :: i,j,jn,k,l,ib

! Intrinsic Functions
      INTRINSIC ABS,REAL,MAX,MIN,SQRT

      iterm = 0
      w = df*t
     
! Initial choice of possibly active parabolas

      k = 0
      l = -1
      jn = (ibun-1)*n
      betr = - large
      DO j = 1,mal - 1
         ib = ibun - 1 + j
         IF (ib > mal) ib = ib - mal
         IF (jn >= mal*n) jn = jn - mal*n
         bet = zero
         r = zero
         alfl = af(ib) - f
         DO i = 1,n
            dx = x(i) - ay(jn+i)
            r = r + dx*dx
            q = ag(jn+i)
            alfl = alfl + dx*q
            bet = bet + d(i)*q
         END DO
         alf = MAX(ABS(alfl),eta*r)
         betr = MAX(betr,bet-alf)
         IF (alf < bet-df) THEN
            k = k + 1
            r = t*bet - alf
            tmparray(k) = alf
            tmparray(ma+k) = bet
            IF (r > w) THEN
               w = r
               l = k
            END IF
         END IF
         jn = jn + n
      END DO
      IF (l <= 0) RETURN
      betr = tmparray(ma+l)
      alfr = tmparray(l)
      alf = alfr
      bet = betr
      alfl = zero
      betl = df
    
     
!     Iteration loop
    
      ns_iteration: DO

         w = bet/df
         IF (ABS(betr-betl)+ABS(alfr-alfl)< -1.0E-4_prec*df) RETURN
         IF (betr-betl  == zero) THEN
            iterm = -6
            RETURN
         END IF
         r = (alfr-alfl)/ (betr-betl)
         IF (ABS(t-w) < ABS(t-r)) r = w
         q = t
         t = r
         IF (ABS(t-q) < 1.0E-3_prec) RETURN
         tmparray(l) = - one
         w = t*bet - alf
         l = -1
         DO j=1,k
            alf = tmparray(j)
            IF (alf < zero) EXIT
            bet = tmparray(ma+j)
            r = t*bet - alf
            IF (r > w) THEN
               w = r
               l = j
            END IF
         END DO

         IF (l <= 0) RETURN
         bet = tmparray(ma+l)
         q = bet - t*df
         IF (Q == zero) RETURN

     
!     New interval selection

         alf = tmparray(l)
         IF (q < zero) THEN
            alfl = alf
            betl = bet
         ELSE
            alfr = alf
            betr = bet
         END IF

      END DO ns_iteration

    END SUBROUTINE nulstep

  END SUBROUTINE tinit
      

!************************************************************************
!*                                                                      *
!*     * SUBROUTINE lls *                                               *
!*                                                                      *
!*     Special line search for limited memory bundle method             *
!*                                                                      *
!************************************************************************

  SUBROUTINE lls(n,x,g,d,xo,t,fo,f,p,alfn,tmin,dnorm,wk,theta,epsl,epsr,&
       eta,iters,nfe,nge,nnk,iterm)

    USE param, ONLY : zero,half,one
    USE obj_fun, ONLY : &
         myf,myg   ! Computation of the value and the subgradient 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) :: &
         g         ! Subgradient of the objective function.

! Scalar Arguments
    REAL(KIND=prec), INTENT(INOUT) :: & 
         epsl,&    ! Linesearch parameter.
         epsr,&    ! Linesearch parameter.
         t, &      ! Stepsize
         p         ! Directional derivative.
    REAL(KIND=prec), INTENT(IN) :: & 
         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.
         nge       ! Number of subgradient evaluations.
    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 Scalars
    INTEGER :: nin ! Number of interpolations.
    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.
         epsa,epst,epslk,epsrk, & ! Line search parameters.
         thdnorm,epsawk,epstwk,epslwk,epsrwk ! Auxiliary scalars.

! 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

    tl = zero
    tu = t
    fl = 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 evaluation at a new point

    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 .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,wk,one-half/(one-epst))
          ELSE
             t = half*(tu+tl)
          END IF
          CYCLE lls_iteration
       END IF

       CALL myg(n,x,g,iterm)
       nge = nge + 1
       IF (iterm /= 0) RETURN

       p = theta*vdot(n,g,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 = qint(tu,fl,fu,wk,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

  END SUBROUTINE lls


!************************************************************************
!*                                                                      *
!*     * SUBROUTINE update *                                            *
!*                                                                      *
!*     Matrix update for D-bundle                                       *
!*                                                                      *
!************************************************************************
    
  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,g,gp,ga,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.
         g, &      ! Current (auxiliary) subgradient of the objective function.
         gp, &     ! Previous subgradient of the objective function.
         u         ! Difference of trial and aggregate gradients.
    REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: &
         ga        ! 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
       ga(i)=lam*g(i) + p*gp(i)
    END DO
      
    alfv = lam*alfn
      
  END SUBROUTINE aggre1


!************************************************************************
!*
!*     * SUBROUTINE aggre2 *
!*
!*     Computation of aggregate values after consequtive null steps.
!*
!************************************************************************
      
  SUBROUTINE aggre2(n,g,gp,ga,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.
         g, &      ! Current (auxiliary) subgradient of the objective function.
         gp        ! Previous subgradient of the objective function.
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
         ga        ! 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(gp-ga) dm (gp-ga), where dm
                   ! presents the L-SR1- approximation of Hessian.
         rrp, &    ! rrp = trans(gp-ga) dm ga - alfv.
         prqr, &   ! prqr = trans(gp-ga) dm (g-ga).
         rrq, &    ! rrq = trans(g-ga) dm ga - alfv + alfn.
         qr, &     ! qr = trans(g-ga) dm (g-ga).
         pq, &     ! pq = trans(g-gp) dm (g-gp).
         qqp, &    ! qqp = trans(g-gp) dm g + 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,gp,ga,tmpn2)
     
      
! Calculation of tmpn3 = trans(gp - ga)dm

    CALL vxdiag(n,diag,tmpn2,tmpn3)

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


! calculation of qr = trans(g - ga) dm (g - ga)

       
    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
       ga(i)=lam1*gp(i)+lam2*g(i)+tmp1*ga(i)
    END DO
    
    alfv = lam2*alfn + tmp1*alfv
    
  END SUBROUTINE aggre2
      

!************************************************************************
!*
!*     * 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 or subgradient 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 == 3) WRITE (6,FMT='(1X,''Error: '' &
               ''The size of the bundle (na) 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,nit,nfe,nge,x,f,wk,qk,iterm,iprint)
    IMPLICIT NONE

! Scalar Arguments
    INTEGER, INTENT(IN) :: & 
         n, &         ! Number of variables 
         nit, &       ! Number of used iterations.
         nfe, &       ! Number of used function evaluations.
         nge, &       ! Number of used subgradient evaluations.
         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,''nit='',I5,2X, &
            ''nfe='',I5,2X,''nge='',I5,2X,''f='',D15.8,2X,''wk='',D11.4,2X, &
            ''qk='',D11.4,2X)') nit,nfe,nge,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,''nit='',I5,2X, &
         ''nfe='',I5,2X,''nge='',I5,2X,''f='',D15.8,2X,''wk='',D11.4,2X, &
         ''qk='',D11.4,2X,''iterm='',I3)') nit,nfe,nge,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 dbundle_mod
