!*************************************************************************
!*                                                                       *
!*     LMBM     - Limited Memory Bundle Method for nonsmooth             *
!*                nonconvex optimization                                 *
!*                                                                       *
!*                                                                       *
!*     by Napsu Karmitsa 2004 (last modified 01.09.2020)                 *
!*                                                                       *
!*                                                                       *
!*     This file is specially modified to solve NN regression problems.  *
!*                                                                       *
!*************************************************************************
!*
!*     Modules included:
!*
!*     lmbm_mod            ! Limited Memory Bundle Method
!*
MODULE lmbm_mod  ! Limited memory bundle method

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

    ! MODULE lmbm_mod includes the following subroutines (S) and functions (F).
    PUBLIC :: &
        optim2        ! S Initializing LMBM subroutine for NN.
    PRIVATE :: &
        lmbm, &       ! S Limited memory bundle subroutine for nonsmooth
                      !   large-scale optimization. Contains:
                      !     S restar  Initialization and reinitialization.
                      !     S dobun   Bundle construction.
        indic1, &     ! S Initialization of indices.
        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 Weak Wolfe line search. Contains:
                      !     F qint    Quadratic interpolation.
        nmlls, &      ! S Nonmonotonic Weak Wolfe line search. Contains:
                      !     F qint    Quadratic interpolation.
        armijo, &     ! S Armijo line search.
        nmarmijo, &   ! S Nonmonotone Armijo line search.
        dlbfgs, &     ! S Computing the search direction by limited memory
                      !   BFGS update. Contains:
                      !     F sclpar  Selection of scaling parameter.
        dlsr1, &      ! S Computing the search direction by limited
                      !   memory SR1 update.
        agbfgs, &     ! S Simplified subgradient aggregation.
        aggsr1, &     ! S Subgradient aggregation.
        agskip, &     ! S Subgradient aggregation using BFGS update.
        wprint, &     ! S Printout the error and warning messages.
        rprint        ! S Printout the results.

CONTAINS

    !SUBROUTINE optim2(x0,x1,fvalue)
    SUBROUTINE optim2(fvalue)

        USE initnn, ONLY : &
            nH, &          ! Current number of hidden nodes.
            maxdim, &      ! Maximum number of variables.
            nft, &         ! Number of input features:
            !nc => nH, &          ! Current number of hidden nodes.
            ns             ! Switch for auxiliary and real NN problem.
            !m              ! Number of variables for function calculations. ! tarvitaanko täällä?
        USE initlmbm, ONLY : &
            n, &           ! Number of variables.
            x, &           ! Vector of variables.
            mcu, &         ! Maximum number of stored corrections.
            mcinit, &      ! Initial maximum number of stored corrections.
            epsl, &        ! Line search parameter.
            defaults, &    ! S  Default values for parameters.
            init_lmbmpar   ! S  Further initialization of parameters.

        REAL(KIND=prec), INTENT(OUT) ::  &
            fvalue

        INTEGER, DIMENSION(4) :: &
            iout           ! Output integer parameters.
                           !   iout(1)   Number of used iterations.
                           !   iout(2)   Number of used function evaluations.
                           !   iout(3)   Number of used subgradient evaluations
                           !   iout(4)   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.
        INTEGER :: i,mc

        ! number of variables
        !n = nH*(nft+2)+1

        mc = mcinit     ! Initial maximum number of stored corrections

        CALL defaults()    !
        CALL init_lmbmpar()    !

        IF (n <= 0) PRINT*,'n<0'
        IF (n > maxdim) PRINT*,'n is greater than maxdim.'
        IF (epsl >= 0.25_prec) PRINT*,'epsl >= 0.25'
        IF (mcu <= 0) PRINT*,'mcu <= 0'

        CALL lmbm(mc,fvalue,iout(1),iout(2),iout(3),iout(4))

        RETURN

    END SUBROUTINE optim2


    !***********************************************************************
    !*                                                                     *
    !*     * SUBROUTINE lmbm *                                             *
    !*                                                                     *
    !*     Limited memory bundle subroutine for nonsmooth optimization.    *
    !*                                                                     *
    !***********************************************************************
      
    SUBROUTINE lmbm(mc,f,nit,nfe,nge,iterm)
        USE param, ONLY : small,large,zero,half,one
        USE initlmbm, ONLY : &
            n, &         ! Number of variables.
            x, &         ! Vector of variables.
            forig, &     ! f at x0.
            na, &        ! Maximum bundle dimension, na >= 2.
            mcu, &       ! Upper limit for maximum number of stored corrections, mcu >= 3.
            inma, &      ! Selection of line search method:
                         !   inma = 0, Armijo line search,
                         !   inma = 1, nonmonotone Armijo line search.
                         !   inma = 2, weak Wolfe line search.
            mnma, &      ! Maximum number of function values used in nonmonotone line search.
            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.
            tolg2, &     ! Tolerance for the second termination criterion.
            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.
        USE objfun, ONLY : &
            myf, &       ! Computation of the value of the objective function.
            myg          ! Computation of the subgradient of the objective function.
        USE subpro, ONLY : &
            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(INOUT) :: &
            mc           ! Maximum number of stored corrections.
        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.
                          !  -1  - Two consecutive restarts.
                          !  -2  - Number of restarts > maximum number of restarts.
                          !  -3  - Failure in function or subgradient calculations
                          !        (assigned by the user).
                          !  -4  - Failure in attaining the demanded accuracy.
                          !  -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.
            tmpn1        ! Auxiliary array.
        REAL(KIND=prec), DIMENSION(n*(mcu+1)) :: &
            sm,um        ! Matrises whose columns are stored differences of
                          ! variables (sm) and subgradients (um).
        REAL(KIND=prec), DIMENSION((mcu+2)*(mcu+1)/2) :: &
            rm, &        ! pper triangular matrix stored columnwise in the
                         ! one-dimensional array.
            umtum        ! Matrix whose columns are stored subgradient differences.
        REAL(KIND=prec), DIMENSION(mcu+1) :: &
            cm, &        ! Diagonal matrix.
            smtgp, &     ! smtgp = trans(sm)*gp.
            umtgp, &     ! umtgp = trans(um)*gp.
            tmpmc1, &    ! Auxiliary array.
            tmpmc2       ! Auxiliary array.
        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.
        REAL(KIND=prec), DIMENSION(mnma) :: &
            fold         ! Old function 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.
            gamma        ! Scaling parameter.
        INTEGER :: i, &
            mcinit, &    ! Initial maximum number of stored corrections.
            mcc, &       ! Current number of stored corrections.
            inew, &      ! Index for the circular arrays.
            ibfgs, &     ! Index of the type of BFGS update.
            isr1, &      ! Index of the type of SR1 update.
            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.
            inewnma, &   ! Index for the circular arrays containing function values.
            ic, &        ! Correction indicator.
            icn, &       ! Correction indicator for null steps.
            iflag, &     ! Index for adaptive version.
            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,DOT_PRODUCT

       
        ! 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.
            rho     = 1.0E-12_prec     ! Correction parameter.
        INTEGER, PARAMETER :: &
            maxeps = 20, &             ! Maximum number of consecutive equal stopping criterions.
            maxnrs = 2000              ! Maximum number of restarts.

        ! Initialization

        inewnma = 1
        nout   = 0
        nit    = 0
        nfe    = 0
        nge    = 0
        ntesf  = 0
        nres   = 1
        ncres  = -1
        nress  = 0
        neps   = 0
        iterm  = 0
        iters  = 1
        nnk    = 0
        isr1   = 0
        alfn   = zero
        alfv   = zero
        mcinit = mc
      
        tmax   = xmax
        xnorm  = large
        fold   = -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)
        PRINT*,'f original =',f
        forig = f
        CALL myg(n,x,g,iterm)
        nfe = nfe + 1
        nge = nge + 1

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

     
        ! Start of the iteration
            
        iteration: DO



            ! Computation of norms

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

                p = DOT_PRODUCT(g,d)
          
            ELSE
                gnorm = DOT_PRODUCT(ga,ga)
                dnorm = SQRT(DOT_PRODUCT(d,d))

                p = DOT_PRODUCT(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,mc,mcc,mcinit,inew,ibun,ibfgs,iters,gp,g,nnk, &
                        alfv,alfn,gamma,d,ic,icn,mal,ncres,iflag)
                    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 <= 1.0E+03_prec*tolg .AND. &
                (mcc > 0 .OR. ibfgs == 2)) THEN
          
                IF(half*gnorm + alfv <= tolg2 .AND. &
                    xnorm <= tolg) THEN
          
                    iterm = 1
                    EXIT iteration
                END IF
       
                IF (mc < mcu .AND. iflag == 0) THEN
                    mc = mc+1
                    iflag = 1
                END IF
            END IF


            IF (nfe >= mfe) THEN
                nout = mfe
                iterm = 4
                EXIT iteration
            END IF

      
            IF (nit >= mit) THEN
                nout = mit
                iterm = 5
                EXIT iteration
            END IF

      
            IF (f <= tolb) THEN
                iterm = 7
                EXIT iteration
            END IF
    
      
            IF (iters == 0) THEN
                IF (ABS(xnorm - pxnorm) <= small) THEN
                    neps = neps + 1
          
                    IF (neps > maxeps) THEN
                        iterm = -4
                        EXIT iteration
                    END IF

                ELSE
                    neps = 0
                END IF

            ELSE
                neps = 0
            END IF


            ! Correction
    
            IF (-p < rho*gnorm .OR. icn == 1) THEN

                xnorm = xnorm + rho*gnorm
                dnorm = SQRT(dnorm*dnorm-2.0_prec*rho*p+rho*rho*gnorm)
         
                IF (iters > 0) THEN
                    DO i=1,n
                        d(i) = d(i)-rho*g(i)
                    END DO

                ELSE
                    DO i=1,n
                        d(i) = d(i)-rho*ga(i)
                    END DO
                    icn = 1
                END IF

                ic = 1
            
            ELSE
                ic = 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,mc,mcc,mcinit,inew,ibun,ibfgs,iters,gp,g,nnk, &
                        alfv,alfn,gamma,d,ic,icn,mal,ncres,iflag)
             
                    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
     
            IF (inma == 0) THEN
                CALL armijo(x,g,d,xo,fo,f,t,p,alfn,xnorm,epsr,iters,nfe,nge,iterm)

            ELSE IF (inma == 1) THEN
                CALL nmarmijo(x,g,d,xo,fo,f,fold,t,p,alfn,xnorm,epsr,iters,nfe,nge,inewnma,iterm)

            ELSE IF (inma == 2) THEN! Line search with directional derivatives which allows null steps
                 ! With this the global convergence can be guaranteed.
                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)
            ELSE ! Nonmonotone line search with directional derivatives which allows null steps
                 ! With this the global convergence can be guaranteed.
                theta = one
                IF (dnorm > lengthd) THEN
                    theta=lengthd/dnorm
                END IF

                CALL nmlls(n,x,g,d,xo,t,fo,f,fold,p,alfn,tmin,dnorm,xnorm, &
                    theta,epsl,epsr,eta,iters,inewnma,nfe,nge,nnk,iterm)

            END IF

            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 copy(n,gp,tmpn1)
                    CALL xdiffy(n,g,gp,u)
             
                    CALL agbfgs(n,mc,mcc,inew,ibfgs,iflag,g,gp,ga,u,d,sm,um, &
                        rm,cm,umtum,alfn,alfv,gamma,ic,rho)
             
                ELSE
                    CALL copy(n,ga,tmpn1)
                    CALL aggsr1(n,mc,mcc,inew,iflag,g,gp,ga,d,alfn,alfv &
                        ,umtum,rm,gamma,smtgp,umtgp,tmpmc1,tmpmc2,sm &
                        ,um,icn,rho)
                    CALL xdiffy(n,g,gp,u)
                END IF
          
                CALL copy(n,xo,x)
                f = fo
          
            ELSE
                IF (nnk /= 0) THEN
                    CALL copy(n,ga,tmpn1)
                ELSE
                    CALL copy(n,gp,tmpn1)
                END IF
                nnk = 0
                CALL xdiffy(n,g,gp,u)
            END IF

     
            ! Serious step initialization

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

      
            ! Direction finding
    
            IF (iters > 0) THEN
         
     
                ! BFGS update and direction determination

                CALL dlbfgs(n,mc,mcc,inew,ibfgs,iflag,d,g,gp,s,u,sm,um,rm, &
                    umtum,cm,smtgp,umtgp,gamma,tmpn1,iscale)

            ELSE

                ! SR1 update and direction determination
             
                CALL dlsr1(n,mc,mcc,inew,isr1,iflag,d,gp,ga,s,u,sm,um,rm, &
                    umtum,cm,smtgp,umtgp,gamma,tmpmc1,tmpmc2,tmpn1,nnk,iprint)
                ibfgs=0
            END IF
    
        END DO iteration
      

