MODULE smdb_mod  ! Splitting Metrics D-Bundle Method

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

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

  PUBLIC :: &
       init_smdb, &  ! S Initialization for SMDB subroutine.
       smdb          ! S SMDB subroutine for nonsmooth large-scale
                     !   optimization. Contains:
                     !     S restar     Initialization and reinitialization.
                     !     S armijo     Armijo line search.
                     !     S nmarmijo   Nonmonotone Armijo line search.
                     !     S aggre1     Simplified subgradient aggregation.
                     !     S aggre2     Subgradient aggregation.
  PRIVATE :: &
       wprint, &     ! S Printout the error and warning messages.
       rprint        ! S Printout the results.

CONTAINS

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

    IMPLICIT NONE

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

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

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

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

! Default values
    IF (mit   <= 0) mit      = 10000           ! Maximum number of iterations.
    IF (mfe   <= 0) mfe      = n*mit           ! Maximum number of function evaluations.
    IF (tolf  <= zero) tolf  = 1.0E-08_prec    ! Tolerance for change of function values.
    IF (tolf2 == zero) tolf2 = 1.0E+04_prec    ! Second tolerance for change of function values.
    IF (tolb  == zero) tolb  = -large + small  ! Tolerance for the function value.
    IF (tolg  <= zero) tolg  = 1.0E-06_prec    ! Tolerance for the termination criterion.
    IF (told  <= zero) told  = one             ! Tolerance for the diagonal elements.
    IF (eta   <  zero) eta   = 1.0E-03_prec    ! Distance measure parameter
    IF (epsl  <= zero) epsl  = 1.0E-04_prec    ! Line search parameter,
    IF (mtesf <= 0) mtesf    = 10              ! Maximum number of iterations with changes 
                                               ! of function values smaller than tolf.

     
  END SUBROUTINE init_smdb


!***********************************************************************
!*                                                                     *
!*     * SUBROUTINE smdb *                                             *
!*                                                                     *
!*     Splitting metrics diagonal bundle subroutine                    *
!*     for nonconvex nonsmooth optimization.                           *
!*                                                                     *
!***********************************************************************
      
  SUBROUTINE smdb(f,nit,nfe,nge,iterm,strtim)
    USE param, ONLY : small,large,zero,half,one  ! Parameters.
    USE exe_time, ONLY : getime  ! Execution time.
    USE initializat, ONLY : &
         n, &         ! Number of variables.
         x, &         ! Vector of variables
         mcu, &       ! Maximum number of stored corrections, mcu >= 1.
         mnma, &      ! Maximum number of function values used in nonmonotone line search.
         inma, &      ! Selection of line search method:
                      !   nma = 0, Armijo line search,
                      !   nma = 1, nonmonotone Armijo line search.
         iprint, &    ! Printout specification (see initializat for more details).
         tolf, &      ! Tolerance for change of function values.
         tolf2, &     ! Second tolerance for change of function values.
         tolb, &      ! Tolerance for the function value.
         tolg, &      ! Tolerance for the first termination criterion. 
         told, &      ! Upper bound for the diagonal elements of diag.
         mintold, &   ! Lower bound for the diagonal elements of diag.
         eta, &       ! Distance measure parameter, eta >= 0. 
         epsl, &      ! Line search parameter, 0 < epsl < 0.25.
         mtesf, &     ! Maximum number of iterations with changes of
                      ! function values smaller than tolf.
         mit, &       ! Maximun number of iterations.
         mfe, &       ! Maximun number of function evaluations.
         time         ! Maximum time
    USE obj_fun, ONLY : &
         myf, &       ! Computation of the value of the objective function. 
         myg          ! Computation of the subgradient of the objective function. 
    USE subpro, ONLY : &
         vxdiag, &    ! Vector is multiplied by a diagonal matrix y:=diag*x.
         xdiffy, &    ! Difference of two vectors.
         copy, &      ! Copying of a vector.
         copy2        ! Copying of two vectors.


    IMPLICIT NONE

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


