C  THIS VERSION: CUTEST 2.2 - 2023-11-29 AT 15:00 GMT.

#include "cutest_modules.h"
#include "cutest_routines.h"

      PROGRAM LBFGSB_main

      USE CUTEST_KINDS_precision
      IMPLICIT NONE

C
C  LBFGSB test driver for problems derived from SIF files.
C
C  Nick Gould, for CGT Productions.
C  September 2004
C  Revised for CUTEst, January 2013
C
      INTEGER ( KIND = ip_ )  I, out, N, M, INPUT, MAXIT, status
      INTEGER ( KIND = ip_ )  IFLAG , INSPEC, IPRINT, lwa, liwa
      INTEGER ( KIND = ip_ )  ISAVE( 44 )
      INTEGER ( KIND = ip_ ) :: io_buffer = 11
      REAL ( KIND = rp_ ) F, GNORM
      REAL ( KIND = rp_ ) PGTOL, FACTR, INFTY, DSAVE( 29 )
      CHARACTER ( LEN = 60 ) :: TASK, CSAVE
      LOGICAL LSAVE( 4 )
      PARAMETER  ( out  = 6 )
      INTEGER ( KIND = ip_ ), ALLOCATABLE, DIMENSION( : ) :: NBD, IWA
      REAL ( KIND = rp_ ), ALLOCATABLE, DIMENSION( : ) :: X, XL, XU
      REAL ( KIND = rp_ ), ALLOCATABLE, DIMENSION( : ) :: G, WA
      PARAMETER ( INPUT = 55, INSPEC = 56 )
      PARAMETER ( INFTY = 1.0E+19_rp_ )
      CHARACTER ( LEN = 10 ) :: PNAME, SPCDAT
      CHARACTER ( LEN = 10 ), ALLOCATABLE, DIMENSION( : )  :: XNAMES
      REAL ( KIND = rp_ ) :: CPU( 4 ), CALLS( 4 )
      EXTERNAL SETULB
C
C  Open the Spec file for the method.
C
      SPCDAT = 'LBFGSB.SPC'
      OPEN ( INSPEC, FILE = SPCDAT, FORM = 'FORMATTED',
     *      STATUS = 'OLD' )
      REWIND INSPEC