900 CONTINUE
      

    ! Printout the final results

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

CONTAINS

    SUBROUTINE restar(n,mc,mcc,mcinit,inew,ibun,ibfgs,iters,gp,g,nnk, &
        alfv,alfn,gamma,d,ic,icn,mal,ncres,iflag)  ! Initialization
      
        !      USE param, ONLY : zero,one  ! given in host
        !      USE lmbm_sub, 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.

        ! Scalar Arguments
        INTEGER, INTENT(IN) :: &
            n, &      ! Number of variables
            mcinit    ! Initial maximum number of stored corrections.
        INTEGER, INTENT(OUT) :: &
            mc, &     ! Current maximum number of stored corrections.
            mcc, &    ! Current number of stored corrections.
            inew, &   ! Index for the circular arrays.
            ibun, &   ! Index for the circular arrays in bundle updating.
            ibfgs, &  ! Index of the type of BFGS update.
            nnk, &    ! Consecutive null steps counter.
            ic, &     ! Correction indicator.
            icn, &    ! Correction indicator for null steps.
            mal, &    ! Current size of the bundle.
            iflag     ! Index for adaptive version.
        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
        mc    = mcinit
        mcc   = 0
        inew  = 1
        ibun  = 1
        ibfgs = 0
        ic    = 0
        icn   = 0
        mal   = 0
        ncres = ncres + 1
        iflag = 0
        
        IF (iters == 0) THEN
            CALL copy(n,gp,g)
            iters = 1
            nnk = 0
            alfv=zero
            alfn=zero
        END IF

        gamma = one
        d=-g
      
    END SUBROUTINE restar
    
      
    SUBROUTINE dobun(n,ma,mal,x,g,f,ay,ag,af,iters,ibun)  ! Bundle construction

        !      USE lmbm_sub, 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 lmbm