! Local Arrays
    REAL(KIND=prec), DIMENSION(n) :: &
         xo, &        ! Previous vector of variables.
         g, &         ! Subgradient of the objective function.
         gp, &        ! Previous subgradient of the objective function.
         ga, &        ! Aggregate subgradient.
         s, &         ! Difference of current and previous variables.
         u, &         ! Difference of current and previous subgradients.
         d, &         ! Direction vector.
         diag, &      ! Diagonal matrix.
         cdiag, &     ! "Convex" diagonal matrix.
         ncdiag, &    ! "Concave" diagonal matrix.
         diagtmp      ! Auxiliary diagonal matrix.
    REAL(KIND=prec), DIMENSION(n*mcu) :: &
         cssmatrix, & ! Matrix whose columns are stored differences of 
                      ! variables^2 when alfn >= 0.
         ncssmatrix, & ! Matrix whose columns are stored differences of 
                      ! variables^2 when alfn < 0.
         csumatrix, & ! Matrix whose columns are stored differences of 
                      ! variables and subgradients when alfn >= 0.
         ncsumatrix   ! Matrix whose columns are stored differences of 
                      ! variables and subgradients when alfn < 0.
    REAL(KIND=prec), DIMENSION(mnma) :: &
         fold         ! Old function values.


! Local Scalars
    REAL(KIND=prec) :: &
         alfn, &      ! Linearization error.
         alfv, &      ! Aggregate linearization error.
         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.
         t, &         ! Stepsize.
         fo, &        ! Previous value of the objective.
         stu, &       ! stu = trans(s)*u. 
         utu, &       ! utu = trans(u)*u. 
         sts, &       ! sts = trans(s)*s. 
         ptmp, &      ! Scalar for convex combination.
         a,b,c
    INTEGER :: i,j, &
         mccc, &      ! Current number of stored "convex" corrections.
         mccnc, &     ! Current number of stored "concave" corrections.
         inewc, &     ! Index for the "convex" circular arrays.
         inewnc, &    ! Index for the "concave" circular arrays.
         inewnma, &   ! Index for the circular arrays containing function values.
         iters, &     ! Null step indicator.
                      !   0  - Null step.
                      !   1  - Serious step.
         nnk, &       ! Consecutive null steps counter.
         nser, &      ! Serious steps counter.
         cnull, &     ! Convex null steps counter.
         ncnull, &    ! Concave null steps counter.
         neps, &      ! Number of consecutive equal stopping criterions.
         ntesf, &     ! Number of tests on function decrease.
         ncres, &     ! Number of restarts.
         nres, &      ! Number of consecutive restarts.
         nout         ! Auxilary printout specification.

! Intrinsic Functions
    INTRINSIC ABS,MAX,SQRT,DOT_PRODUCT

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

    inewnma = 1
    nser    = 0
    cnull   = 0
    ncnull  = 0
    nout    = 0
    nit     = 0
    nfe     = 0
    nge     = 0
    ntesf   = 0
    nres    = 1
    ncres   = -1
    neps    = 0
    iterm   = 0
    iters   = 1
    nnk     = 0
    alfn    = zero
    alfv    = zero
    stu     = one
    sts     = one
    utu     = one
    a       = one
    fold    = -large
      
    xnorm   = large

    epsr   = 0.25_prec+small
    IF (epsl+epsl >= epsr) THEN
       epsr = epsl+epsl + small
       IF (epsr >= half) THEN
          CALL wprint(iterm,-2)
       END IF
    END IF
            
     
! Computation of the value and the subgradient of the objective
! function and the search direction for the first iteration

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

    IF (iterm /= 0) GO TO 900
    CALL restar(n,mccc,inewc,mccnc,inewnc,iters,gp,g,nnk,alfv,alfn,d,diag,cdiag,ncdiag,ncres)
     
! Start of the iteration
            
    iteration: DO


! Computational time

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


! Computation of norms

       IF (iters > 0) THEN
          gnorm = 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,-3)
             
             CALL restar(n,mccc,inewc,mccnc,inewnc,iters,gp,g,nnk,alfv,alfn,d,diag,cdiag,ncdiag,ncres)

             IF (ncres > maxnrs) THEN
                nout = maxnrs
                iterm = -2
                EXIT iteration
             END IF
             
             CYCLE iteration
          END IF
          nout = -1
          iterm = -1
          EXIT iteration
       END IF
       
       