C
C  Read input Spec data.
C
C     M        : the maximum number of variable metric corrections
C     MAXIT    : the maximum number of iterations,
C     IPRINT   : print level (<0,none,=0,one line/iteration,>1,more detail)
C     FACTR    : the function accuracy tolerence (see hint below)
C     PGTOL    : the required norm of the projected gradient
C
C  Hint - the iteration will stop when
C
C         (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
C
C   where epsmch is the machine precision, which is automatically
C   generated by the code. Typical values for factr: 1.d+12 for
C   low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely
C   high accuracy.

      READ ( INSPEC, 1000 ) M, MAXIT, IPRINT, FACTR, PGTOL
C
C  Close input file.
C
      CLOSE ( INSPEC )
C
C  Open the relevant file.
C
      OPEN ( INPUT, FILE = 'OUTSDIF.d', FORM = 'FORMATTED',
     *       STATUS = 'OLD' )
C
C  Check to see if there is sufficient room
C
      CALL CUTEST_udimen_r( status, INPUT, N )
      IF ( status /= 0 ) GO TO 910

      liwa = 3 * n
!     lwa  = 2 * m * n + 4 * n + 12 * m * m + 12 * m ! Version 2
      lwa  = 2 * m * n + 5 * n + 11 * m * m + 8 * m  ! Version 3
      ALLOCATE( NBD( n ), IWA( liwa ), X( n ), XL( n ), XU( n ), G( n ),
     +          WA( lwa ), XNAMES( n ), STAT = status )
      IF ( status /= 0 ) GO TO 990
C
C  Set up SIF data.
C
      CALL CUTEST_usetup_r( status, INPUT, out, io_buffer, 
     *                     N, X, XL, XU )
      IF ( status /= 0 ) GO TO 910
C
C  Set bound constraint status
C
      DO 10 I = 1, N
         IF( XL( I ) .LE. - INFTY ) THEN
            IF( XU( I ) .GE. INFTY ) THEN
               NBD( I ) = 0
            ELSE
               NBD( I ) = 3
            END IF
         ELSE
            IF( XU( I ) .GE. INFTY ) THEN
               NBD( I ) = 1
            ELSE
               NBD( I ) = 2
            END IF
         END IF
   10 CONTINUE
C
C  Obtain variable names.
C
      CALL CUTEST_unames_r( status, N, PNAME, XNAMES )
      IF ( status /= 0 ) GO TO 910
C
C  Set up algorithmic input data.
C
      IFLAG = 0
C
C  Optimization loo
C
      TASK = 'START'
   30 CONTINUE
C
C  Call the optimizer
C
         CALL SETULB( N, M, X, XL, XU, NBD, F, G, FACTR, PGTOL, WA,
     +                IWA, TASK, IPRINT, CSAVE, LSAVE, ISAVE, DSAVE )
C
C  Evaluate the function, f, and gradient, G
C
         IF (TASK( 1: 2 ) .EQ. 'FG' ) THEN
            CALL CUTEST_uofg_r( status, N, X, F, G, .TRUE. )
            IF ( status /= 0 ) GO TO 910
            GO TO 30
C
C  Test for convergence
C
         ELSE IF ( TASK( 1: 4 ) .EQ. 'CONV' ) THEN
            IFLAG = 0
         ELSE IF (  TASK( 1: 4 ) .EQ. 'ABNO' ) THEN
            IFLAG = 1
            WRITE( out, "( ' Abnormal exit ' )" )
         ELSE IF (  TASK( 1: 5 ) .EQ. 'ERROR' ) THEN
            IFLAG = 2
            WRITE( out, "( ' Error exit ' )" )
         ELSE IF (  TASK( 1: 5 ) .EQ. 'NEW_X' ) THEN
            IF ( ISAVE( 30 ) .GT. MAXIT ) THEN
              IFLAG = 3
              WRITE( out,
     *           "( ' Maximum number of iterations exceeded ' )" )
            ELSE
               GO TO 30
            END IF
         END IF
C
C  Terminal exit.
C
      CALL CUTEST_ureport_r( status, CALLS, CPU )
      IF ( status /= 0 ) GO TO 910
      GNORM = DSAVE( 13 )
      WRITE ( out, 2010 ) F, GNORM
      DO 120 I = 1, N
         WRITE( out, 2020 ) XNAMES( I ), XL( I ), X( I ), XU( I ),
     *                       G( I )
  120 CONTINUE
      WRITE ( out, 2000 ) PNAME, N, INT( CALLS(1) ), INT( CALLS(2) ),
     *                    IFLAG, F, CPU(1), CPU(2)
      CLOSE( INPUT  )
      CALL CUTEST_uterminate_r( status )
      STOP

  910 CONTINUE
      WRITE( out, "( ' CUTEst error, status = ', i0, ', stopping' )")
     *   status
      STOP

  990 CONTINUE
      WRITE( out, "( ' Allocation error, status = ', I0 )" ) status
      STOP
C
C  Non-executable statements.
C
 1000 FORMAT( 3( I10, / ), ( D10.3, / ), D10.3 )
 2000 FORMAT( /, 24('*'), ' CUTEst statistics ', 24('*') //
     *    ,' Package used            :  L-BFGS-B',     /
     *    ,' Problem                 :  ', A10,    /
     *    ,' # variables             =      ', I10 /
     *    ,' # objective functions   =      ', I10 /
     *    ,' # objective gradients   =      ', I10 /
     *     ' Exit code               =      ', I10 /
     *    ,' Final f                 = ', E15.7 /
     *    ,' Set up time             =      ', 0P, F10.2, ' seconds' /
     *     ' Solve time              =      ', 0P, F10.2, ' seconds' //
     *     66('*') / )
 2010 FORMAT( ' Final objective function value   = ', 1P, D12.4,
     *        /, ' Final norm of projected gradient = ', 1P, D12.4,
     *        //, '                XL           X        XU',
     *           '           G ' )
 2020 FORMAT(  1X, A10, 1P, 4D12.4 )
      END

      SUBROUTINE REORDA( NC, NNZ, IRN, JCN, A, IP, IW )
      USE CUTEST_KINDS_precision
      INTEGER ( KIND = ip_ ) NC, NNZ
      INTEGER ( KIND = ip_ ) IRN( NNZ  ), JCN( NNZ )
      INTEGER ( KIND = ip_ ) IW( NC + 1 ), IP( NC + 1 )
      REAL ( KIND = rp_ )  A( NNZ )

C  Sort a sparse matrix from arbitrary order to column order

C  Nick Gould
C  7th November, 1990

      INTEGER ( KIND = ip_ ) I, J, K, L, IC, NCP1, ITEMP, JTEMP,  LOCAT
      REAL ( KIND = rp_ ) ANEXT , ATEMP

C  Initialize the workspace as zero

      NCP1       = NC + 1
      DO 10 J    = 1, NCP1
         IW( J ) = 0
   10 CONTINUE

C  Pass 1. Count the number of elements in each column

      DO 20 K   = 1, NNZ
        J       = JCN( K )
        IW( J ) = IW( J ) + 1
   20 CONTINUE

C  Put the positions where each column begins in
C  a compressed collection into IP and IW

      IP( 1 )       = 1
      DO 30 J       = 2, NCP1
        IP( J )     = IW( J - 1 ) + IP( J - 1 )
        IW( J - 1 ) = IP( J - 1 )
   30 CONTINUE

C  Pass 2. Reorder the elements into column order.
C          Fill in each column in turn

      DO 70 IC = 1, NC

C  Consider the next unfilled position in column IC

        DO 60 K = IW( IC ), IP( IC + 1 ) - 1

C  The entry should be placed in column J

          I       = IRN( K )
          J       = JCN( K )
          ANEXT   = A( K )
          DO 40 L = 1, NNZ

C  See if the entry is already in place

             IF ( J .EQ. IC ) GO TO 50
             LOCAT = IW( J )

C  As a new entry is placed in column J, increase the pointer
C  IW( J ) by one

             IW( J  ) = LOCAT + 1

C  Record details of the entry which currently occupies location LOCAT

             ITEMP = IRN( LOCAT )
             JTEMP = JCN( LOCAT )
             ATEMP = A( LOCAT )

C  Move the new entry to its correct place

             IRN( LOCAT ) = I
             JCN( LOCAT ) = J
             A( LOCAT )   = ANEXT

C  Make the displaced entry the new entry

             I          = ITEMP
             J          = JTEMP
             ANEXT      = ATEMP
   40     CONTINUE

C  Move the new entry to its correct place

   50     CONTINUE
          JCN( K ) = J
          IRN( K ) = I
          A( K )   = ANEXT
   60   CONTINUE
   70 CONTINUE

      RETURN

C  End of REORDA

      END

      SUBROUTINE TIMER( TTIME )
      USE CUTEST_KINDS_precision

C  CPU timer

      REAL ( KIND = rp_ ) TTIME

      CALL CPU_TIME( TTIME )
      RETURN
      END