!************************************************************************
!*                                                                      *
!*     * 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 objfun, 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.
    USE initlmbm, ONLY : &
        maxint => maxnin ! Maximum number of interpolations.

    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.

    ! 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*DOT_PRODUCT(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 nmlls *                                             *
!*                                                                      *
!*     Nonmonotonic weak Wolfe-type line search                         *
!*                                                                      *
!************************************************************************

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

    USE param, ONLY : zero,half,one
    USE objfun, 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.
    USE initlmbm, ONLY : &
        mnma, &    ! Maximum number of function values used in line search.
        maxint => maxnin   ! Maximum number of interpolations.

    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.
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
        fold           ! Old function values.

    ! 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) :: &
        inewnma, &! index for array.
        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.
        ftmp, &   ! Maximum function value from mnma last iterations.
        epsa,epst,epslk,epsrk, & ! Line search parameters.
        thdnorm,epsawk,epstwk,epslwk,epsrwk ! Auxiliary scalars.

    ! Intrinsic Functions
    INTRINSIC ABS,MAX,MAXVAL,DOT_PRODUCT



    ! 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

    ! Updating circular array fold

    IF (iters==1) THEN
        fold(inewnma) = fo
        inewnma = inewnma + 1
        IF (inewnma > mnma) inewnma = 1
    END IF
    ftmp = MAXVAL(fold)


    ! 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 <= ftmp - t*epstwk) THEN
            tl = t
            fl = f
        ELSE
            tu = t
            fu = f
        END IF


        ! Additional interpolation

        IF (f > ftmp .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*DOT_PRODUCT(g,d)
        alfn = MAX(ABS(fo-f+p*t),eta*(t*thdnorm)**2)


        ! Serious step

        IF (f <= ftmp - 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 nmlls



!************************************************************************
!*                                                                      *
!*     * SUBROUTINE armijo *                                            *
!*                                                                      *
!*     Armijo line search for limited memory bundle method              *
!*                                                                      *
!************************************************************************


    ! Armijo line search
    SUBROUTINE armijo(x,g,d,xo,fo,f,tin,p,alfn,wk,epsr,iters,nfe,nge,iterm)
        USE param, ONLY : two,one,zero  ! given in host
        USE initlmbm, ONLY : &
            n, &         ! Number of variables.         ! given in host
            epsl, &      ! Descent parameter.           ! given in host
            eta, &       ! Distance measure parameter.  ! given in host
            maxnin         ! Maximum number of interpolations.
        USE objfun, ONLY : &
            myf, &       ! Computation of the value f = f1 - f2.
            myg          ! Computation of the gradient of the first DC-component,
                         ! given in host
        IMPLICIT NONE

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

        ! Scalar Arguments
        REAL(KIND=prec), INTENT(INOUT) :: &
            epsr           ! Linesearch parameter, used only in null step test to print an error message
        REAL(KIND=prec), INTENT(IN) :: &
            fo, &          ! Previous value of the objective function.
            tin, &         ! Initial stepsize.
            wk             ! Stopping parameter.
        REAL(KIND=prec), INTENT(OUT) :: &
            f, &           ! Value of the objective function.
            p, &           ! Directional derivative.
            alfn           ! Linearization error.
        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 arrays
        REAL(KIND=prec), DIMENSION(n) :: &
            xu             ! trial point

        ! Local Scalars
        REAL(KIND=prec) :: &
            fu,fu1,fu2, &  ! Auxiliary values of the objective functions
            a_orig, &      ! Armijo parameter.
            inv_a_orig, &  ! Inverse of Armijo parameter.
            alpha, &       ! Armijo step size.
            dnorm, &       ! |d|
            epslwk,epsrwk  ! Auxiliary scalars.
        INTEGER :: i,nin   ! Number of interpolations.

        ! Intrinsic Functions
        INTRINSIC ABS,MAX,DOT_PRODUCT

        ! Initialization
        nin = 0

        epslwk   = epsl*wk
        epsrwk   = epsr*wk

        alpha  = two*tin
        a_orig = two
        inv_a_orig = one/a_orig

        dnorm = DOT_PRODUCT(d,d)

        ! Function evaluation at a new point

        x = xo + d

        CALL myf(n,x,f,iterm)
        nfe = nfe + 1
        IF (iterm /= 0) THEN
            f=fo
            RETURN
        END IF



        ! Serious step

        iters = 1

        IF (f <= fo - epslwk) THEN ! descent condition

            xu = xo + alpha*d
            CALL myf(n,xu,fu,iterm)
            nfe = nfe + 1
            IF (iterm /= 0) RETURN

            IF(fu > fo - alpha*epslwk .OR. maxnin==0) THEN ! Serious step with t=tin

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

                p = DOT_PRODUCT(g,d)
                alfn = MAX(ABS(fo-f+p*tin),eta*(tin*dnorm)**2)

                RETURN
            END IF

            inct: DO ! Increase t

                nin=nin+1
                f=fu
                x=xu
                alpha = alpha * a_orig
                xu = xo + alpha*d

                CALL myf(n,xu,fu,iterm)
                nfe = nfe + 1
                IF (iterm /= 0) RETURN

                IF (fu > fo - alpha*epslwk .OR. nin >= maxnin) THEN ! Serious step

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

                    p = alpha * inv_a_orig * DOT_PRODUCT(g,d)
                    alfn = MAX(ABS(fo-f+p),eta*(alpha * inv_a_orig*dnorm)**2)

                    RETURN
                END IF
            END DO inct

        ELSE

            alpha = tin

            dect: DO ! Decrease t
                nin=nin+1
                IF (nin > maxnin) EXIT dect

                alpha = alpha * inv_a_orig !
                x = xo + alpha * d
                CALL myf(n,x,f,iterm)
                nfe = nfe + 1
                IF (iterm /= 0) RETURN

                IF (f <= fo - alpha*epslwk) THEN ! Serious step

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

                    p = alpha * DOT_PRODUCT(g,d)
                    alfn = MAX(ABS(fo-f+p),eta*(alpha*dnorm)**2)

                    RETURN
                END IF
            END DO dect

        END IF


        iters = 0 ! Null step

        ! original d
        !       x = xo + d
        !      CALL myf(n,x,f,iterm)
        !      nfe = nfe + 1


        CALL myg(n,x,g,iterm) ! x = xo + alpha*d
        nge = nge + 1
        IF (iterm /= 0) RETURN

        !       p = DOT_PRODUCT(g,d) ! original d
        p = alpha * DOT_PRODUCT(g,d)
        alfn = MAX(ABS(fo-f+p),eta*(alpha*dnorm)**2)


    END SUBROUTINE armijo


!************************************************************************
!*                                                                      *
!*     * SUBROUTINE nmarmijo *                                          *
!*                                                                      *
!*     Nonmonotone Armijo line search for limited memory bundle method  *
!*                                                                      *
!************************************************************************

    ! Nonmonotone Armijo line search.
    SUBROUTINE nmarmijo(x,g,d,xo,fo,f,fold,tin,p,alfn,wk,epsr,iters,nfe,nge,inewnma,iterm)
        USE param, ONLY : two,one,zero
        USE initlmbm, ONLY : &
            n, &           ! Number of variables.
            epsl, &        ! Descent parameter.
            eta, &         ! Distance measure parameter.
            mnma, &        ! Maximum number of function values used in line search.
            maxnin         ! Maximum number of interpolations.
        USE objfun, ONLY : &
            myf, &         ! Computation of the value f = f1 - f2.
            myg            ! Computation of the subgradient,
            !              ! given in host
        IMPLICIT NONE

        ! Array Arguments
        REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
            d, &           ! Direction vector.
            xo             ! Previous vector of variables.
        REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: &
            x, &           ! Vector of variables.
            g              ! Subgradient of the objective function.
        REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
            fold           ! Old function values.

        ! Scalar Arguments
        REAL(KIND=prec), INTENT(INOUT) :: &
            epsr           ! Linesearch parameter, used only in null step test to print error
        REAL(KIND=prec), INTENT(IN) :: &
            fo, &          ! Previous value of the objective function.
            tin, &         ! Initial stepsize.
            wk             ! Stopping parameter.
        REAL(KIND=prec), INTENT(OUT) :: &
            f, &           ! Value of the objective function.
            p, &           ! Directional derivative.
            alfn           ! Linearization error.
        INTEGER, INTENT(INOUT) :: &
            nfe, &         ! Number of function evaluations.
            nge, &         ! Number of subgradient evaluations.
            iters, &       ! Null step indicator.
                           !   0  - Null step.
                           !   1  - Serious step.
            inewnma        ! Index for circular arrays.
        INTEGER, INTENT(OUT) :: &
            iterm          ! Cause of termination:
                           !   0  - Everything is ok.
                           !  -3  - Failure in function or subgradient calculations.

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

        ! Local Scalars
        REAL(KIND=prec) :: &
            fu,fu1,fu2, &  ! Auxiliary values of the objective function
            a_orig, &      ! Armijo parameter.
            inv_a_orig, &  ! Inverse of Armijo parameter.
            alpha, &       ! Armijo step size.
            dnorm, &       ! |d|
            epslwk,epsrwk  ! Auxiliary scalars.
        INTEGER :: i,nin   ! Number of interpolations.

        ! Intrinsic Functions
        INTRINSIC ABS,MAX,MAXVAL,DOT_PRODUCT

        ! Initialization

        nin = 0

        epslwk   = epsl*wk
        epsrwk   = epsr*wk
        dnorm    = DOT_PRODUCT(d,d)

        alpha  = two*tin
        a_orig = two
        inv_a_orig = one/a_orig

        ! Updating circular array fold

        IF (iters==1) THEN
            fold(inewnma) = fo
            inewnma = inewnma + 1
            IF (inewnma > mnma) inewnma = 1
        END IF


        ! Function evaluation at a new point

        x = xo + d

        CALL myf(n,x,f,iterm)
        nfe = nfe + 1
        IF (iterm /= 0) THEN
            f=fo
            RETURN
        END IF



        ! Serious step

        iters = 1

        IF (f <= MAXVAL(fold) - epslwk) THEN ! Nonmonotone descent condition

            xu = xo + alpha*d
            CALL myf(n,xu,fu,iterm)
            nfe = nfe + 1
            IF (iterm /= 0) RETURN

            IF(fu > fo - alpha*epslwk .OR. maxnin==0) THEN ! Serious step with t=tin

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

                p = DOT_PRODUCT(g,d)
                alfn = MAX(ABS(fo-f+p*tin),eta*(tin*dnorm)**2)

                RETURN
            END IF

            inct: DO ! Increase t

                nin=nin+1
                f=fu
                x=xu
                alpha = alpha * a_orig
                xu = xo + alpha*d

                CALL myf(n,xu,fu,iterm)
                nfe = nfe + 1
                IF (iterm /= 0) RETURN

                IF (fu > fo - alpha*epslwk .OR. nin >= maxnin) THEN ! Serious step

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

                    p = alpha * inv_a_orig * DOT_PRODUCT(g,d)
                    alfn = MAX(ABS(fo-f+p),eta*(alpha * inv_a_orig*dnorm)**2)

                    RETURN
                END IF
            END DO inct

        ELSE

            alpha = one

            dect: DO ! Decrease t
                nin=nin+1
                IF (nin > maxnin) EXIT dect

                alpha = alpha * inv_a_orig !
                x = xo + alpha * d
                CALL myf(n,x,f,iterm)
                nfe = nfe + 1
                IF (iterm /= 0) RETURN

                IF (f <= MAXVAL(fold) - alpha*epslwk) THEN ! Serious step

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

                    p = alpha * DOT_PRODUCT(g,d)
                    alfn = MAX(ABS(fo-f+p),eta*(alpha*dnorm)**2)

                    RETURN
                END IF
            END DO dect

        END IF


        iters = 0 ! Null step

        ! original d
        !       x = xo + d
        !      CALL myf(n,x,f,iterm)
        !      nfe = nfe + 1


        CALL myg(n,x,g,iterm) ! x = xo + alpha*d
        nge = nge + 1
        IF (iterm /= 0) RETURN

        !       p = DOT_PRODUCT(g,d) ! original d
        p = alpha * DOT_PRODUCT(g,d)
        alfn = MAX(ABS(fo-f+p),eta*(alpha*dnorm)**2)


    END SUBROUTINE nmarmijo




!************************************************************************
!*                                                                      *
!*     * SUBROUTINE dlbfgs *                                            *
!*                                                                      *
!*     Matrix update and computation of the search direction d = -dm*g  *
!*     by the limited memory BFGS update.                               *
!*                                                                      *
!************************************************************************
    
SUBROUTINE dlbfgs(n,mc,mcc,inew,ibfgs,iflag,d,g,gp,s,u,sm,um,rm, &
    umtum,cm,smtgp,umtgp,gamma,tmpn1,iscale)
      
    USE param, ONLY : zero,small,one,half ! half is used at sclpar
    USE subpro, ONLY : &
        xdiffy, & ! Difference of two vectors.
        xsumy, &  ! Sum of two vectors.
        scdiff, & ! Difference of the scaled vector and a vector.
        scsum, &  ! Sum of a vector and the scaled vector.
        vxdiag, & ! Multiplication of a vector and a diagonal matrix.
        symax, &  ! Multiplication of a dense symmetric matrix by a vector.
        cwmaxv, & ! Multiplication of a vector by a dense rectangular matrix.
        rwaxv2, & ! Multiplication of two rowwise stored dense rectangular
                  ! matrices A and B by vectors x and y.
        trlieq, & ! Solving x from linear equation L*x=y or trans(L)*x=y.
        copy2     ! Copying of two vectors.
    IMPLICIT NONE

    ! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
        g, &      ! Current subgradient of the objective function.
        gp        ! Previous subgradient of the objective function.
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
        sm, &     ! Matrix whose columns are stored corrections.
        um, &     ! Matrix whose columns are stored subgradient differences.
        rm, &     ! Upper triangular matrix.
        cm, &     ! Diagonal matrix.
        umtum, &  ! Matrix umtum = trans(um) * um.
        smtgp, &  ! Vector smtgp = trans(sm)*gp.
        umtgp, &  ! vector umtgp = trans(um)*gp.
        s, &      ! Difference of current and previous variables.
        u, &      ! Difference of current and previous subgradients.
        tmpn1     ! Auxiliary array. On input: previous aggregate subgradient.
    REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: &
        d         ! Direction vector.

    ! Scalar Arguments
    REAL(KIND=prec), INTENT(INOUT) :: & 
        gamma     ! Scaling parameter.
    INTEGER, INTENT(IN) :: & 
        n, &      ! Number of variables
        mc, &     ! Declared number of stored corrections.
        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.
                   !   5  - Preliminary scaling with STS/STU.
                   !   6  - No scaling.      
    INTEGER, INTENT(INOUT) :: & 
        mcc, &    ! Current number of stored corrections.
        inew, &   ! Index for circular arrays.
        iflag     ! Index for adaptive version:
                   !   0  - Maximum number of stored corrections
                   !        has not been changed at this iteration.
                   !   1  - Maximum number of stored corrections
                   !        has been changed at this iteration.
    INTEGER, INTENT(OUT) :: & 
        ibfgs     ! Index of the type of BFGS update:
                   !   1  - BFGS update: the corrections are stored.
                   !   2  - BFGS update: the corrections are not stored.
                   !   3  - BFGS update is skipped.
      
    ! Local Arrays
    REAL(KIND=prec), DIMENSION(mcc+1) :: &
        tmpmc1,tmpmc2,tmpmc3,tmpmc4

    ! Local Scalars
    REAL(KIND=prec) :: &
        stu, &    ! stu = trans(s)*u.
        sts       ! sts = trans(s)*s.
    INTEGER :: i,j,k, &
        mcnew, &  ! Current size of vectors.
        iold, &   ! Index of the oldest corrections.
        iflag2, & ! Index for adaptive version.
        ierr      ! Error indicator

    ! Intrinsic Functions
    INTRINSIC SQRT,MIN,MAX


    ierr = 0
    ibfgs = 0
    iflag2 = 0
    stu = DOT_PRODUCT(s,u)
    sts = DOT_PRODUCT(s,s)


    ! Positive definiteness

    IF (stu > zero) THEN
        IF (-DOT_PRODUCT(d,u)-DOT_PRODUCT(tmpn1,s) < -small) THEN
          
     
            ! Update matrices
         
            ibfgs = 1

            ! Initialization of indices.

            CALL indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,ibfgs)
            
     
            ! Update sm and um

            CALL copy2(n,s,sm((inew-1)*n+1:),u,um((inew-1)*n+1:))
           
     
            ! Computation of trans(sm)*g and trans(um)*g

            IF (inew >= mcnew) THEN
                CALL rwaxv2(n,mcnew,sm((inew-mcnew)*n+1:),&
                    um((inew-mcnew)*n+1:),g,g,tmpmc1(iold:),tmpmc2(iold:))
            ELSE
                CALL rwaxv2(n,inew,sm,um,g,g,tmpmc1,tmpmc2)
                CALL rwaxv2(n,mcnew-inew,sm((iold-1)*n+1:),&
                    um((iold-1)*n+1:),g,g,tmpmc1(iold:),tmpmc2(iold:))
            END IF
            
            
            ! Computation of trans(sm)*u and trans(um)*u

            IF (inew >= mcnew) THEN
                DO i=iold,inew-1
                    tmpmc3(i) = tmpmc1(i) - smtgp(i)
                    smtgp(i)  = tmpmc1(i)
                    tmpmc4(i) = tmpmc2(i) - umtgp(i)
                    umtgp(i)  = tmpmc2(i)
                END DO
            ELSE
                DO i=1,inew-1
                    tmpmc3(i) = tmpmc1(i) - smtgp(i)
                    smtgp(i)  = tmpmc1(i)
                    tmpmc4(i) = tmpmc2(i) - umtgp(i)
                    umtgp(i)  = tmpmc2(i)
                END DO
                DO i=iold,mcnew+1
                    tmpmc3(i) = tmpmc1(i) - smtgp(i)
                    smtgp(i)  = tmpmc1(i)
                    tmpmc4(i) = tmpmc2(i) - umtgp(i)
                    umtgp(i)  = tmpmc2(i)
                END DO
            END IF
            tmpmc3(inew) = tmpmc1(inew) - DOT_PRODUCT(s,gp)
            smtgp(inew)  = tmpmc1(inew)
            tmpmc4(inew) = tmpmc2(inew) - DOT_PRODUCT(u,gp)
            umtgp(inew)  = tmpmc2(inew)
            
         
            ! Update rm and umtum

            IF (mcc >= mc .AND. iflag2 /= 1) THEN
                DO i=1,mcnew-1
                    j=(i-1)*i/2+1
                    k=i*(i+1)/2+2
                    CALL copy2(i,rm(k:),rm(j:),umtum(k:),umtum(j:))
                END DO
            END IF
                      
            IF (inew >= mcnew) THEN
                CALL copy2(mcnew,tmpmc3(iold:),rm((mcnew-1)*mcnew/2+1:),&
                    tmpmc4(iold:),umtum((mcnew-1)*mcnew/2+1:))
            ELSE
                CALL copy2(mcnew-inew,tmpmc3(iold:),rm((mcnew-1)*mcnew/2+1:)&
                    ,tmpmc4(iold:),umtum((mcnew-1)*mcnew/2+1:))
                CALL copy2(inew,tmpmc3,rm((mcnew-1)*mcnew/2+mcnew-inew+1:),&
                    tmpmc4,umtum((mcnew-1)*mcnew/2+mcnew-inew+1:))
            END IF
            

            ! Update CM

            cm(inew) = stu
            
            ! Computation of gamma

            gamma = sclpar(mcc,iscale,sts,stu,tmpmc4(inew))
            
            inew = inew + 1
            IF (inew > mc + 1) inew = 1
            IF (iflag == 0 .AND. mcc < mc + 1) mcc = mcc + 1
          
        ELSE

            
            ! BFGS update, corrections are not saved.
     
            ibfgs = 2

            ! Initialization of indices.

            CALL indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,ibfgs)
          
            ! Update sm and um
          
            CALL copy2(n,s,sm((inew-1)*n+1:),u,um((inew-1)*n+1:))
            
            ! Computation of trans(sm)*g and trans(um)*g

            CALL rwaxv2(n,mcnew,sm,um,g,g,tmpmc1,tmpmc2)

            ! Computation of trans(sm)*u and trans(um)*u
          
            IF (iold /= 1) THEN
                DO i=1,inew-1
                    tmpmc3(i) = tmpmc1(i) - smtgp(i)
                    smtgp(i)  = tmpmc1(i)
                    tmpmc4(i) = tmpmc2(i) - umtgp(i)
                    umtgp(i)  = tmpmc2(i)
                END DO
                DO i=iold,mcnew
                    tmpmc3(i) = tmpmc1(i) - smtgp(i)
                    smtgp(i)  = tmpmc1(i)
                    tmpmc4(i) = tmpmc2(i) - umtgp(i)
                    umtgp(i)  = tmpmc2(i)
                END DO
            ELSE
                DO i=1,mcnew-1
                    tmpmc3(i) = tmpmc1(i) - smtgp(i)
                    smtgp(i)  = tmpmc1(i)
                    tmpmc4(i) = tmpmc2(i) - umtgp(i)
                    umtgp(i)  = tmpmc2(i)
                END DO
            END IF
            tmpmc3(inew) = tmpmc1(inew) - DOT_PRODUCT(s,gp)
            smtgp(inew)  = tmpmc1(inew)
            tmpmc4(inew) = tmpmc2(inew) - DOT_PRODUCT(u,gp)
            umtgp(inew)  = tmpmc2(inew)


            ! Update rm and umtum

            IF (iold /= 1) THEN
                CALL copy2(mcnew-inew,tmpmc3(iold:),&
                    rm((mcnew-1)*mcnew/2+1:),tmpmc4(iold:),&
                    umtum((mcnew-1)*mcnew/2+1:))
                CALL copy2(inew,tmpmc3,rm((mcnew-1)*mcnew/2+mcnew-inew+1:),tmpmc4,&
                    umtum((mcnew-1)*mcnew/2+mcnew-inew+1:))
            ELSE
                CALL copy2(mcnew,tmpmc3,rm((mcnew-1)*mcnew/2+1:)&
                    ,tmpmc4,umtum((mcnew-1)*mcnew/2+1:))
            END IF
            

            ! Update cm

            cm(inew) = stu
            

            ! Computation of gamma

            gamma = sclpar(mcc,iscale,sts,stu,tmpmc4(inew))
               
        END IF
    ELSE
         
     
        ! BFGS update is skipped

        ibfgs = 3

        IF (mcc == 0) THEN
            iflag = 0
            d=-g
            RETURN
        END IF
         

        !     Initialization of indices.

        CALL indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,ibfgs)


        !     Computation of gamma

        IF (iscale >= 4) gamma = one
               
         
        !     Computation of trans(sm)*g and trans(um)*g and the two intermediate values

        IF (iold <= 2) THEN
            CALL rwaxv2(n,mcnew,sm((iold-1)*n+1:),um((iold-1)*n+1:),g,g,&
                smtgp(iold:),umtgp(iold:))
        ELSE
            CALL rwaxv2(n,inew-1,sm,um,g,g,smtgp,umtgp)
            CALL rwaxv2(n,mcnew-inew+1,sm((iold-1)*n+1:),&
                um((iold-1)*n+1:),g,g,smtgp(iold:),umtgp(iold:))
        END IF
    END IF


    ! Computation of two intermediate values tmpmc1 and tmpmc2

    IF (iold == 1 .OR. ibfgs == 2) THEN
        CALL trlieq(mcnew,mcnew,iold,rm,tmpmc1,smtgp,1,ierr)
        CALL symax(mcnew,mcnew,iold,umtum,tmpmc1,tmpmc3)
        CALL vxdiag(mcnew,cm,tmpmc1,tmpmc2)
        CALL scsum(mcnew,gamma,tmpmc3,tmpmc2,tmpmc2)
        CALL scsum(mcnew,-gamma,umtgp,tmpmc2,tmpmc3)
        CALL trlieq(mcnew,mcnew,iold,rm,tmpmc2,tmpmc3,0,ierr)

    ELSE IF (iflag == 0) THEN
        CALL trlieq(mcnew,mc+1,iold,rm,tmpmc1,smtgp,1,ierr)
        CALL symax(mcnew,mc+1,iold,umtum,tmpmc1,tmpmc3)
        CALL vxdiag(mc+1,cm,tmpmc1,tmpmc2)
        CALL scsum(mc+1,gamma,tmpmc3,tmpmc2,tmpmc2)
        CALL scsum(mc+1,-gamma,umtgp,tmpmc2,tmpmc3)
        CALL trlieq(mcnew,mc+1,iold,rm,tmpmc2,tmpmc3,0,ierr)

    ELSE
        CALL trlieq(mcnew,mc,iold,rm,tmpmc1,smtgp,1,ierr)
        CALL symax(mcnew,mc,iold,umtum,tmpmc1,tmpmc3)
        CALL vxdiag(mc,cm,tmpmc1,tmpmc2)
        CALL scsum(mc,gamma,tmpmc3,tmpmc2,tmpmc2)
        CALL scsum(mc,-gamma,umtgp,tmpmc2,tmpmc3)
        CALL trlieq(mcnew,mc,iold,rm,tmpmc2,tmpmc3,0,ierr)
    END IF

      
    ! Computation of the search direction d

    IF (iold == 1 .OR. ibfgs == 2) THEN
        CALL cwmaxv(n,mcnew,um,tmpmc1,d)
        CALL xdiffy(n,d,g,d)
        CALL cwmaxv(n,mcnew,sm,tmpmc2,tmpn1)
        CALL scdiff(n,gamma,d,tmpn1,d)
    ELSE 
        CALL cwmaxv(n,inew-1,um,tmpmc1,d)
        CALL cwmaxv(n,mcnew-inew+1,um((iold-1)*n+1:),tmpmc1(iold:),tmpn1)
        CALL xsumy(n,d,tmpn1,d)
        CALL xdiffy(n,d,g,d)
        CALL cwmaxv(n,inew-1,sm,tmpmc2,tmpn1)
        CALL scdiff(n,gamma,d,tmpn1,d)
        CALL cwmaxv(n,mcnew-inew+1,sm((iold-1)*n+1:),tmpmc2(iold:),tmpn1)
        CALL xdiffy(n,d,tmpn1,d)
    END IF