! Stopping criterion

       nit = nit + 1
       pxnorm = xnorm


       IF (iters > 0) THEN
         CALL vxdiag(n,cdiag,g,diag)  !
         xnorm = DOT_PRODUCT(g,diag) + 2.0_prec*alfv

       ELSE
         CALL vxdiag(n,cdiag,ga,diag)  !
         xnorm = DOT_PRODUCT(ga,diag) + 2.0_prec*alfv

       END IF

     
! Tests for termination

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


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

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

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

          ELSE
             neps = 0
          END IF

       ELSE
          neps = 0
       END IF


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

       CALL rprint(nit,nser,cnull,ncnull,nfe,nge,x,f,xnorm,iterm)
      
     
!     Preparation of line search

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


       t = one 

! Line search

       IF (inma == 0) THEN
         CALL armijo(x,g,d,xo,fo,f,p,alfn,xnorm,epsr,iters,nfe,nge,iterm)
       ELSE
         CALL nmarmijo(x,g,d,xo,fo,f,fold,p,alfn,xnorm,epsr,iters,nfe,nge,inewnma,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
      


! Computation of variables and gradients differences

       CALL xdiffy(n,x,xo,s)
       CALL xdiffy(n,g,gp,u)
 

! Computation of aggregate values

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

          IF (nnk == 1) THEN
            CALL aggre1(g,gp,ga,u,cdiag,alfn,alfv,dnorm)
             
          ELSE
            CALL aggre2(g,gp,ga,cdiag,alfn,alfv,dnorm)

          END IF
          
          CALL copy(n,xo,x)
          f = fo
          
       ELSE
          nnk = 0
       END IF

     

      
! Direction finding


       IF (iters > 0) THEN  ! Serious step update and "convex" direction determination
          nser = nser + 1
       
! Serious step initialization

          alfv = zero

! Herskovits diagonal update

          IF (alfn >= zero) THEN
            CALL update(mccc,inewc,s,u,cssmatrix,csumatrix)  ! convex update with sts and stu
            DO i = 1, n
              stu = zero
              sts = zero
              DO j = 1, mccc
                stu = stu + csumatrix((j-1)*n+i)
                sts = sts + cssmatrix((j-1)*n+i)
              END DO
              IF (stu > small) THEN
                cdiag(i) = sts/stu
                IF (cdiag(i) < mintold) cdiag(i) = mintold
                IF (cdiag(i) > told) cdiag(i) = told
              ELSE

                 cdiag(i) = told

              END IF
            END DO
          ELSE

            CALL update(mccnc,inewnc,s,u,ncssmatrix,ncsumatrix)  ! concave update
            DO i = 1, n
              stu = zero
              sts = zero
              DO j = 1, mccnc
                stu = stu + ncsumatrix((j-1)*n+i)
                sts = sts + ncssmatrix((j-1)*n+i)
              END DO
              IF (-stu > small) THEN
                ncdiag(i) = sts/stu
                IF (ncdiag(i) > -mintold) ncdiag(i) = -mintold
                IF (ncdiag(i) < -told) ncdiag(i) = -told
              ELSE

                ncdiag(i) = -told

              END IF
            END DO

          END IF

! Compute "convex" direction (we have a serious step, no need for "concave" direction).
          CALL vxdiag(n,cdiag,g,d)  !
          d = -one*d



       ELSE  ! Update and direction determination for null steps.
     

! Herskovits's diagonal update

         IF (nnk == 1) THEN  ! the first null step

           IF (alfn >= zero) THEN
             CALL update(mccc,inewc,s,u,cssmatrix,csumatrix)  ! convex update with sts and stu
             DO i = 1, n
               stu = zero
               sts = zero
               DO j = 1, mccc
                 stu = stu + csumatrix((j-1)*n+i)
                 sts = sts + cssmatrix((j-1)*n+i)
               END DO
               IF (stu > small) THEN
                 cdiag(i) = sts/stu
                 IF (cdiag(i) < mintold) cdiag(i) = mintold
                 IF (cdiag(i) > told) cdiag(i) = told
               ELSE

                 cdiag(i) = told

               END IF
             END DO
           ELSE


             CALL update(mccnc,inewnc,s,u,ncssmatrix,ncsumatrix)  ! concave update

             DO i = 1, n
               stu = zero
               sts = zero
               DO j = 1, mccnc
                 stu = stu + ncsumatrix((j-1)*n+i)
                 sts = sts + ncssmatrix((j-1)*n+i)
               END DO
               IF (-stu > small) THEN
                 ncdiag(i) = sts/stu
                 IF (ncdiag(i) > -mintold) ncdiag(i) = -mintold
                 IF (ncdiag(i) < -told) ncdiag(i) = -told
               ELSE

                 ncdiag(i) = -told

               END IF
             END DO
           END IF


         ELSE ! consecutive null step

           IF (alfn >= zero) THEN
           ! Skip the update or compute

! If you skip comment from here ... (preliminary results say skip always)
!                   DO i = 1, n
!                      stu = zero
!                      sts = zero
!                      DO j = 1, inewc-1
!                         stu = stu + csumatrix((j-1)*n+i)
!                         sts = sts + cssmatrix((j-1)*n+i)
!                      END DO
!                      DO j = inewc+1, mccc
!                         stu = stu + csumatrix((j-1)*n+i)
!                         sts = sts + cssmatrix((j-1)*n+i)
!                      END DO
!                      stu = stu + s(i)*u(i)
!                      sts = sts + s(i)*s(i)
!
!                    IF (stu > small) THEN
!                         diagtmp(i) = sts/stu
!                         IF (diagtmp(i) < mintold) diagtmp(i) = mintold
!                         IF (diagtmp(i) > told) THEN
!                            diagtmp(i) = told
!                         END IF
!                      ELSE
!                         diagtmp(i) = told
!
!                      END IF
!
!                   END DO
!                   a=zero
!                   b=zero
!                   DO i=1,n
!                     c=ga(i)*ga(i)
!                     a=a+c*diagtmp(i)
!                     b=b+c*cdiag(i)
!                   END DO
!                   IF (a <= b) THEN
!                     CALL update(mccc,inewc,s,u,cssmatrix,csumatrix)  ! convex update with sts and stu
!                     cdiag=diagtmp
!                   END IF
! ... to here
           ELSE

                
             CALL update(mccnc,inewnc,s,u,ncssmatrix,ncsumatrix)  ! concave update

             DO i = 1, n
               stu = zero
               sts = zero
               DO j = 1, mccnc
                 stu = stu + ncsumatrix((j-1)*n+i)
                 sts = sts + ncssmatrix((j-1)*n+i)
               END DO
               IF (-stu > small) THEN
                 ncdiag(i) = sts/stu
                 IF (ncdiag(i) > -mintold) ncdiag(i) = -mintold
                 IF (ncdiag(i) < -told) ncdiag(i) = -told
               ELSE
                         
                 ncdiag(i) = -told

               END IF
             END DO
           END IF

         END IF
             
         IF (alfn >= zero) THEN !convex null step
           cnull = cnull+1

           CALL vxdiag(n,cdiag,ga,d)  !
           d = -one*d

         ELSE !concave null step

           ncnull = ncnull+1

           ptmp = zero
             DO I=1,n

               IF (ABS(ncdiag(i)-cdiag(i)) > small) THEN
                 ptmp = MAX(ncdiag(i)/(ncdiag(i)-cdiag(i)),ptmp)
               END IF
                   
             END DO

             IF (ptmp > one) ptmp = one
                
             DO I=1,n
               diag(i) = ptmp*cdiag(i)+(one-ptmp)*ncdiag(i)
             END DO

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

           END IF

       END IF


    END DO iteration
      

900 CONTINUE
      

! Printout the final results

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

    CALL wprint(iterm,nout)
    CALL rprint(nit,nser,cnull,ncnull,nfe,nge,x,f,xnorm,iterm)

  CONTAINS


    SUBROUTINE restar(n,mccc,inewc,mccnc,inewnc,iters,gp,g,nnk,alfv,alfn,d,diag,cdiag,ncdiag,ncres)  ! Initialization
      
!      USE param, ONLY : zero,one  ! given in host
!      USE subpro, ONLY : copy     ! given in host
      IMPLICIT NONE

! Array Arguments
      REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
           gp        ! Basic subgradient of the objective function.
      REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
           g         ! Current (auxiliary) subgradient of the objective function.
      REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: &
           d, &      ! Search direction.
           diag, &   ! diag = I.
           cdiag, &  ! cdiag = I.
           ncdiag    ! ncdiag = -I.

! Scalar Arguments
      INTEGER, INTENT(IN) :: &
           n         ! Number of variables (given in host)
      INTEGER, INTENT(OUT) :: &
           mccc, &   !
           mccnc, &  !
           inewc, &  !
           inewnc, & !
           nnk       ! Consecutive null steps counter.
      INTEGER, INTENT(INOUT) :: &
           iters, &  ! Null step indicator.
                     !   0  - Null step.
                     !   1  - Serious step.
           ncres     ! Number of restarts.
      REAL(KIND=prec), INTENT(OUT) :: &
           alfn, &   ! Linearization error.
           alfv      ! Aggregate linearization error.


! Restart
      mccc   = 0
      mccnc  = 0
      inewc  = 1
      inewnc = 1
      ncres = ncres + 1
        
      IF (iters == 0) THEN
         CALL copy(n,gp,g)
         iters = 1
         nnk = 0
         alfv=zero
         alfn=zero
      END IF

      d = -one*g
      diag = one
      cdiag = one
      ncdiag = -one
      
    END SUBROUTINE restar


    SUBROUTINE armijo(x,g,d,xo,fo,f,p,alfn,wk,epsr,iters,nfe,nge,iterm) ! Armijo line search

!      USE param, ONLY : one,zero  ! given in host
      USE initializat, 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.
         sigma     ! Tolerance for alfn.
!      USE obj_fun, ONLY : &
!         myf,myg   ! Computation of the value and the subgradient of
!                   ! the objective function.      ! 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.
         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, &          ! Auxiliary values of the objective function
         a_orig, &      ! Armijo parameter.
         inv_a_orig, &  ! Inverse of Armijo parameter.
         alpha, &       ! Armijo step size.
         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  = 2.0_prec
      a_orig = alpha
      inv_a_orig = one/a_orig

! 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=1

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

            p = DOT_PRODUCT(g,d)
            alfn = fo-f+p ! Linearization error
            IF (alfn > -sigma .AND. alfn < zero) alfn=zero

            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 = fo-f+p ! Linearization error
               IF (alfn > -sigma .AND. alfn < zero) alfn=zero

               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 <= 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 = fo-f+p ! Linearization error
               IF (alfn > -sigma .AND. alfn < zero) alfn=zero

               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 = fo-f+p ! Linearization error

       IF (alfn > -sigma .AND. alfn < zero) alfn=zero

!       IF (p-alfn < -epsrwk) PRINT*,'not a good null step',p-alfn, -epsrwk
!       IF (p-max(abs(alfn),eta*alpha*alpha*DOT_PRODUCT(d,d)) < -epsrwk) PRINT*,'not a good null step',p-max(abs(alfn),eta*alpha*alpha*DOT_PRODUCT(d,d)), -epsrwk

    END SUBROUTINE armijo


    SUBROUTINE nmarmijo(x,g,d,xo,fo,f,fold,p,alfn,wk,epsr,iters,nfe,nge,inewnma,iterm) ! Nonmonotone Armijo line search.

!      USE param, ONLY : one,zero
      USE initializat, ONLY : &
!         n, &      ! Number of variables.         ! given in host
!         epsl, &   ! Descent parameter.           ! given in host
!         eta, &    ! Distance measure parameter.  ! given in host
!         mnma, &   ! Maximum number of function values used in line search.  ! given in host
         maxnin, & ! Maximum number of interpolations.
         sigma     ! Tolerance for alfn.
!    USE obj_fun, ONLY : &
!         myf,myg   ! Computation of the value and the subgradient of
                   ! the objective function.       ! 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.
         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) :: &
!         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, &          ! Auxiliary values of the objective function
         a_orig, &      ! Armijo parameter.
         inv_a_orig, &  ! Inverse of Armijo parameter.
         alpha, &       ! Armijo step size.
         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

    alpha  = 2.0_prec
    a_orig = alpha
    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=1

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

             p = DOT_PRODUCT(g,d)
             alfn = fo-f+p ! Linearization error
             IF (alfn > -sigma .AND. alfn < zero) alfn=zero

             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 = fo-f+p ! Linearization error
                IF (alfn > -sigma .AND. alfn < zero) alfn=zero

                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 = fo-f+p ! Linearization error
                IF (alfn > -sigma .AND. alfn < zero) alfn=zero

                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 = fo-f+p ! Linearization error

       IF (alfn > -sigma .AND. alfn < zero) alfn=zero

!       IF (p-alfn < -epsrwk) PRINT*,'not a good null step',p-alfn, -epsrwk
!       IF (p-max(abs(alfn),eta*alpha*alpha*DOT_PRODUCT(d,d)) < -epsrwk) PRINT*,'not a good null step',p-max(abs(alfn),eta*alpha*alpha*DOT_PRODUCT(d,d)), -epsrwk

    END SUBROUTINE nmarmijo


    SUBROUTINE update(mcc,inew,s,u,ssmatrix,sumatrix) ! Matrix update for SMDB
      