CONTAINS

    FUNCTION sclpar(mcc,iscale,sts,stu,utu) RESULT(spar) ! Calculation of the scaling parameter.
      
        ! USE param, ONLY : small,one,half ! given in host
        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) :: &
            mcc, &     ! Current number of stored corrections.
            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.
                       !   5  - Preliminary 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,4)
                IF (utu < SQRT(small)) THEN
                    spar = one
                    RETURN
                ELSE
                    spar = stu/utu
                END IF
    
            ! Scaling parameter = STS/STU

            CASE(1,3,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
            
        IF (MCC == 0) THEN
            IF (spar < 0.01_prec) spar=0.01_prec
            IF (spar > 100.0_prec) spar=100.0_prec

        ELSE

            SELECT CASE(iscale)

                ! Interval scaling
                CASE(2)
                    IF (spar < 0.6_prec .OR. spar > 6.0_prec) spar = one

                CASE(3)
                    IF (spar < half .OR. spar > 5.0_prec) spar = one
               
                ! Preliminary scaling
                CASE(4,5)
                    spar = one

                ! Scaling at every iteration
                CASE DEFAULT
                CONTINUE
        END SELECT

    END IF

    IF (spar < 1.0E+03_prec*small) spar = 1.0E+03_prec*small
         
END FUNCTION sclpar

END SUBROUTINE dlbfgs

      
!************************************************************************
!*                                                                      *
!*     * SUBROUTINE dlsr1 *                                             *
!*                                                                      *
!*     Matrix update and computation of the search direction d = -dm*ga *
!*     by the limited memory SR1 update.                                *
!*                                                                      *
!************************************************************************

SUBROUTINE dlsr1(n,mc,mcc,inew,isr1,iflag,d,gp,ga,s,u,sm,um,rm,&
    umtum,cm,smtgp,umtgp,gamma,tmpmc1,tmpmc2,tmpn1,nnk,iprint)
      
    USE param, ONLY : zero,small,one
    USE subpro, ONLY : &
        scalex, & ! Scaling a vector.
        xdiffy, & ! Difference of two vectors.
        scdiff, & ! Difference of the scaled vector and a vector.
        xsumy, &  ! Sum of two vectors.
        cwmaxv, & ! Multiplication of a vector by a dense rectangular matrix.
        rwaxv2, & ! Multiplication of two rowwise stored dense rectangular
                  ! matrices A and B by vectors x and y.
        calq, &   ! Solving x from linear equation A*x=y.
        copy, &   ! Copying of a vector.
        copy2     ! Copying of two vectors.
    IMPLICIT NONE

    ! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
        ga, &      ! Current aggregate subgradient of the objective function.
        gp         ! Basic subgradient of the objective function.
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
        sm, &      ! Matrix whose columns are stored corrections.
        um, &      ! Matrix whose columns are stored subgradient differences.
        rm, &      ! Upper triangular matrix.
        cm, &      ! Diagonal matrix.
        umtum, &   ! Matrix umtum = trans(um) * um.
        smtgp, &   ! Vector smtgp = trans(sm)*gp.
        umtgp, &   ! vector umtgp = trans(um)*gp.
        s, &       ! Difference of current and previous variables.
        u, &       ! Difference of current and previous subgradients.
        tmpn1      ! Auxiliary array. On input: previous aggregate subgradient.
    REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: &
        tmpmc1, &  ! Auxiliary array. On output: trans(sm)*ga.
        tmpmc2, &  ! Auxiliary array. On output: trans(um)*ga.
        d          ! Direction vector.

    ! Scalar Arguments
    REAL(KIND=prec), INTENT(INOUT) :: & 
        gamma      ! Scaling parameter.
    INTEGER, INTENT(IN) :: & 
        n, &       ! Number of variables
        mc, &      ! Declared number of stored corrections.
        nnk, &     ! Consecutive null steps counter.
        iprint     ! Printout specification.
    INTEGER, INTENT(INOUT) :: & 
        mcc, &     ! Current number of stored corrections.
        inew, &    ! Index for circular arrays.
        iflag      ! Index for adaptive version:
                   !   0  - Maximum number of stored corrections
                   !        has not been changed at this iteration.
                   !   1  - Maximum number of stored corrections
                   !        has been changed at this iteration.
    INTEGER, INTENT(OUT) :: & 
        isr1       ! Index of the type of L-SR1 update:
                   !   1  - SR1 update: the corrections are stored.
                   !   3  - SR1 update is skipped.
      
    ! Local Arrays
    REAL(KIND=prec), DIMENSION(n) :: tmpn2
    REAL(KIND=prec), DIMENSION((mcc+1)*(mcc+2)/2) :: tmpmat
    REAL(KIND=prec), DIMENSION(mcc+1) :: tmpmc3,tmpmc4,tmpmc5,tmpmc6

    ! Local Scalars
    REAL(KIND=prec) :: &
        stu, &    ! stu = trans(s)*u.
        a, &      ! a = trans(ga) dm_(k-1) ga.
        b         ! b = trans(ga) dm_k ga.
    INTEGER :: i,j,k, &
        mcnew, &  ! Current size of vectors.
        iold, &   ! Index of the oldest corrections.
        iflag2    ! Index for adaptive version.


    iflag2 = 0
    isr1 = 0 
      

    ! Initialization of indices
      
    CALL indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,3)
      
         
    ! Computation of gamma

    gamma = one


    ! Computation of trans(sm)*ga and trans(um)*ga

    IF (iold <= 2) THEN
        CALL rwaxv2(n,mcnew,sm((iold-1)*n+1:),um((iold-1)*n+1:),ga,ga,&
            tmpmc1(iold:),tmpmc2(iold:))
    ELSE
        CALL rwaxv2(n,inew-1,sm,um,ga,ga,tmpmc1,tmpmc2)
        CALL rwaxv2(n,mcnew-inew+1,sm((iold-1)*n+1:),um((iold-1)*n+1:),&
            ga,ga,tmpmc1(iold:),tmpmc2(iold:))
    END IF


    ! Positive definiteness

    IF (-DOT_PRODUCT(d,u) - DOT_PRODUCT(tmpn1,s) >= -small) THEN
      
     
        ! SR1 update is skipped

        isr1 = 3
       
        IF (mcc == 0) THEN
            iflag = 0
            d=-ga
            RETURN
        END IF

    ELSE
      
        stu = DOT_PRODUCT(s,u)
        
        tmpmc1(inew) = DOT_PRODUCT(s,ga)
        tmpmc2(inew) = DOT_PRODUCT(u,ga)
      

        !     Convergence conditions

        IF ((nnk == 1 .OR. mcc < mc) .OR. &
            (iflag == 1 .AND. (inew == 1 .OR. inew == mc))) THEN

            ! SR1 Update

            isr1 = 1
 

            ! Initialization of indices

            CALL indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,1)

            if (iflag2 == 1 .and. iold == 2) then
                tmpmc1(inew) = tmpmc1(1)
                tmpmc2(inew) = tmpmc2(1)
            end if


            ! Update sm and um

            CALL copy2(n,s,sm((inew-1)*n+1:),u,um((inew-1)*n+1:))

         
            ! Update trans(sm)*gp and trans(um)*gp

            smtgp(inew) = DOT_PRODUCT(s,gp)
            umtgp(inew) = DOT_PRODUCT(u,gp)

     
            ! Computation of trans(sm)*u and trans(um)*u

            IF (iold <= 2) THEN
                CALL rwaxv2(n,mcnew-1,sm((iold-1)*n+1:),um((iold-1)*n+1:),&
                    u,u,tmpmc3(iold:),tmpmc4(iold:))
            ELSE
                CALL rwaxv2(n,inew-1,sm,um,u,u,tmpmc3,tmpmc4)
                CALL rwaxv2(n,mcnew-inew,sm((iold-1)*n+1:),um((iold-1)*n+1:),&
                    u,u,tmpmc3(iold:),tmpmc4(iold:))
            END IF

            tmpmc3(inew) = stu
            tmpmc4(inew) = DOT_PRODUCT(u,u)

         
            ! Update rm and umtum

            IF (mcc >= mc .AND. iflag2 /= 1) THEN
                DO i=1,mcnew-1
                    j=(i-1)*i/2+1
                    k=i*(i+1)/2+2
                    CALL copy2(i,rm(k:),rm(j:),umtum(k:),umtum(j:))
                END DO
            END IF


            IF (inew >= mcnew) THEN
                CALL copy2(mcnew,tmpmc3(iold:),rm((mcnew-1)*mcnew/2+1:),&
                    tmpmc4(iold:),umtum((mcnew-1)*mcnew/2+1:))
            ELSE
                CALL copy2(mcnew-inew,tmpmc3(iold:),&
                    rm((mcnew-1)*mcnew/2+1:),tmpmc4(iold:),&
                    umtum((mcnew-1)*mcnew/2+1:))
                CALL COPY2(inew,tmpmc3,rm((mcnew-1)*mcnew/2+mcnew-inew+1:),tmpmc4,&
                    umtum((mcnew-1)*mcnew/2+mcnew-inew+1:))
            END IF

      
            ! Update CM

            cm(inew) = stu

            inew = inew + 1
            IF (inew > mc + 1) inew = 1
            IF (iflag == 0 .AND. mcc < mc + 1) mcc = mcc + 1

        ELSE


            ! Calculation of matrix (umtum-rm-trans(rm)+cm) from previous iteration
         
            DO  i=1,mcnew*(mcnew+1)/2
                tmpmat(i)= gamma * umtum(i) - rm(i)
            END DO

     
            ! Computation of tmpmat*tmpmc4 = gamma*trans(um)*ga-trans(sm)*ga

            IF (iold == 1) THEN
                CALL scdiff(mcnew,gamma,tmpmc2,tmpmc1,tmpmc5)
                CALL calq(mcnew,mcnew,iold,tmpmat,tmpmc4,tmpmc5,0)

            ELSE IF (iflag == 0) THEN
                CALL scdiff(mc+1,gamma,tmpmc2,tmpmc1,tmpmc5)
                CALL calq(mcnew,mc+1,iold,tmpmat,tmpmc4,tmpmc5,0)

            ELSE
                CALL scdiff(mc,gamma,tmpmc2,tmpmc1,tmpmc5)
                CALL calq(mcnew,mc,iold,tmpmat,tmpmc4,tmpmc5,0)
            END IF


            ! Computation of a = -trans(ga)*dm_(k-1)*ga
      
            IF (iold <= 2) THEN
                CALL scalex(mcnew,gamma,tmpmc4(iold:),tmpmc3(iold:))
                CALL cwmaxv(n,mcnew,sm((iold-1)*n+1:),tmpmc4(iold:),tmpn1)
                CALL scdiff(n,-gamma,ga,tmpn1,tmpn2)
                CALL cwmaxv(n,mcnew,um((iold-1)*n+1:),tmpmc3(iold:),tmpn1)
                CALL xsumy(n,tmpn2,tmpn1,tmpn2)
             
            ELSE
                CALL scalex(mcc,gamma,tmpmc4,tmpmc3)
                CALL cwmaxv(n,inew-1,sm,tmpmc4,tmpn1)
                CALL scdiff(n,-gamma,ga,tmpn1,tmpn2)
                CALL cwmaxv(n,mcnew-inew+1,sm((iold-1)*n+1:),tmpmc4(iold:),&
                    tmpn1)
                CALL xdiffy(n,tmpn2,tmpn1,tmpn2)
                CALL cwmaxv(n,inew-1,um,tmpmc3,tmpn1)
                CALL xsumy(n,tmpn2,tmpn1,tmpn2)
                CALL cwmaxv(n,mcnew-inew+1,um((iold-1)*n+1:),tmpmc3(iold:),&
                    tmpn1)
                CALL xsumy(n,tmpn2,tmpn1,tmpn2)
            END IF
          
            a = DOT_PRODUCT(ga,tmpn2)
          
            IF (iflag == 0) THEN
                mcnew = mc
                iold = inew + 2
                IF (iold > mc+1) iold = iold - mc - 1
            ELSE
                mcnew = mc - 1
                iold = inew + 2
                IF (iold > mc) iold = iold - mc
            END IF
      

            ! Calculation of the new canditate for search direction
            ! Updates are not necessarily saved

            ! Update sm and um

            CALL copy2(n,s,sm((inew-1)*n+1:),u,um((inew-1)*n+1:))

     
            ! Computation of trans(sm)*u and trans(um)*u

            IF (iold == 1 .OR. iold == 2) THEN
                CALL rwaxv2(n,mcnew-1,sm((iold-1)*n+1:),um((iold-1)*n+1:),u,u,&
                    tmpmc3(iold:),tmpmc4(iold:))
            ELSE
                CALL rwaxv2(n,inew-1,sm,um,u,u,tmpmc3,tmpmc4)
                CALL rwaxv2(n,mcnew-inew,sm((iold-1)*n+1:),um((iold-1)*n+1:),u,u,&
                    tmpmc3(iold:),tmpmc4(iold:))
            END IF
       
            tmpmc3(inew) = stu
            tmpmc4(inew) = DOT_PRODUCT(u,u)


            ! Calculation of matrix (umtum-rm-trans(rm)+cm) without updating
            ! matrices rm, umtum and cm
      
            DO i=1,mcnew*(mcnew+1)/2
                tmpmat(i)= gamma * umtum(i) - rm(i)
            END DO

            DO i=1,mcnew-1
                j=(i-1)*i/2+1
                k=i*(i+1)/2+2
                CALL copy(i,tmpmat(k:),tmpmat(j:))
            END DO
         
            CALL scdiff(mcnew+1,gamma,tmpmc4,tmpmc3,tmpmc5)
         
            IF (inew >= mcnew) THEN
                CALL copy(mcnew,tmpmc5(iold:),tmpmat((mcnew-1)*mcnew/2+1:))
            ELSE
                CALL copy(mcnew-inew,tmpmc5(iold:),tmpmat((mcnew-1)*mcnew/2+1:))
                CALL copy(inew,tmpmc5,tmpmat((mcnew-1)*mcnew/2+mcnew-inew+1:))
            END IF
      
            IF (iflag == 0) THEN
                CALL scdiff(mc+1,gamma,tmpmc2,tmpmc1,tmpmc5)
                CALL calq(mcnew,mc+1,iold,tmpmat,tmpmc5,tmpmc5,iprint)

            ELSE
                CALL scdiff(mc,gamma,tmpmc2,tmpmc1,tmpmc5)
                CALL calq(mcnew,mc,iold,tmpmat,tmpmc5,tmpmc5,iprint)
            END IF


            ! Calculation of the new canditate for search direction d = -dm_k*ga
            ! and computation of b = -trans(ga)*dm_k*ga
      
            IF (iold <= 2) THEN
                CALL scalex(mcnew,gamma,tmpmc5(iold:),tmpmc6(iold:))
                CALL cwmaxv(n,mcnew,sm((iold-1)*n+1:),tmpmc5(iold:),tmpn1)
                CALL scdiff(n,-gamma,ga,tmpn1,d)
                CALL cwmaxv(n,mcnew,um((iold-1)*n+1:),tmpmc6(iold:),tmpn1)
                CALL xsumy(n,d,tmpn1,d)
         
            ELSE
                CALL scalex(mcnew+1,gamma,tmpmc5,tmpmc6)
                CALL cwmaxv(n,inew,sm,tmpmc5,tmpn1)
                CALL scdiff(n,-gamma,ga,tmpn1,d)
                CALL cwmaxv(n,mcnew-inew,sm((iold-1)*n+1:),tmpmc5(iold:),&
                    tmpn1)
                CALL xdiffy(n,d,tmpn1,d)
                CALL cwmaxv(n,inew,um,tmpmc6,tmpn1)
                CALL xsumy(n,d,tmpn1,d)
                CALL cwmaxv(n,mcnew-inew,um((iold-1)*n+1:),tmpmc6(iold:),&
                    tmpn1)
                CALL xsumy(n,d,tmpn1,d)
            END IF

            b = DOT_PRODUCT(ga,d)


            ! Checking the convergence conditions

            IF (b - a < zero) THEN
                isr1 = 3
                CALL copy(n,tmpn2,d)
            
            ELSE

                isr1 = 1
         
     
                ! Update trans(sm)*gp and trans(um)*gp
     
                smtgp(inew) = DOT_PRODUCT(s,gp)
                umtgp(inew) = DOT_PRODUCT(u,gp)

                     
                ! Update rm and umtum

                DO i=1,mcnew-1
                    j=(i-1)*i/2+1
                    k=i*(i+1)/2+2
                    CALL copy2(i,rm(k:),rm(j:),umtum(k:),umtum(j:))
                END DO

                IF (inew >= mcnew) THEN
                    CALL copy2(mcnew,tmpmc3(iold:),rm((mcnew-1)*mcnew/2+1:),tmpmc4(iold:),&
                        umtum((mcnew-1)*mcnew/2+1:))
                ELSE
                    CALL copy2(mcnew-inew,tmpmc3(iold:),rm((mcnew-1)*mcnew/2+1:),tmpmc4(iold:),&
                        umtum((mcnew-1)*mcnew/2+1:))
                    CALL copy2(inew,tmpmc3,rm((mcnew-1)*mcnew/2+mcnew-inew+1:),tmpmc4,&
                        umtum((mcnew-1)*mcnew/2+mcnew-inew+1:))
                END IF
            

                ! Update cm

                cm(inew) = stu
                     
                inew = inew + 1
                IF (inew > mc + 1) inew = 1
                IF (iflag == 0 .AND. mcc < mc + 1) mcc = mcc + 1
            
            END IF
         
            RETURN
         
        END IF
    END IF
      
    DO i=1,mcnew*(mcnew+1)/2
        tmpmat(i)= gamma * umtum(i) - rm(I)
    END DO
      
     
    ! Computation of tmpmat*tmpmc4 = gamma*trans(um)*ga-trans(sm)*ga

    IF (iold == 1) THEN
        CALL scdiff(mcnew,gamma,tmpmc2,tmpmc1,tmpmc4)
        CALL calq(mcnew,mcnew,iold,tmpmat,tmpmc4,tmpmc4,iprint)
    ELSE IF (iflag == 0) THEN
        CALL scdiff(mc+1,gamma,tmpmc2,tmpmc1,tmpmc4)
        CALL calq(mcnew,mc+1,iold,tmpmat,tmpmc4,tmpmc4,iprint)
    ELSE
        CALL scdiff(mc,gamma,tmpmc2,tmpmc1,tmpmc4)
        CALL calq(mcnew,mc,iold,tmpmat,tmpmc4,tmpmc4,iprint)
    END IF
      

    ! Computation of the search direction d
      
    IF (iold <= 2) THEN
        CALL scalex(mcnew,gamma,tmpmc4(iold:),tmpmc3(iold:))
        CALL cwmaxv(n,mcnew,sm((iold-1)*n+1:),tmpmc4(iold:),tmpn1)
        CALL scdiff(n,-gamma,ga,tmpn1,d)
        CALL cwmaxv(n,mcnew,um((iold-1)*n+1:),tmpmc3(iold:),tmpn1)
        CALL xsumy(n,d,tmpn1,d)
    ELSE
        CALL scalex(mcc,gamma,tmpmc4,tmpmc3)
        CALL cwmaxv(n,inew-1,sm,tmpmc4,tmpn1)
        CALL scdiff(n,-gamma,ga,tmpn1,d)
        CALL cwmaxv(n,mcnew-inew+1,sm((iold-1)*n+1:),tmpmc4(iold:),&
            tmpn1)
        CALL xdiffy(n,d,tmpn1,d)
        CALL cwmaxv(n,inew-1,um,tmpmc3,tmpn1)
        CALL xsumy(n,d,tmpn1,d)
        CALL cwmaxv(n,mcnew-inew+1,um((iold-1)*n+1:),tmpmc3(iold:),&
            tmpn1)
        CALL xsumy(n,d,tmpn1,d)
    END IF
      