!      USE subpro, ONLY : &
!         copy2     ! Copying of two vectors.  ! Given in host
!      USE initializat, ONLY : &
!         n, &      ! Number of variables.     ! Given in host
!         mcu       ! Upper limit for maximum number of stored corrections.  ! Given in host
      IMPLICIT NONE

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

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

! Local Scalars
      INTEGER :: i

     
! Compute new values

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

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


    SUBROUTINE aggre1(g,gp,ga,u,diag,alfn,alfv,dnorm)  ! Computation of aggregate values after first null step.

!      USE param, ONLY : zero,half,one     ! Given in host
!      USE subpro, ONLY : &                ! Given in host
!         vxdiag    ! Vector multiplied by a diagonal matrix.
!      USE initializat, ONLY : &           ! Given in host
!         n, &      ! Number of variables.
!         eta       ! Distance measure parameter, eta >= 0.

      IMPLICIT NONE

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

! Scalar Arguments
      REAL(KIND=prec), INTENT(OUT) :: &
         alfv      ! Aggregate locality measure.
      REAL(KIND=prec), INTENT(IN) :: &
         dnorm     ! sqrt(norm(d)).
      REAL(KIND=prec), INTENT(INOUT) :: &
         alfn      ! Locality measure.

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

! Local Scalars
      REAL(KIND=prec) :: &
         tmp, &    ! temporary storing of alfn.
         p, &      ! trans(gp)*diag*u - alfn.
         q, &      ! q = trans(u)*diag*u, where diag is the diagonal inverse
                   ! approximation of the Hessian.
         lam       ! Multiplier used to calculate aggregate values.
      INTEGER :: i