END SUBROUTINE dlsr1

!************************************************************************
!*                                                                      *
!*     * SUBROUTINE indic1 *                                            *
!*                                                                      *
!*     Initialization of indices.                                       *
!*                                                                      *
!************************************************************************

SUBROUTINE indic1(mc,mcc,mcnew,inew,iold,iflag,iflag2,itype)

    IMPLICIT NONE

    ! Scalar Arguments

    INTEGER, INTENT(IN) :: & 
        mc, &     ! Declared number of stored corrections.
        mcc, &    ! Current number of stored corrections.
        itype     ! Type of Initialization:
                   !   1  - corrections are stored,
                   !   2  - corrections are not stored,
                   !   3  - update is skipped.
    INTEGER, INTENT(INOUT) :: & 
        inew, &   ! Index for circular arrays.
        iflag, &  ! Index for adaptive version:
                  !   0  - Maximum number of stored corrections
                  !        has not been changed at this iteration.
                  !   1  - Maximum number of stored corrections
                  !        has been changed at this iteration.
        iflag2    ! Index for adaptive version.
                   !   0  - iflag has not been changed.
                   !   1  - iflag has been changed.
    INTEGER, INTENT(OUT) :: & 
        mcnew, &  ! Current size of vectors.
        iold      ! Index of the oldest corrections.
      
    IF (itype == 1) THEN
        IF (mcc < mc) THEN
            mcnew = mcc + 1
            iold = 1
            iflag = 0
        ELSE
            IF (iflag == 0) THEN
                mcnew = mc
                iold = inew + 2
                IF (iold > mc+1) iold = iold - mc - 1
            ELSE
                IF (inew == 1) THEN
                    inew = mc + 1
                    mcnew = mc
                    iold = 2
                    iflag = 0
                    iflag2 = 1
                ELSE IF (inew == mc) THEN
                    mcnew = mc
                    iold = 1
                    iflag = 0
                    iflag2 = 1
                ELSE
                    mcnew = mc - 1
                    iold = inew + 2
                    IF (iold > mc) iold = iold - mc
                END IF
            END IF
        END IF
      
    ELSE IF (itype == 2) THEN
        IF (mcc < mc) THEN
            mcnew = mcc + 1
            iold = 1
            iflag = 0
        ELSE
            IF (iflag == 0) THEN
                mcnew = mc + 1
                iold = inew + 1
                IF (iold > mc + 1) iold = 1
            ELSE
                mcnew = mc
                iold = inew + 1
                IF (iold > mc) iold = 1
            END IF
        END IF

    ELSE 
        IF (mcc < mc) THEN
            mcnew = mcc
            iold = 1
            iflag = 0
        ELSE
            IF (iflag == 0) THEN
                mcnew = mc
                iold = inew + 1
                IF (iold > mc + 1) iold = 1
            ELSE
                mcnew = mc - 1
                iold = inew + 1
                IF (iold > mc) iold = 1
            END IF
        END IF
    END IF
END SUBROUTINE indic1
  
!************************************************************************
!*
!*     * SUBROUTINE agbfgs *
!*
!*     Computation of aggregate values by the limited memory BFGS update.
!*
!************************************************************************

SUBROUTINE agbfgs(n,mc,mcc,inew,ibfgs,iflag,g,gp,ga,u,d,sm,um, &
    rm,cm,umtum,alfn,alfv,gamma,ic,rho)

    USE param, ONLY : zero,half,one
    USE subpro, ONLY : &
        symax, &  ! Multiplication of a dense symmetric matrix by a vector.
        rwaxv2, & ! Multiplication of two rowwise stored dense rectangular
                  ! matrices A and B by vectors x and y.
        trlieq, & ! Solving x from linear equation L*x=y or trans(L)*x=y.
        vdot      ! Dot product of vectors

    IMPLICIT NONE

    ! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
        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.
        sm, &     ! Matrix whose columns are stored corrections.
        um, &     ! Matrix whose columns are stored subgradient differences.
        rm, &     ! Upper triangular matrix.
        umtum, &  ! Matrix umtum = trans(um) * um.
        cm        ! Diagonal matrix.
    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.
        rho       ! Correction parameter.
    INTEGER, INTENT(IN) :: & 
        n, &      ! Number of variables
        mc, &     ! Declared number of stored corrections.
        mcc, &    ! Current number of stored corrections.
        inew, &   ! Index for circular arrays.
        ibfgs, &  ! Index of the type of BFGS update.
        ic        ! Correction indicator.
    INTEGER, INTENT(INOUT) :: & 
        iflag     ! Index for adaptive version:
                   !   0  - Maximum number of stored corrections
                   !        has not been changed at this iteration.
                   !   1  - Maximum number of stored corrections
                   !        has been changed at this iteration.

    ! Local arrays
    REAL(KIND=prec), DIMENSION(mcc+1) :: tmpmc1, tmpmc2

    ! 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, &
        mcnew, &  ! Current size of vectors.
        iold, &   ! Index of the oldest corrections.
        ierr      ! Error indicador.

    ! Intrinsic Functions
    INTRINSIC MAX,MIN,SIGN

    ierr = 0

    IF (mcc < mc) THEN
        IF (ibfgs == 2) THEN
            mcnew = mcc + 1
        ELSE
            mcnew = mcc
        END IF
        iold = 1

    ELSE
        IF (iflag == 0) THEN
            IF (ibfgs == 2) THEN
                mcnew = mc + 1
            ELSE
                mcnew = mc
            END IF
            iold = inew + 1
            IF (iold > mc+1) iold = 1

        ELSE
            IF (ibfgs == 2) THEN
                mcnew = mc
            ELSE
                mcnew = mc - 1
            END IF
            iold = inew + 1
            IF (iold > mc) iold = 1
        END IF
    END IF
      
      
    ! Computation of trans(d) * u - alfn

    p = DOT_PRODUCT(d,u) - alfn
    q = DOT_PRODUCT(u,u)

    IF (ic == 1) THEN
        w = rho * q
    ELSE
        w = zero
    END IF
         
     
    ! Computation of the product trans(u)*dm*u

    IF (mcc > 0 .OR. ibfgs == 2) THEN

        IF (iold == 1 .OR. ibfgs == 2) THEN
            CALL rwaxv2(n,mcnew,sm,um,u,u,tmpmc1,tmpmc2)
            CALL trlieq(mcnew,mcnew,iold,rm,tmpmc1,tmpmc1,1,ierr)

            q = q - 2.0_prec*vdot(mcnew,tmpmc2,tmpmc1)
            q = gamma*q
            
            DO i=1,mcnew
                tmpmc2(i) = cm(i)*tmpmc1(i)
            END DO

            q = q + vdot(mcnew,tmpmc1,tmpmc2)

            CALL symax(mcnew,mcnew,iold,umtum,tmpmc1,tmpmc2)

            q = q + gamma*vdot(mcnew,tmpmc1,tmpmc2)

        ELSE
            CALL rwaxv2(n,inew-1,sm,um,u,u,tmpmc1,tmpmc2)
            CALL rwaxv2(n,mcc-inew,sm((iold-1)*n+1:),um((iold-1)*n+1:), &
                u,u,tmpmc1(iold:),tmpmc2(iold:))
            CALL trlieq(mcnew,mcc,iold,rm,tmpmc1,tmpmc1,1,ierr)

            q = q - 2.0_prec*(vdot(mcc-inew,tmpmc2(iold:),tmpmc1(iold:)) + &
                vdot(inew-1,tmpmc2,tmpmc1))
            q = gamma*q

            DO i=1,mcc
                tmpmc2(i) = cm(i)*tmpmc1(i)
            END DO

            q = q + vdot(mcc-inew,tmpmc1(iold:),tmpmc2(iold:)) + &
                vdot(inew-1,tmpmc1,tmpmc2)

            CALL symax(mcnew,mcc,iold,umtum,tmpmc1,tmpmc2)

            q = q + gamma*(vdot(mcc-inew,tmpmc1(iold:),tmpmc2(iold:)) + &
                vdot(inew-1,tmpmc1,tmpmc2))
        END IF

    END IF
    
    q = q + w
    
    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 agbfgs