! Intrinsic Functions
      INTRINSIC ABS,MAX,MIN,SIGN,DOT_PRODUCT

      tmp=alfn
      alfn=MAX(ABS(alfn), eta*dnorm*dnorm)

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

      CALL vxdiag(n,diag,u,tmpn)
      p = - DOT_PRODUCT(tmpn,gp) - alfn
      q = DOT_PRODUCT(u,tmpn)

    
      lam = half + SIGN(half,p)

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

! Computation of the aggregate values

      p = one - lam
      DO i=1,n
        ga(i)=lam*g(i) + p*gp(i)
      END DO
      
      alfv = lam*alfn

      alfn=tmp
      
    END SUBROUTINE aggre1


    SUBROUTINE aggre2(g,gp,ga,diag,alfn,alfv,dnorm)  ! Computation of aggregate values after consequtive null steps.
      
!     USE param, ONLY : zero,one,small    ! Given in host
!     USE subpro, ONLY : &                 ! Given in host
!         vxdiag, &  ! Vector multiplied by a diagonal matrix.
!         xdiffy     ! Difference of two vectors.
!     USE initializat, ONLY : &            ! Given in host
!         n, &       ! Number of variables.
!         eta        ! Distance measure parameter, eta >= 0.
      IMPLICIT NONE

! Array Arguments
      REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: &
         diag, &    ! Diagonal variable metric update.
         g, &       ! Current (auxiliary) subgradient of the objective function.
         gp         ! Previous subgradient of the objective function.
      REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: &
         ga         ! Aggregate subgradient of the objective function.

! Scalar Arguments
      REAL(KIND=prec), INTENT(INOUT) :: &
         alfn, &    ! Locality measure.
         alfv       ! Aggregate locality measure.
      REAL(KIND=prec), INTENT(IN) :: &
         dnorm      ! sqrt(norm(d)).
      
! Local arrays
      REAL(KIND=prec), DIMENSION(n) :: tmpn2,tmpn3,tmpn4

! Local Scalars
      REAL(KIND=prec) :: &
         tmp, &    ! temporary storage for alfn.
         pr, &     ! pr = trans(gp-ga) diag (gp-ga).
         rrp, &    ! rrp = trans(gp-ga) diag ga - alfv.
         prqr, &   ! prqr = trans(gp-ga) diag (g-ga).
         rrq, &    ! rrq = trans(g-ga) diag ga - alfv + alfn.
         qr, &     ! qr = trans(g-ga) diag (g-ga).
         pq, &     ! pq = trans(g-gp) diag (g-gp).
         qqp, &    ! qqp = trans(g-gp) diag g + alfn.
         lam1, &   ! Multiplier used to calculate aggregate values.
         lam2, &   ! Multiplier used to calculate aggregate values.
         tmp1, &   ! Auxiliary scalar.
         tmp2, &   ! Auxiliary scalar.
         tmp3      ! Auxiliary scalar.
      INTEGER :: i

      
! Intrinsic Functions
      INTRINSIC ABS,MIN,MAX,DOT_PRODUCT

      tmp=alfn
      alfn=MAX(ABS(alfn), eta*dnorm*dnorm)

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

      CALL vxdiag(n,diag,tmpn2,tmpn3)

         
      pr = DOT_PRODUCT(tmpn3,tmpn2)
      rrp = DOT_PRODUCT(tmpn3,ga)
      CALL xdiffy(n,g,ga,tmpn4)
      prqr = DOT_PRODUCT(tmpn3,tmpn4)


! calculation of qr = trans(g - ga) diag (g - ga) and rrq = trans(g-ga) diag ga + alfn - alfv.
       
      CALL vxdiag(n,diag,tmpn4,tmpn3)
      rrq = DOT_PRODUCT(tmpn3,ga)
      qr = DOT_PRODUCT(tmpn4,tmpn3)
      
      pq = qr - prqr - prqr + pr
      qqp = pq + prqr + rrq - pr - rrp + alfn
      rrp = rrp - alfv
      rrq = rrq + alfn - alfv

! computation of multipliers lam1 and lam2

      IF (pr > zero .AND. qr > zero) THEN
        tmp1 = rrq/qr
        tmp2 = prqr/qr
        tmp3 = pr - prqr*tmp2
        IF (tmp3 /= zero) THEN
          lam1 = (tmp1*prqr - rrp)/tmp3
          lam2 = -tmp1 - lam1*tmp2
          IF (lam1*(lam1 - one) < zero .AND. &
               lam2*(lam1 + lam2 - one) < zero) THEN
            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

            alfn=tmp

            RETURN

          END IF

        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))
      tmp3 = (lam2*qr + rrq+rrq)*lam2
      tmp1 = zero
      IF (alfv >= zero) tmp1 = one
      IF (pr > zero) tmp1 = MIN(one,MAX(zero,-rrp/pr))
      tmp2 = (tmp1*pr + rrp+rrp)*tmp1
      IF (tmp2 < tmp3) THEN
        tmp3 = tmp2
        lam1 = tmp1
        lam2 = zero
      END IF

      IF (qqp*(qqp - pq) < zero) THEN
        IF (qr + rrq + rrq - qqp*qqp/pq < tmp3) THEN
          lam1 = qqp/pq
          lam2 = one - lam1
        END IF
      END IF


      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

      alfn=tmp

    
    END SUBROUTINE aggre2
      
  END SUBROUTINE smdb