!************************************************************************
!*
!*     * SUBROUTINE aggsr1 *
!*
!*     Computation of aggregate values by the limited memory SR1 update.
!*
!************************************************************************
      
SUBROUTINE aggsr1(n,mc,mcc,inew,iflag,g,gp,ga,d,alfn,alfv, &
    umtum,rm,gamma,smtgp,umtgp,smtga,umtga,sm,um,icn,rho)
      
    USE param, ONLY : zero,one,small
    USE subpro, ONLY : &
        vdot, &   ! Dot product.
        scalex, & ! Scaling a vector.
        xsumy, &  ! Sum of two vectors.
        xdiffy, & ! Difference of two vectors.
        scsum, &  ! Sum of a vector and the scaled vector.
        scdiff, & ! Difference of the scaled vector and a vector.
        rwaxv2, & ! Multiplication of two rowwise stored dense rectangular
                  ! matrices A and B by vectors x and y.
        cwmaxv, & ! Multiplication of a vector by a dense rectangular matrix.
        lineq, &  ! Solver for linear equation.
        calq      ! Solving x from linear equation A*x=y.
    IMPLICIT NONE

    ! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
        d, &      ! Direction vector.
        g, &      ! Current (auxiliary) subgradient of the objective function.
        gp, &     ! Previous subgradient of the objective function.
        sm, &     ! Matrix whose columns are stored corrections.
        um, &     ! Matrix whose columns are stored subgradient differences.
        rm, &     ! Upper triangular matrix.
        umtum, &  ! Matrix umtum = trans(um) * um.
        smtgp, &  ! Vector smtgp = trans(sm)*gp.
        umtgp, &  ! vector umtgp = trans(um)*gp.
        smtga, &  ! vector smtga = trans(sm)*ga.
        umtga     ! vector umtga = trans(um)*ga.
    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.
        rho       ! Correction parameter.
    INTEGER, INTENT(IN) :: & 
        n, &      ! Number of variables
        mc, &     ! Declared number of stored corrections.
        mcc, &    ! Current number of stored corrections.
        inew, &   ! Index for circular arrays.
        icn       ! Correction indicator.
    INTEGER, INTENT(INOUT) :: & 
        iflag     ! Index for adaptive version:
                   !   0  - Maximum number of stored corrections
                   !        has not been changed at this iteration.
                   !   1  - Maximum number of stored corrections
                   !        has been changed at this iteration.
      
    ! Local arrays
    REAL(KIND=prec), DIMENSION(n) :: tmpn2,tmpn3,tmpn4
    REAL(KIND=prec), DIMENSION((mcc+1)*(mcc)/2) :: tmpmat
    REAL(KIND=prec), DIMENSION(mcc+1) :: tmpmc3, tmpmc4

    ! 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.
        w, &      ! Correction.
        tmp1, &   ! Auxiliary scalar.
        tmp2      ! Auxiliary scalar.
    INTEGER :: i, &
        mcnew, &  ! Current size of vectors.
        iold, &   ! Index of the oldest corrections.
        ierr      ! Error indicador.

      
    ! Intrinsic Functions
    INTRINSIC MIN,MAX,DOT_PRODUCT


    ierr = 0
      
    IF (mcc < mc) THEN
        iold = 1
        mcnew = mcc
    ELSE IF (iflag == 0) THEN
        mcnew = mc
        iold = inew + 1
        IF (iold > mc+1) iold = 1
    ELSE
        mcnew = mc - 1
        iold = inew + 1
        IF (iold > mc) iold = 1
    END IF
      
    CALL xdiffy(n,gp,ga,tmpn2)
     
      
    ! Calculation of tmpn3 = trans(gp - ga)dm

    IF (mcc > 0) THEN

        DO i=1,mcnew*(mcnew+1)/2
            tmpmat(i)= gamma * umtum(i) - rm(i)
        END DO

        IF (iold == 1) THEN
            CALL xdiffy(mcnew,umtgp,umtga,tmpmc4)
            CALL scdiff(mcnew,gamma,tmpmc4,smtgp,tmpmc4)
            CALL xsumy(mcnew,tmpmc4,smtga,tmpmc4)
          
            CALL calq(mcnew,mcnew,iold,tmpmat,tmpmc3,tmpmc4,0)
            CALL scalex(mcnew,gamma,tmpmc3,tmpmc4)
          
            CALL cwmaxv(n,mcnew,sm,tmpmc3,tmpn4)
            CALL scsum(n,gamma,tmpn2,tmpn4,tmpn3)
            CALL cwmaxv(n,mcnew,um,tmpmc4,tmpn4)
            CALL xdiffy(n,tmpn3,tmpn4,tmpn3)

        ELSE
            CALL xdiffy(mcc,umtgp,umtga,tmpmc4)
            CALL scdiff(mcc,gamma,tmpmc4,smtgp,tmpmc4)
            CALL xsumy(mcc,tmpmc4,smtga,tmpmc4)

            CALL calq(mcnew,mcc,iold,tmpmat,tmpmc3,tmpmc4,0)
            CALL scalex(mcc,gamma,tmpmc3,tmpmc4)
          
            CALL cwmaxv(n,inew-1,sm,tmpmc3,tmpn4)
            CALL scsum(n,gamma,tmpn2,tmpn4,tmpn3)
            CALL cwmaxv(n,mcnew-inew+1,sm((iold-1)*n+1:),tmpmc3(iold:)&
                ,tmpn4)
            CALL xsumy(n,tmpn3,tmpn4,tmpn3)
            CALL cwmaxv(n,inew-1,um,tmpmc4,tmpn4)
            CALL xdiffy(n,tmpn3,tmpn4,tmpn3)
            CALL cwmaxv(n,mcnew-inew+1,um((iold-1)*n+1:),tmpmc4(iold:)&
                ,tmpn4)
            CALL xdiffy(n,tmpn3,tmpn4,tmpn3)
        END IF
       
        IF (icn == 1) THEN
            CALL scsum(n,rho,tmpn2,tmpn3,tmpn3)
        END IF
         
        pr = DOT_PRODUCT(tmpn3,tmpn2)
        rrp = DOT_PRODUCT(tmpn3,ga)
        CALL xdiffy(n,g,ga,tmpn4)
        prqr = DOT_PRODUCT(tmpn3,tmpn4)
        rrq = -DOT_PRODUCT(tmpn4,d)

    ELSE

        pr = DOT_PRODUCT(tmpn2,tmpn2)
        rrp = DOT_PRODUCT(tmpn2,ga)
        CALL xdiffy(n,g,ga,tmpn4)
        prqr = DOT_PRODUCT(tmpn2,tmpn4)
        rrq = -DOT_PRODUCT(tmpn4,d)
    END IF

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

    qr = DOT_PRODUCT(tmpn4,tmpn4)
    IF (icn == 1) THEN
        w = rho*qr
    ELSE
        w = zero
    END IF
      
    IF (mcc > 0) THEN
        qr = gamma*qr

        IF (iold == 1) THEN
            CALL rwaxv2(n,mcnew,sm,um,tmpn4,tmpn4,tmpmc4,tmpmc3)
            CALL scsum(mcnew,-gamma,tmpmc3,tmpmc4,tmpmc4)
            CALL lineq(mcnew,mcnew,iold,tmpmat,tmpmc3,tmpmc4,ierr)
            
            qr = qr - vdot(mcnew,tmpmc4,tmpmc3) + w

        ELSE
            CALL rwaxv2(n,inew-1,sm,um,tmpn4,tmpn4,tmpmc4,tmpmc3)
            CALL rwaxv2(n,mcnew-inew+1,sm((iold-1)*n+1:),&
                um((iold-1)*n+1:),tmpn4,tmpn4,tmpmc4(iold:),tmpmc3(iold:))
            CALL scsum(mcc,-gamma,tmpmc3,tmpmc4,tmpmc4)
            CALL lineq(mcnew,mcc,iold,tmpmat,tmpmc3,tmpmc4,ierr)
          
            qr = qr - vdot(mcc-inew,tmpmc4(iold:),tmpmc3(iold:)) -&
                vdot(inew-1,tmpmc4,tmpmc3) + w
        END IF

    END IF
      
    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
        w = pr - prqr*tmp2
        IF (w /= zero) THEN
            lam1 = (tmp1*prqr - rrp)/w
            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))
    w = (lam2*qr + rrq+rrq)*lam2
    !    w = (lam2*qr + 2.0_prec*rrq)*lam2
    tmp1 = zero
    IF (alfv >= zero) tmp1 = one
    IF (pr > zero) tmp1 = MIN(one,MAX(zero,-rrp/pr))
    !    tmp2 = (tmp1*pr + 2.0_prec*rrp)*tmp1
    tmp2 = (tmp1*pr + rrp+rrp)*tmp1
    IF (tmp2 < w) THEN
        w = tmp2
        lam1 = tmp1
        lam2 = zero
    END IF
    
    IF (qqp*(qqp - pq) < zero) THEN
        IF (qr + rrq + rrq - qqp*qqp/pq < W) 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 aggsr1
      
!************************************************************************
!*
!*     * SUBROUTINE agskip *
!*
!*     Computation of aggregate values after consecutive null steps
!*     by the limited memory BFGS update.
!*
!************************************************************************
      