!************************************************************************
!*
!*     * SUBROUTINE wprint *
!*
!*     Printout the warning and error messages.
!*
!************************************************************************
      
  SUBROUTINE wprint(iterm,nout)
    USE initializat, ONLY : &
         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
    IMPLICIT NONE

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


    IF (iprint >= 0) THEN

! Initial error messages

       IF (iterm == -5) THEN
          IF (nout == 1) WRITE (6,FMT='(1X,''Error: '' &
               ''Number of variables (n) is too small, iterm='',I3)') iterm
          IF (nout == 2) WRITE (6,FMT='(1X,''Error: '' &
               ''The maximum number of stored corrections (mcu) '' &
               ''is too small, iterm='',I3)') iterm
          IF (nout == 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 == -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.'')')
             RETURN
          END IF
         

! Printout the final results
            
          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 == 4) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''Number of function evaluations > '',I5)') nout
          IF (iterm == 5) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''Number of iterations > '',I5)') nout
          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 == 8) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''Failure in attaining the demanded accuracy.'')')
          IF (iterm == -1) WRITE (6,FMT='(1X,''Abnormal exit: '' &
               ''Two consecutive restarts.'')')
          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.'')')
       END IF
    END IF
      
  END SUBROUTINE wprint

            
!************************************************************************
!*
!*     * SUBROUTINE rprint *
!*      
!*     Printout the (final) results.
!*
!************************************************************************
      
  SUBROUTINE rprint(nit,nser,cnull,ncnull,nfe,nge,x,f,wk,iterm)
    USE initializat, ONLY : &
         n, &         ! Number of variables
         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
    IMPLICIT NONE

! Scalar Arguments
    INTEGER, INTENT(IN) :: & 
         nit, &       ! Number of used iterations.
         nser, &      ! Number of serious steps.
         cnull, &     ! Number of convex null steps.
         ncnull, &    ! Number of nonconvex null steps.
         nfe, &       ! Number of used function evaluations.
         nge, &       ! Number of used subgradient evaluations.
         iterm        ! Cause of termination:
                      !   1  - The problem has been solved with desired accuracy.
                      !   2  - Changes in function values < tolf in mtesf
                      !        subsequent iterations.
                      !   3  - Changes in function value < tolf*small*MAX(|f_k|,|f_(k-1)|,1),
                      !        where small is the smallest positive number such that 
                      !        1.0 + small > 1.0.
                      !   4  - Number of function calls > mfe.
                      !   5  - Number of iterations > mit.
                      !   6  - Time limit exceeded. 
                      !   7  - f < tolb.
                      !   8  - Failure in attaining the demanded accuracy.
                      !  -1  - Two consecutive restarts.
                      !  -2  - Number of restarts > maximum number of restarts.
                      !  -3  - Failure in function or subgradient calculations 
                      !        (assigned by the user).
                      !  -4  
                      !  -5  - Invalid input parameters.


    REAL(KIND=prec), INTENT(IN) :: &
         f, &         ! Value of the objective function.
         wk           ! Value of the 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 &
            )') nit,nfe,nge,f,wk
       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,''nser='',I5,2X,''cnull='',I5,2X, &
         ''ncnull='',I5,2X,''nfe='',I5,2X,''nge='',I5,2X,''f='',D15.8,2X,''wk='',D11.4,2X, &
         ''iterm='',I3)') nit,nser,cnull,ncnull,nfe,nge,f,wk,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 smdb_mod