SUBROUTINE agskip(n,mc,mcc,inew,iflag,g,gp,ga,d,u,alfn,alfv, &
    umtum,rm,cm,gamma,smtgp,umtgp,smtga,umtga,sm,um,icn,rho)
      
    USE param, ONLY : zero,half,one,small
    USE subpro, ONLY : &
        xdiffy, & ! Difference of two vectors.
        symax, &  ! Multiplication of a dense symmetric matrix by a vector.
        rwaxv2, & ! Multiplication of two rowwise stored dense rectangular
                  ! matrices A and B by vectors x and y.
        trlieq, & ! Solving x from linear equation l*x=y or trans(l)*x=y.
        vdot      ! Dot product
    IMPLICIT NONE

    ! Array Arguments
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
        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.
        sm, &     ! Matrix whose columns are stored corrections.
        um, &     ! Matrix whose columns are stored subgradient differences.
        rm, &     ! Upper triangular matrix.
        cm, &     ! Diagonal matrix.
        umtum, &  ! Matrix umtum = trans(um) * um.
        smtgp, &  ! Vector smtgp = trans(sm)*gp.
        umtgp, &  ! vector umtgp = trans(um)*gp.
        smtga, &  ! vector smtga = trans(sm)*ga.
        umtga     ! vector umtga = trans(um)*ga.
    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.
        rho       ! Correction parameter.
    INTEGER, INTENT(IN) :: & 
        n, &      ! Number of variables
        mc, &     ! Declared number of stored corrections.
        mcc, &    ! Current number of stored corrections.
        inew, &   ! Index for circular arrays.
        icn       ! Correction indicator.
    INTEGER, INTENT(INOUT) :: & 
        iflag     ! Index for adaptive version:
                   !   0  - Maximum number of stored corrections
                   !        has not been changed at this iteration.
                   !   1  - Maximum number of stored corrections
                   !        has been changed at this iteration.
      
    ! Local arrays
    REAL(KIND=prec), DIMENSION(n) :: tmpn2
    REAL(KIND=prec), DIMENSION(mcc+1) :: tmpmc3, tmpmc4

    ! 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.
        w, &      ! Correction.
        tmp1, &   ! Auxiliary scalar.
        tmp2      ! Auxiliary scalar.
    INTEGER :: i, &
        mcnew, &  ! Current size of vectors.
        iold, &   ! Index of the oldest corrections.
        ierr      ! Error indicador.
    
    ! Intrinsic Functions
    INTRINSIC MIN,MAX,DOT_PRODUCT

    ierr = 0
           
    IF (mcc < mc) THEN
        iold = 1
        mcnew = mcc
    ELSE
        IF (iflag == 0) THEN
            mcnew = mc
            iold = inew + 1
            IF (iold > mc+1) iold = 1
        ELSE
            mcnew = mc - 1
            iold = inew + 1
            IF (iold > mc) iold = 1
        END IF
    END IF

      
    ! Calculation of pq = trans(g-gp) dm (g-gp) = trans(u) dm u.

    pq = DOT_PRODUCT(u,u)

    IF (icn == 1) THEN
        w = rho * pq
    ELSE
        w = zero
    END IF
    
    IF (mcc > 0) THEN

        IF (iold == 1) THEN
            CALL rwaxv2(n,mcnew,sm,um,u,u,tmpmc3,tmpmc4)
            CALL trlieq(mcnew,mcnew,iold,rm,tmpmc3,tmpmc3,1,ierr)
          
            pq = pq - 2.0_prec*vdot(mcnew,tmpmc4,tmpmc3)
            pq = gamma*pq
          
            DO i=1,mcnew
                tmpmc4(i) = cm(i)*tmpmc3(i)
            END DO

            pq = pq + vdot(mcnew,tmpmc3,tmpmc4)
          
            CALL symax(mcnew,mcnew,iold,umtum,tmpmc3,tmpmc4)
          
            pq = pq + gamma*vdot(mcnew,tmpmc3,tmpmc4)

        ELSE
            CALL rwaxv2(n,inew-1,sm,um,u,u,tmpmc3,tmpmc4)
            CALL rwaxv2(n,mcc-inew,sm((iold-1)*n+1:),um((iold-1)*n+1:),&
                u,u,tmpmc3(iold:),tmpmc4(iold:))
            CALL trlieq(mcnew,mcc,iold,rm,tmpmc3,tmpmc3,1,ierr)

            pq = pq - 2.0_prec*(vdot(mcc-inew,tmpmc4(iold:),tmpmc3(iold:)) + &
                vdot(inew-1,tmpmc4,tmpmc3))
            pq = gamma*pq
          
            DO i=1,mcc
                tmpmc4(i) = cm(i)*tmpmc3(i)
            END DO

            pq = pq + vdot(mcc-inew,tmpmc3(iold:),tmpmc4(iold:)) + &
                vdot(inew-1,tmpmc3,tmpmc4)

            CALL symax(mcnew,mcc,iold,umtum,tmpmc3,tmpmc4)
          
            pq = pq + gamma*(vdot(mcc-inew,tmpmc3(iold:),tmpmc4(iold:)) &
                + vdot(inew-1,tmpmc3,tmpmc4))
        END IF

    END IF

    pq = pq + w
      

    ! Calculation of pr = trans(gp-ga) dm (gp-ga).
      
    CALL xdiffy(n,gp,ga,tmpn2)
    pr = DOT_PRODUCT(tmpn2,tmpn2)

    IF (icn == 1) THEN
        w = rho * pr
    ELSE
        w = zero
    END IF

    IF (mcc > 0) THEN
       
        IF (iold == 1) THEN
            DO i=1, mcnew
                tmpmc3(i)=smtgp(i)-smtga(i)
                tmpmc4(i)=umtgp(i)-umtga(i)
            END DO
            CALL trlieq(mcnew,mcnew,iold,rm,tmpmc3,tmpmc3,1,ierr)
             
            pr = pr - 2.0_prec*vdot(mcnew,tmpmc4,tmpmc3)
            pr = gamma*pr
            
            DO i=1,mcnew
                tmpmc4(i) = cm(i)*tmpmc3(i)
            END DO

            pr = pr + vdot(mcnew,tmpmc3,tmpmc4)

            CALL symax(mcnew,mcnew,iold,umtum,tmpmc3,tmpmc4)

            pr = pr + gamma*vdot(mcnew,tmpmc3,tmpmc4)

        ELSE
            DO i=1, mcc
                tmpmc3(i)=smtgp(i)-smtga(i)
                tmpmc4(i)=umtgp(i)-umtga(i)
            END DO
            CALL trlieq(mcnew,mcc,iold,rm,tmpmc3,tmpmc3,1,ierr)

            pr = pr - 2.0_prec*(vdot(mcc-inew,tmpmc4(iold:),tmpmc3(iold:)) + &
                vdot(inew-1,tmpmc4,tmpmc3))
            pr = gamma*pr

            DO  i=1,mcc
                tmpmc4(i) = cm(i)*tmpmc3(i)
            END DO

            pr = pr + vdot(mcc-inew,tmpmc3(iold:),tmpmc4(iold:)) + &
                vdot(inew-1,tmpmc3,tmpmc4)

            CALL symax(mcnew,mcc,iold,umtum,tmpmc3,tmpmc4)

            pr = pr + gamma*(vdot(mcc-inew,tmpmc3(iold:),tmpmc4(iold:)) &
                + vdot(inew-1,tmpmc3,tmpmc4))
        END IF

    END IF

    pr = pr + w

      
    ! Calculation of rrp = trans(gp-ga) dm ga - alfv.
      
    rrp = - DOT_PRODUCT(tmpn2,d) - alfv
      

    ! Calculation of qr = trans(g-ga) dm (g-ga).

    CALL xdiffy(n,g,ga,tmpn2)
    qr = DOT_PRODUCT(tmpn2,tmpn2)

    IF (icn == 1) THEN
        w = rho * qr
    ELSE
        w = zero
    END IF

    IF (mcc > 0) THEN

        IF (iold == 1) THEN
            CALL rwaxv2(n,mcnew,sm,um,tmpn2,tmpn2,tmpmc3,tmpmc4)
            CALL trlieq(mcnew,mcnew,iold,rm,tmpmc3,tmpmc3,1,ierr)

            qr = qr - 2.0_prec*vdot(mcnew,tmpmc4,tmpmc3)
            qr = gamma*qr
            
            DO i=1,mcnew
                tmpmc4(i) = cm(i)*tmpmc3(i)
            END DO

            qr = qr + vdot(mcnew,tmpmc3,tmpmc4)
          
            CALL symax(mcnew,mcnew,iold,umtum,tmpmc3,tmpmc4)
          
            qr = qr + gamma*vdot(mcnew,tmpmc3,tmpmc4)

        ELSE
            CALL rwaxv2(n,inew-1,sm,um,tmpn2,tmpn2,tmpmc3,tmpmc4)
            CALL rwaxv2(n,mcc-inew,sm((iold-1)*n+1:),um((iold-1)*n+1:),&
                tmpn2,tmpn2,tmpmc3(iold:),tmpmc4(iold:))
            CALL trlieq(mcnew,mcc,iold,rm,tmpmc3,tmpmc3,1,ierr)

            qr = qr - 2.0_prec*(vdot(mcc-inew,tmpmc4(iold:),tmpmc3(iold:)) + &
                vdot(inew-1,tmpmc4,tmpmc3))
            qr = gamma*qr
          
            DO i=1,mcc
                tmpmc4(i) = cm(i)*tmpmc3(i)
            END DO
          
            qr = qr + vdot(mcc-inew,tmpmc3(iold:),tmpmc4(iold:)) + &
                vdot(inew-1,tmpmc3,tmpmc4)
          
            CALL symax(mcnew,mcc,iold,umtum,tmpmc3,tmpmc4)
          
            qr = qr + gamma*(vdot(mcc-inew,tmpmc3(iold:),tmpmc4(iold:)) &
                +vdot(inew-1,tmpmc3,tmpmc4))
        END IF
       
    END IF
    
    qr = qr + w
      

    ! Calculation of rrq = trans(g-ga) dm ga - alfv + alfn.

    rrq = - DOT_PRODUCT(tmpn2,d) - alfv + alfn

     
    ! Calculation of prqr = trans(gp-ga) dm (g-ga).
      
    prqr = half*(qr - pq + pr)

     
    ! Calculation of qqp = trans(g-gp) dm g + alfn.

    qqp = pq + prqr + rrq - pr - rrp

     
    ! Computation of multipliers lam1 and lam2

    IF (pr > zero .AND. qr > zero) THEN
        tmp1 = rrq/qr
        tmp2 = prqr/qr
        w = pr - prqr*tmp2
        IF (w /= zero) THEN

            lam1 = (tmp1*prqr - rrp)/w
            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

    lam1 = zero
    lam2 = zero
    IF (alfn <= alfv) lam2 = one
    IF (qr > zero) lam2 = MIN(one,MAX(zero,-rrq/qr))
    w = (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 < w) THEN
        w = tmp2
        lam1 = tmp1
        lam2 = zero
    END IF
      
    IF (qqp*(qqp - pq) < zero) THEN
        IF (qr + rrq + rrq - qqp*qqp/pq < w) 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 agskip

!************************************************************************
!*
!*     * 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.
                      !  -1  - Two consecutive restarts.
                      !  -2  - Number of restarts > maximum number of restarts.
                      !  -3  - Failure in function or subgradient calculations 
                      !        (assigned by the user).
                      !  -4  - Failure in attaining the demanded accuracy.
                      !  -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 == -4) 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.
                      !  -1  - Two consecutive restarts.
                      !  -2  - Number of restarts > maximum number of restarts.
                      !  -3  - Failure in function or subgradient calculations 
                      !        (assigned by the user).
                      !  -4  - Failure in attaining the demanded accuracy.
                      !  -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 lmbm_mod
