LCOV - code coverage report
Current view: top level - src/motion - cp_lbfgs.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 975 1281 76.1 %
Date: 2024-04-20 06:29:22 Functions: 23 24 95.8 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief LBFGS-B routine (version 3.0, April 25, 2011)
      10             : !> \note
      11             : !>      L-BFGS-B (version 3.0, April 25, 2011) converted to Fortran 90 module
      12             : !> \par History
      13             : !>      02.2005 Update to the new version 2.4 and deleting the blas part of
      14             : !>              the code (Teodoro Laino)
      15             : !>      11.2012 New version 3.0 converted to Fortran 90 (Matthias Krack)
      16             : !>      12.2020 Implementation of Space Group Symmetry (Pierre-André Cazade)
      17             : !> \author Fawzi Mohamed (first version)
      18             : ! **************************************************************************************************
      19             : MODULE cp_lbfgs
      20             :    USE bibliography,                    ONLY: Byrd1995,&
      21             :                                               cite_reference
      22             :    USE cp_files,                        ONLY: open_file
      23             :    USE kinds,                           ONLY: dp
      24             :    USE machine,                         ONLY: m_walltime
      25             :    USE space_groups,                    ONLY: spgr_apply_rotations_coord,&
      26             :                                               spgr_apply_rotations_force
      27             :    USE space_groups_types,              ONLY: spgr_type
      28             : #include "../base/base_uses.f90"
      29             : 
      30             :    IMPLICIT NONE
      31             : 
      32             :    PRIVATE
      33             : 
      34             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_lbfgs'
      35             : 
      36             :    PUBLIC :: setulb
      37             : 
      38             : CONTAINS
      39             : 
      40             : !===========   L-BFGS-B (version 3.0.  April 25, 2011  =================
      41             : !
      42             : !     This is a modified version of L-BFGS-B.
      43             : !
      44             : !     Major changes are described in the accompanying paper:
      45             : !
      46             : !         Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778:
      47             : !         L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constraine
      48             : !         Optimization"  (2011). To appear in  ACM Transactions on
      49             : !         Mathematical Software,
      50             : !
      51             : !     The paper describes an improvement and a correction to Algorithm 7
      52             : !     It is shown that the performance of the algorithm can be improved
      53             : !     significantly by making a relatively simple modication to the subs
      54             : !     minimization phase. The correction concerns an error caused by the
      55             : !     of routine dpmeps to estimate machine precision.
      56             : !
      57             : !     The total work space **wa** required by the new version is
      58             : !
      59             : !                  2*m*n + 11m*m + 5*n + 8*m
      60             : !
      61             : !     the old version required
      62             : !
      63             : !                  2*m*n + 12m*m + 4*n + 12*m
      64             : !
      65             : !
      66             : !            J. Nocedal  Department of Electrical Engineering and
      67             : !                        Computer Science.
      68             : !                        Northwestern University. Evanston, IL. USA
      69             : !
      70             : !
      71             : !           J.L Morales  Departamento de Matematicas,
      72             : !                        Instituto Tecnologico Autonomo de Mexico
      73             : !                        Mexico D.F. Mexico.
      74             : !
      75             : !                        March  2011
      76             : !
      77             : !=======================================================================
      78             : ! **************************************************************************************************
      79             : !> \brief          This subroutine partitions the working arrays wa and iwa, and
      80             : !>                 then uses the limited memory BFGS method to solve the bound
      81             : !>                 constrained optimization problem by calling mainlb.
      82             : !>                 (The direct method will be used in the subspace minimization.)
      83             : !> \param n        n is the dimension of the problem.
      84             : !> \param m        m is the maximum number of variable metric corrections
      85             : !>                 used to define the limited memory matrix.
      86             : !> \param x        On entry x is an approximation to the solution.
      87             : !>                 On exit x is the current approximation.
      88             : !> \param lower_bound  the lower bound on x.
      89             : !> \param upper_bound  the upper bound on x.
      90             : !> \param nbd      nbd represents the type of bounds imposed on the
      91             : !>                 variables, and must be specified as follows:
      92             : !>                 nbd(i)=0 if x(i) is unbounded,
      93             : !>                        1 if x(i) has only a lower bound,
      94             : !>                        2 if x(i) has both lower and upper bounds, and
      95             : !>                        3 if x(i) has only an upper bound.
      96             : !> \param f        On first entry f is unspecified.
      97             : !>                 On final exit f is the value of the function at x.
      98             : !> \param g        On first entry g is unspecified.
      99             : !>                 On final exit g is the value of the gradient at x.
     100             : !> \param factr    factr >= 0 is specified by the user.  The iteration
     101             : !>                 will stop when
     102             : !>
     103             : !>                 (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
     104             : !>
     105             : !>                 where epsmch is the machine precision, which is automatically
     106             : !>                 generated by the code. Typical values for factr: 1.d+12 for
     107             : !>                 low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely
     108             : !>                 high accuracy.
     109             : !> \param pgtol    pgtol >= 0 is specified by the user.  The iteration
     110             : !>                 will stop when
     111             : !>
     112             : !>                 max{|proj g_i | i = 1, ..., n} <= pgtol
     113             : !>
     114             : !>                 where pg_i is the ith component of the projected gradient.
     115             : !> \param wa       working array
     116             : !> \param iwa      integer working array
     117             : !> \param task     is a working string of characters of length 60 indicating
     118             : !>                 the current job when entering and quitting this subroutine.
     119             : !> \param iprint   iprint is a variable that must be set by the user.
     120             : !>                 It controls the frequency and type of output generated:
     121             : !>                 iprint<0    no output is generated;
     122             : !>                 iprint=0    print only one line at the last iteration;
     123             : !>                 0<iprint<99 print also f and |proj g| every iprint iterations;
     124             : !>                 iprint=99   print details of every iteration except n-vectors;
     125             : !>                 iprint=100  print also the changes of active set and final x;
     126             : !>                 iprint>100  print details of every iteration including x and g;
     127             : !>                 When iprint > 0, the file iterate.dat will be created to
     128             : !>                 summarize the iteration.
     129             : !> \param csave    is a working string of characters
     130             : !> \param lsave    lsave is a working array
     131             : !>                 On exit with 'task' = NEW_X, the following information is available:
     132             : !>                 If lsave(1) = .true.  then  the initial X has been replaced by
     133             : !>                               its projection in the feasible set
     134             : !>                 If lsave(2) = .true.  then  the problem is constrained;
     135             : !>                 If lsave(3) = .true.  then  each variable has upper and lower bounds;
     136             : !> \param isave    isave is a working array
     137             : !>                 On exit with 'task' = NEW_X, the following information is available:
     138             : !>                 isave(22) = the total number of intervals explored in the
     139             : !>                         search of Cauchy points;
     140             : !>                 isave(26) = the total number of skipped BFGS updates before the current iteration;
     141             : !>                 isave(30) = the number of current iteration;
     142             : !>                 isave(31) = the total number of BFGS updates prior the current iteration;
     143             : !>                 isave(33) = the number of intervals explored in the search of
     144             : !>                             Cauchy point in the current iteration;
     145             : !>                 isave(34) = the total number of function and gradient evaluations;
     146             : !>                 isave(36) = the number of function value or gradient
     147             : !>                             evaluations in the current iteration;
     148             : !>                 if isave(37) = 0  then the subspace argmin is within the box;
     149             : !>                 if isave(37) = 1  then the subspace argmin is beyond the box;
     150             : !>                 isave(38) = the number of free variables in the current iteration;
     151             : !>                 isave(39) = the number of active constraints in the current iteration;
     152             : !>                 n + 1 - isave(40) = the number of variables leaving the set of
     153             : !>                                     active constraints in the current iteration;
     154             : !>                 isave(41) = the number of variables entering the set of active
     155             : !>                             constraints in the current iteration.
     156             : !> \param dsave    dsave is a working array of dimension 29.
     157             : !>                 On exit with 'task' = NEW_X, the following information is available:
     158             : !>                 dsave(1) = current 'theta' in the BFGS matrix;
     159             : !>                 dsave(2) = f(x) in the previous iteration;
     160             : !>                 dsave(3) = factr*epsmch;
     161             : !>                 dsave(4) = 2-norm of the line search direction vector;
     162             : !>                 dsave(5) = the machine precision epsmch generated by the code;
     163             : !>                 dsave(7) = the accumulated time spent on searching for Cauchy points;
     164             : !>                 dsave(8) = the accumulated time spent on subspace minimization;
     165             : !>                 dsave(9) = the accumulated time spent on line search;
     166             : !>                 dsave(11) = the slope of the line search function at the current point of line search;
     167             : !>                 dsave(12) = the maximum relative step length imposed in line search;
     168             : !>                 dsave(13) = the infinity norm of the projected gradient;
     169             : !>                 dsave(14) = the relative step length in the line search;
     170             : !>                 dsave(15) = the slope of the line search function at the starting point of the line search;
     171             : !>                 dsave(16) = the square of the 2-norm of the line search direction vector.
     172             : !> \param trust_radius ...
     173             : !> \param spgr ...
     174             : !> \par History
     175             : !>      12.2020 Implementation of Space Group Symmetry [pcazade]
     176             : !> \author         NEOS, November 1994. (Latest revision June 1996.)
     177             : !>                 Optimization Technology Center.
     178             : !>                 Argonne National Laboratory and Northwestern University.
     179             : !>                 Written by
     180             : !>                             Ciyou Zhu
     181             : !>                 in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
     182             : ! **************************************************************************************************
     183        4617 :    SUBROUTINE setulb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, wa, iwa, &
     184             :                      task, iprint, csave, lsave, isave, dsave, trust_radius, spgr)
     185             : 
     186             :       INTEGER, INTENT(in)                                :: n, m
     187             :       REAL(KIND=dp), INTENT(inout)                       :: x(n)
     188             :       REAL(KIND=dp)                                      :: lower_bound(n), upper_bound(n)
     189             :       INTEGER                                            :: nbd(n)
     190             :       REAL(KIND=dp)                                      :: f, g(n)
     191             :       REAL(KIND=dp), INTENT(in)                          :: factr, pgtol
     192             :       REAL(KIND=dp)                                      :: wa(2*m*n + 5*n + 11*m*m + 8*m)
     193             :       INTEGER                                            :: iwa(3*n)
     194             :       CHARACTER(LEN=60)                                  :: task
     195             :       INTEGER                                            :: iprint
     196             :       CHARACTER(LEN=60)                                  :: csave
     197             :       LOGICAL                                            :: lsave(4)
     198             :       INTEGER                                            :: isave(44)
     199             :       REAL(KIND=dp)                                      :: dsave(29)
     200             :       REAL(KIND=dp), INTENT(in)                          :: trust_radius
     201             :       TYPE(spgr_type), OPTIONAL, POINTER                 :: spgr
     202             : 
     203             :       INTEGER                                            :: i, ld, lr, lsnd, lss, lsy, lt, lwa, lwn, &
     204             :                                                             lws, lwt, lwy, lxp, lz
     205             : 
     206             : !     References:
     207             : !
     208             : !       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
     209             : !       memory algorithm for bound constrained optimization'',
     210             : !       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
     211             : !
     212             : !       [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
     213             : !       limited memory FORTRAN code for solving bound constrained
     214             : !       optimization problems'', Tech. Report, NAM-11, EECS Department,
     215             : !       Northwestern University, 1994.
     216             : !
     217             : !       (Postscript files of these papers are available via anonymous
     218             : !        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
     219             : !
     220             : !                           *  *  *
     221             : 
     222        4617 :       IF (task == 'START') THEN
     223         192 :          CALL cite_reference(Byrd1995)
     224         192 :          isave(1) = m*n
     225         192 :          isave(2) = m**2
     226         192 :          isave(3) = 4*m**2
     227             :          ! ws      m*n
     228         192 :          isave(4) = 1
     229             :          ! wy      m*n
     230         192 :          isave(5) = isave(4) + isave(1)
     231             :          ! wsy     m**2
     232         192 :          isave(6) = isave(5) + isave(1)
     233             :          ! wss     m**2
     234         192 :          isave(7) = isave(6) + isave(2)
     235             :          ! wt      m**2
     236         192 :          isave(8) = isave(7) + isave(2)
     237             :          ! wn      4*m**2
     238         192 :          isave(9) = isave(8) + isave(2)
     239             :          ! wsnd    4*m**2
     240         192 :          isave(10) = isave(9) + isave(3)
     241             :          ! wz      n
     242         192 :          isave(11) = isave(10) + isave(3)
     243             :          ! wr      n
     244         192 :          isave(12) = isave(11) + n
     245             :          ! wd      n
     246         192 :          isave(13) = isave(12) + n
     247             :          ! wt      n
     248         192 :          isave(14) = isave(13) + n
     249             :          ! wxp     n
     250         192 :          isave(15) = isave(14) + n
     251             :          ! wa      8*m
     252         192 :          isave(16) = isave(15) + n
     253             :       END IF
     254        4617 :       lws = isave(4)
     255        4617 :       lwy = isave(5)
     256        4617 :       lsy = isave(6)
     257        4617 :       lss = isave(7)
     258        4617 :       lwt = isave(8)
     259        4617 :       lwn = isave(9)
     260        4617 :       lsnd = isave(10)
     261        4617 :       lz = isave(11)
     262        4617 :       lr = isave(12)
     263        4617 :       ld = isave(13)
     264        4617 :       lt = isave(14)
     265        4617 :       lxp = isave(15)
     266        4617 :       lwa = isave(16)
     267             : 
     268             :       !in case we use a trust radius we set the boundaries to be one times the trust radius away from the current positions
     269             :       !the original implementation only allowed for boundaries that remain constant during the optimization.
     270             :       !This way of including a trust radius seems to work,
     271             :       !but the change of the boundaries during optimization might introduce some not yet discovered problems.
     272        4617 :       IF (trust_radius >= 0) THEN
     273       31794 :          DO i = 1, n
     274       31752 :             lower_bound(i) = x(i) - trust_radius
     275       31752 :             upper_bound(i) = x(i) + trust_radius
     276       31794 :             nbd(i) = 2
     277             :          END DO
     278             :       END IF
     279             : 
     280             :       ! passes spgr to mainlb
     281             :       CALL mainlb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, &
     282             :                   wa(lws), wa(lwy), wa(lsy), wa(lss), wa(lwt), &
     283             :                   wa(lwn), wa(lsnd), wa(lz), wa(lr), wa(ld), wa(lt), wa(lxp), &
     284             :                   wa(lwa), &
     285             :                   iwa(1), iwa(n + 1), iwa(2*n + 1), task, iprint, &
     286        4617 :                   csave, lsave, isave(22), dsave, spgr=spgr)
     287             : 
     288        4617 :       RETURN
     289             : 
     290             :    END SUBROUTINE setulb
     291             : 
     292             : ! **************************************************************************************************
     293             : !> \brief        This subroutine solves bound constrained optimization problems by
     294             : !>               using the compact formula of the limited memory BFGS updates.
     295             : !> \param n      n is the number of variables
     296             : !> \param m      m is the maximum number of variable metric
     297             : !>               corrections allowed in the limited memory matrix.
     298             : !> \param x      On entry x is an approximation to the solution.
     299             : !>               On exit x is the current approximation.
     300             : !> \param lower_bound  lower_bound is the lower bound of x.
     301             : !> \param upper_bound  upper_bound is the upper bound of x.
     302             : !> \param nbd    nbd represents the type of bounds imposed on the
     303             : !>               variables, and must be specified as follows:
     304             : !>               nbd(i)=0 if x(i) is unbounded,
     305             : !>               1 if x(i) has only a lower bound,
     306             : !>               2 if x(i) has both lower and upper bounds,
     307             : !>               3 if x(i) has only an upper bound.
     308             : !> \param f      On first entry f is unspecified.
     309             : !>               On final exit f is the value of the function at x.
     310             : !> \param g      On first entry g is unspecified.
     311             : !>               On final exit g is the value of the gradient at x.
     312             : !> \param factr  factr >= 0 is specified by the user.  The iteration
     313             : !>               will stop when
     314             : !>
     315             : !>               (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
     316             : !>
     317             : !>               where epsmch is the machine precision, which is automatically
     318             : !>               generated by the code.
     319             : !> \param pgtol  pgtol >= 0 is specified by the user.  The iteration
     320             : !>               will stop when
     321             : !>
     322             : !>                 max{|proj g_i | i = 1, ..., n} <= pgtol
     323             : !>
     324             : !>               where pg_i is the ith component of the projected gradient.
     325             : !> \param ws     ws, wy, sy, and wt are working arrays used to store the following
     326             : !>               information defining the limited memory BFGS matrix:
     327             : !>               ws stores S, the matrix of s-vectors;
     328             : !> \param wy     stores Y, the matrix of y-vectors;
     329             : !> \param sy     stores S'Y;
     330             : !> \param ss     stores S'S;
     331             : !> \param wt     stores the Cholesky factorization of (theta*S'S+LD^(-1)L');
     332             : !>               see eq. (2.26) in [3].
     333             : !> \param wn     wn is a working array of dimension 2m x 2m
     334             : !>               used to store the LEL^T factorization of the indefinite matrix
     335             : !>               K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
     336             : !>                   [L_a -R_z           theta*S'AA'S ]
     337             : !>
     338             : !>               where     E = [-I  0]
     339             : !>                             [ 0  I]
     340             : !> \param snd    is a working array of dimension 2m x 2m
     341             : !>               used to store the lower triangular part of
     342             : !>               N = [Y' ZZ'Y   L_a'+R_z']
     343             : !>                   [L_a +R_z  S'AA'S   ]
     344             : !> \param z      z(n),r(n),d(n),t(n), xp(n),wa(8*m) are working arrays
     345             : !>               z  is used at different times to store the Cauchy point and
     346             : !>               the Newton point.
     347             : !> \param r      working array
     348             : !> \param d      working array
     349             : !> \param t      workign array
     350             : !> \param xp     xp is a workng array used to safeguard the projected Newton direction
     351             : !> \param wa     working array
     352             : !> \param index  In subroutine freev, index is used to store the free and fixed
     353             : !>               variables at the Generalized Cauchy Point (GCP).
     354             : !> \param iwhere iwhere is an integer working array of dimension n used to record
     355             : !>               the status of the vector x for GCP computation.
     356             : !>               iwhere(i)=0 or -3 if x(i) is free and has bounds,
     357             : !>                         1       if x(i) is fixed at l(i), and l(i) .ne. u(i)
     358             : !>                         2       if x(i) is fixed at u(i), and u(i) .ne. l(i)
     359             : !>                         3       if x(i) is always fixed, i.e.,  u(i)=x(i)=l(i)
     360             : !>                        -1       if x(i) is always free, i.e., no bounds on it.
     361             : !> \param indx2  indx2 is a working array. Within subroutine cauchy, indx2 corresponds to the array iorder.
     362             : !>               In subroutine freev, a list of variables entering and leaving
     363             : !>               the free set is stored in indx2, and it is passed on to
     364             : !>               subroutine formk with this information
     365             : !> \param task   task is a working string of characters indicating
     366             : !>               the current job when entering and leaving this subroutine.
     367             : !> \param iprint is an variable that must be set by the user.
     368             : !>               It controls the frequency and type of output generated:
     369             : !>               iprint<0    no output is generated;
     370             : !>               iprint=0    print only one line at the last iteration;
     371             : !>               0<iprint<99 print also f and |proj g| every iprint iterations;
     372             : !>               iprint=99   print details of every iteration except n-vectors;
     373             : !>               iprint=100  print also the changes of active set and final x;
     374             : !>               iprint>100  print details of every iteration including x and g;
     375             : !>               When iprint > 0, the file iterate.dat will be created to summarize the iteration.
     376             : !> \param csave  csave is a working string of characters
     377             : !> \param lsave  lsave is a logical working array
     378             : !> \param isave  isave is an integer working array
     379             : !> \param dsave  is a double precision working array
     380             : !> \param spgr ...
     381             : !> \par History
     382             : !>      12.2020 Implementation of Space Group Symmetry [pcazade]
     383             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
     384             : !>               Optimization Technology Center.
     385             : !>               Argonne National Laboratory and Northwestern University.
     386             : !>               Written by
     387             : !>                           Ciyou Zhu
     388             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
     389             : ! **************************************************************************************************
     390        4617 :    SUBROUTINE mainlb(n, m, x, lower_bound, upper_bound, nbd, f, g, factr, pgtol, ws, wy, &
     391        4617 :                      sy, ss, wt, wn, snd, z, r, d, t, xp, wa, &
     392        4617 :                      index, iwhere, indx2, task, &
     393             :                      iprint, csave, lsave, isave, dsave, spgr)
     394             :       INTEGER, INTENT(in)                                :: n, m
     395             :       REAL(KIND=dp), INTENT(inout)                       :: x(n)
     396             :       REAL(KIND=dp), INTENT(in)                          :: lower_bound(n), upper_bound(n)
     397             :       INTEGER                                            :: nbd(n)
     398             :       REAL(KIND=dp) :: f, g(n), factr, pgtol, ws(n, m), wy(n, m), sy(m, m), ss(m, m), wt(m, m), &
     399             :          wn(2*m, 2*m), snd(2*m, 2*m), z(n), r(n), d(n), t(n), xp(n), wa(8*m)
     400             :       INTEGER                                            :: INDEX(n), iwhere(n), indx2(n)
     401             :       CHARACTER(LEN=60)                                  :: task
     402             :       INTEGER                                            :: iprint
     403             :       CHARACTER(LEN=60)                                  :: csave
     404             :       LOGICAL                                            :: lsave(4)
     405             :       INTEGER                                            :: isave(23)
     406             :       REAL(KIND=dp)                                      :: dsave(29)
     407             :       TYPE(spgr_type), OPTIONAL, POINTER                 :: spgr
     408             : 
     409             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     410             : 
     411             :       CHARACTER(LEN=3)                                   :: word
     412             :       INTEGER                                            :: col, head, i, iback, ifun, ileave, info, &
     413             :                                                             itail, iter, itfile, iupdat, iword, k, &
     414             :                                                             nact, nenter, nfgv, nfree, nintol, &
     415             :                                                             nseg, nskip
     416             :       LOGICAL                                            :: boxed, constrained, first, &
     417             :                                                             keep_space_group, updatd, wrk, &
     418             :                                                             x_projected
     419             :       REAL(KIND=dp) :: cachyt, cpu1, cpu2, ddot, ddum, dnorm, dr, dtd, epsmch, fold, g_inf_norm, &
     420             :          gd, gdold, lnscht, rr, sbtime, step_max, stp, theta, time, time1, time2, tol, xstep
     421             : 
     422             : !     References:
     423             : !
     424             : !       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
     425             : !       memory algorithm for bound constrained optimization'',
     426             : !       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
     427             : !
     428             : !       [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
     429             : !       Subroutines for Large Scale Bound Constrained Optimization''
     430             : !       Tech. Report, NAM-11, EECS Department, Northwestern University,
     431             : !       1994.
     432             : !
     433             : !       [3] R. Byrd, J. Nocedal and R. Schnabel "Representations of
     434             : !       Quasi-Newton Matrices and their use in Limited Memory Methods'',
     435             : !       Mathematical Programming 63 (1994), no. 4, pp. 129-156.
     436             : !
     437             : !       (Postscript files of these papers are available via anonymous
     438             : !        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
     439             : !
     440             : !                           *  *  *
     441             : 
     442        4617 :       keep_space_group = .FALSE.
     443        4617 :       IF (PRESENT(spgr)) THEN
     444        4425 :          IF (ASSOCIATED(spgr)) keep_space_group = spgr%keep_space_group
     445             :       END IF
     446             : 
     447        4617 :       IF (task == 'START') THEN
     448             : 
     449         192 :          epsmch = EPSILON(one)
     450             : 
     451         192 :          CALL timer(time1)
     452             : 
     453             : !        Initialize counters and scalars when task='START'.
     454             : 
     455             : !           for the limited memory BFGS matrices:
     456         192 :          col = 0
     457         192 :          head = 1
     458         192 :          theta = one
     459         192 :          iupdat = 0
     460         192 :          updatd = .FALSE.
     461         192 :          iback = 0
     462         192 :          itail = 0
     463         192 :          iword = 0
     464         192 :          nact = 0
     465         192 :          ileave = 0
     466         192 :          nenter = 0
     467         192 :          fold = zero
     468         192 :          dnorm = zero
     469         192 :          cpu1 = zero
     470         192 :          gd = zero
     471         192 :          step_max = zero
     472         192 :          g_inf_norm = zero
     473         192 :          stp = zero
     474         192 :          gdold = zero
     475         192 :          dtd = zero
     476             : 
     477             : !           for operation counts:
     478         192 :          iter = 0
     479         192 :          nfgv = 0
     480         192 :          nseg = 0
     481         192 :          nintol = 0
     482         192 :          nskip = 0
     483         192 :          nfree = n
     484         192 :          ifun = 0
     485             : !           for stopping tolerance:
     486         192 :          tol = factr*epsmch
     487             : 
     488             : !           for measuring running time:
     489         192 :          cachyt = 0
     490         192 :          sbtime = 0
     491         192 :          lnscht = 0
     492             : 
     493             : !           'word' records the status of subspace solutions.
     494         192 :          word = '---'
     495             : 
     496             : !           'info' records the termination information.
     497         192 :          info = 0
     498             : 
     499         192 :          itfile = 8
     500         192 :          IF (iprint >= 1) THEN
     501             : !                                open a summary file 'iterate.dat'
     502           0 :             CALL open_file(file_name='iterate.dat', unit_number=itfile, file_action='WRITE', file_status='UNKNOWN')
     503             :          END IF
     504             : 
     505             : !        Check the input arguments for errors.
     506             : 
     507         192 :          CALL errclb(n, m, factr, lower_bound, upper_bound, nbd, task, info, k)
     508         192 :          IF (task(1:5) == 'ERROR') THEN
     509             :             CALL prn3lb(n, x, f, task, iprint, info, itfile, &
     510             :                         iter, nfgv, nintol, nskip, nact, g_inf_norm, &
     511             :                         zero, nseg, word, iback, stp, xstep, k, &
     512           0 :                         cachyt, sbtime, lnscht)
     513           0 :             RETURN
     514             :          END IF
     515             : 
     516         192 :          CALL prn1lb(n, m, lower_bound, upper_bound, x, iprint, itfile, epsmch)
     517             : 
     518             : !        Initialize iwhere & project x onto the feasible set.
     519             : 
     520         192 :          CALL active(n, lower_bound, upper_bound, nbd, x, iwhere, iprint, x_projected, constrained, boxed)
     521             :          ! applies rotation matrices to coordinates
     522         192 :          IF (keep_space_group) THEN
     523           0 :             CALL spgr_apply_rotations_coord(spgr, x)
     524             :          END IF
     525             : 
     526             : !        The end of the initialization.
     527         192 :          task = 'FG_START'
     528             : !        return to the driver to calculate f and g; reenter at 111.
     529             :      CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
     530             :                         iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
     531         192 :                          cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
     532         192 :          RETURN
     533             :       ELSE
     534             :          ! applies rotation matrices to coordinates
     535        4425 :          IF (keep_space_group) THEN
     536           2 :             CALL spgr_apply_rotations_coord(spgr, x)
     537           2 :             CALL spgr_apply_rotations_force(spgr, g)
     538             :          END IF
     539             : 
     540             : !          restore local variables.
     541             : 
     542        4425 :          x_projected = lsave(1)
     543        4425 :          constrained = lsave(2)
     544        4425 :          boxed = lsave(3)
     545        4425 :          updatd = lsave(4)
     546             : 
     547        4425 :          nintol = isave(1)
     548        4425 :          itfile = isave(3)
     549        4425 :          iback = isave(4)
     550        4425 :          nskip = isave(5)
     551        4425 :          head = isave(6)
     552        4425 :          col = isave(7)
     553        4425 :          itail = isave(8)
     554        4425 :          iter = isave(9)
     555        4425 :          iupdat = isave(10)
     556        4425 :          nseg = isave(12)
     557        4425 :          nfgv = isave(13)
     558        4425 :          info = isave(14)
     559        4425 :          ifun = isave(15)
     560        4425 :          iword = isave(16)
     561        4425 :          nfree = isave(17)
     562        4425 :          nact = isave(18)
     563        4425 :          ileave = isave(19)
     564        4425 :          nenter = isave(20)
     565             : 
     566        4425 :          theta = dsave(1)
     567        4425 :          fold = dsave(2)
     568        4425 :          tol = dsave(3)
     569        4425 :          dnorm = dsave(4)
     570        4425 :          epsmch = dsave(5)
     571        4425 :          cpu1 = dsave(6)
     572        4425 :          cachyt = dsave(7)
     573        4425 :          sbtime = dsave(8)
     574        4425 :          lnscht = dsave(9)
     575        4425 :          time1 = dsave(10)
     576        4425 :          gd = dsave(11)
     577        4425 :          step_max = dsave(12)
     578        4425 :          g_inf_norm = dsave(13)
     579        4425 :          stp = dsave(14)
     580        4425 :          gdold = dsave(15)
     581        4425 :          dtd = dsave(16)
     582             : 
     583             : !        After returning from the driver go to the point where execution
     584             : !        is to resume.
     585             : 
     586        4425 :          IF (task(1:4) == 'STOP') THEN
     587           0 :             IF (task(7:9) == 'CPU') THEN
     588             : !                                          restore the previous iterate.
     589           0 :                CALL dcopy(n, t, 1, x, 1)
     590           0 :                CALL dcopy(n, r, 1, g, 1)
     591           0 :                f = fold
     592             :             END IF
     593           0 :             CALL timer(time2)
     594           0 :             time = time2 - time1
     595             :             CALL prn3lb(n, x, f, task, iprint, info, itfile, &
     596             :                         iter, nfgv, nintol, nskip, nact, g_inf_norm, &
     597             :                         time, nseg, word, iback, stp, xstep, k, &
     598           0 :                         cachyt, sbtime, lnscht)
     599             :      CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
     600             :                         iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
     601           0 :                             cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
     602           0 :             RETURN
     603             :          END IF
     604             :       END IF
     605             : 
     606        4425 :       IF (.NOT. (task(1:5) == 'FG_LN' .OR. task(1:5) == 'NEW_X')) THEN
     607             : 
     608             : !     Compute f0 and g0.
     609         191 :          nfgv = 1
     610             : 
     611             : !     Compute the infinity norm of the (-) projected gradient.
     612             : 
     613         191 :          CALL projgr(n, lower_bound, upper_bound, nbd, x, g, g_inf_norm)
     614             : 
     615         191 :          IF (iprint >= 1) THEN
     616           0 :             WRITE (*, 1002) iter, f, g_inf_norm
     617           0 :             WRITE (itfile, 1003) iter, nfgv, g_inf_norm, f
     618             :          END IF
     619         191 :          IF (g_inf_norm <= pgtol) THEN
     620             : !                                terminate the algorithm.
     621           0 :             task = 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     622           0 :             CALL timer(time2)
     623           0 :             time = time2 - time1
     624             :             CALL prn3lb(n, x, f, task, iprint, info, itfile, &
     625             :                         iter, nfgv, nintol, nskip, nact, g_inf_norm, &
     626             :                         time, nseg, word, iback, stp, xstep, k, &
     627           0 :                         cachyt, sbtime, lnscht)
     628             :      CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
     629             :                         iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
     630           0 :                             cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
     631           0 :             RETURN
     632             :          END IF
     633             :       END IF
     634             : 
     635             :       first = .TRUE.
     636             :       DO WHILE (.TRUE.)
     637        4425 :       IF (.NOT. first .OR. .NOT. (task(1:5) == 'FG_LN' .OR. task(1:5) == 'NEW_X')) THEN
     638        2048 :          IF (iprint >= 99) WRITE (*, 1001) iter + 1
     639        2048 :          iword = -1
     640             : !
     641        2048 :          IF (.NOT. constrained .AND. col > 0) THEN
     642             : !                                            skip the search for GCP.
     643        1844 :             CALL dcopy(n, x, 1, z, 1)
     644        1844 :             wrk = updatd
     645        1844 :             nseg = 0
     646             :          ELSE
     647             : 
     648             : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     649             : !
     650             : !     Compute the Generalized Cauchy Point (GCP).
     651             : !
     652             : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     653             : 
     654         204 :             CALL timer(cpu1)
     655             :             CALL cauchy(n, x, lower_bound, upper_bound, nbd, g, indx2, iwhere, t, d, z, &
     656             :                         m, wy, ws, sy, wt, theta, col, head, &
     657             :                         wa(1), wa(2*m + 1), wa(4*m + 1), wa(6*m + 1), nseg, &
     658         204 :                         iprint, g_inf_norm, info, epsmch)
     659             :             ! applies rotation matrices to coordinates
     660         204 :             IF (keep_space_group) THEN
     661           1 :                CALL spgr_apply_rotations_coord(spgr, z)
     662             :             END IF
     663         204 :             IF (info /= 0) THEN
     664             : !            singular triangular system detected; refresh the lbfgs memory.
     665           0 :                IF (iprint >= 1) WRITE (*, 1005)
     666           0 :                info = 0
     667           0 :                col = 0
     668           0 :                head = 1
     669           0 :                theta = one
     670           0 :                iupdat = 0
     671           0 :                updatd = .FALSE.
     672           0 :                CALL timer(cpu2)
     673           0 :                cachyt = cachyt + cpu2 - cpu1
     674             :                first = .FALSE.
     675             :                CYCLE
     676             :             END IF
     677         204 :             CALL timer(cpu2)
     678         204 :             cachyt = cachyt + cpu2 - cpu1
     679         204 :             nintol = nintol + nseg
     680             : 
     681             : !        Count the entering and leaving variables for iter > 0;
     682             : !        find the index set of free and active variables at the GCP.
     683             : 
     684             :             CALL freev(n, nfree, index, nenter, ileave, indx2, &
     685         204 :                        iwhere, wrk, updatd, constrained, iprint, iter)
     686         204 :             nact = n - nfree
     687             : 
     688             :          END IF
     689             : 
     690             : !     If there are no free variables or B=theta*I, then
     691             : !                                        skip the subspace minimization.
     692             : 
     693        2048 :          IF (.NOT. (nfree == 0 .OR. col == 0)) THEN
     694             : 
     695             : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     696             : !
     697             : !     Subspace minimization.
     698             : !
     699             : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     700             : 
     701        1853 :             CALL timer(cpu1)
     702             : 
     703             : !     Form  the LEL^T factorization of the indefinite
     704             : !       matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
     705             : !                     [L_a -R_z           theta*S'AA'S ]
     706             : !       where     E = [-I  0]
     707             : !                     [ 0  I]
     708             : 
     709        1853 :             IF (wrk) CALL formk(n, nfree, index, nenter, ileave, indx2, iupdat, &
     710        1853 :                                 updatd, wn, snd, m, ws, wy, sy, theta, col, head, info)
     711        1853 :             IF (info /= 0) THEN
     712             : !          nonpositive definiteness in Cholesky factorization;
     713             : !          refresh the lbfgs memory and restart the iteration.
     714           0 :                IF (iprint >= 1) WRITE (*, 1006)
     715           0 :                info = 0
     716           0 :                col = 0
     717           0 :                head = 1
     718           0 :                theta = one
     719           0 :                iupdat = 0
     720           0 :                updatd = .FALSE.
     721           0 :                CALL timer(cpu2)
     722           0 :                sbtime = sbtime + cpu2 - cpu1
     723             :                first = .FALSE.
     724             :                CYCLE
     725             :             END IF
     726             : 
     727             : !        compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x)
     728             : !                                                   from 'cauchy').
     729             :             CALL cmprlb(n, m, x, g, ws, wy, sy, wt, z, r, wa, index, &
     730        1853 :                         theta, col, head, nfree, constrained, info)
     731             :             ! applies rotation matrices to coordinates
     732        1853 :             IF (keep_space_group) THEN
     733           0 :                CALL spgr_apply_rotations_force(spgr, r)
     734             :             END IF
     735        1853 :             IF (info == 0) THEN
     736             : 
     737             : !     call the direct method.
     738             : 
     739             :                CALL subsm(n, m, nfree, index, lower_bound, upper_bound, nbd, z, r, xp, ws, wy, &
     740        1853 :                           theta, x, g, col, head, iword, wa, wn, iprint, info)
     741             :                ! applies rotation matrices to coordinates
     742        1853 :                IF (keep_space_group) THEN
     743           0 :                   CALL spgr_apply_rotations_coord(spgr, z)
     744           0 :                   CALL spgr_apply_rotations_force(spgr, r)
     745             :                END IF
     746             :             END IF
     747        1853 :             IF (info /= 0) THEN
     748             : !          singular triangular system detected;
     749             : !          refresh the lbfgs memory and restart the iteration.
     750           0 :                IF (iprint >= 1) WRITE (*, 1005)
     751           0 :                info = 0
     752           0 :                col = 0
     753           0 :                head = 1
     754           0 :                theta = one
     755           0 :                iupdat = 0
     756           0 :                updatd = .FALSE.
     757           0 :                CALL timer(cpu2)
     758           0 :                sbtime = sbtime + cpu2 - cpu1
     759             :                first = .FALSE.
     760             :                CYCLE
     761             :             END IF
     762             : 
     763        1853 :             CALL timer(cpu2)
     764        1853 :             sbtime = sbtime + cpu2 - cpu1
     765             :          END IF
     766             : 
     767             : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     768             : !
     769             : !     Line search and optimality tests.
     770             : !
     771             : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     772             : 
     773             : !     Generate the search direction d:=z-x.
     774             :          ! applies rotation matrices to coordinates
     775        2048 :          IF (keep_space_group) THEN
     776           1 :             CALL spgr_apply_rotations_coord(spgr, x)
     777           1 :             CALL spgr_apply_rotations_coord(spgr, z)
     778             :          END IF
     779      926426 :          DO i = 1, n
     780      926426 :             d(i) = z(i) - x(i)
     781             :          END DO
     782        2048 :          CALL timer(cpu1)
     783             :       END IF
     784        6282 :       IF (.NOT. first .OR. .NOT. (task(1:5) == 'NEW_X')) THEN
     785             :          ! applies rotation matrices to coordinates
     786        4425 :          IF (keep_space_group) THEN
     787           2 :             CALL spgr_apply_rotations_coord(spgr, x)
     788           2 :             CALL spgr_apply_rotations_coord(spgr, z)
     789           2 :             CALL spgr_apply_rotations_force(spgr, d)
     790           2 :             CALL spgr_apply_rotations_force(spgr, g)
     791           2 :             CALL spgr_apply_rotations_force(spgr, r)
     792             :          END IF
     793             :          CALL lnsrlb(n, lower_bound, upper_bound, nbd, x, f, fold, gd, gdold, g, d, r, t, z, stp, dnorm, &
     794             :                      dtd, xstep, step_max, iter, ifun, iback, nfgv, info, task, &
     795        4425 :                      boxed, constrained, csave, isave(22), dsave(17))
     796             :          ! applies rotation matrices to coordinates
     797        4425 :          IF (keep_space_group) THEN
     798           2 :             CALL spgr_apply_rotations_coord(spgr, x)
     799           2 :             CALL spgr_apply_rotations_force(spgr, g)
     800             :          END IF
     801        4425 :          IF (info /= 0 .OR. iback >= 20) THEN
     802             : !          restore the previous iterate.
     803           0 :             CALL dcopy(n, t, 1, x, 1)
     804           0 :             CALL dcopy(n, r, 1, g, 1)
     805           0 :             f = fold
     806           0 :             IF (col == 0) THEN
     807             : !             abnormal termination.
     808           0 :                IF (info == 0) THEN
     809           0 :                   info = -9
     810             : !                restore the actual number of f and g evaluations etc.
     811           0 :                   nfgv = nfgv - 1
     812           0 :                   ifun = ifun - 1
     813           0 :                   iback = iback - 1
     814             :                END IF
     815           0 :                task = 'ABNORMAL_TERMINATION_IN_LNSRCH'
     816           0 :                iter = iter + 1
     817           0 :                CALL timer(time2)
     818           0 :                time = time2 - time1
     819             :                CALL prn3lb(n, x, f, task, iprint, info, itfile, &
     820             :                            iter, nfgv, nintol, nskip, nact, g_inf_norm, &
     821             :                            time, nseg, word, iback, stp, xstep, k, &
     822           0 :                            cachyt, sbtime, lnscht)
     823             :      CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
     824             :                         iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
     825           0 :                                cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
     826           0 :                RETURN
     827             :             ELSE
     828             : !             refresh the lbfgs memory and restart the iteration.
     829           0 :                IF (iprint >= 1) WRITE (*, 1008)
     830           0 :                IF (info == 0) nfgv = nfgv - 1
     831           0 :                info = 0
     832           0 :                col = 0
     833           0 :                head = 1
     834           0 :                theta = one
     835           0 :                iupdat = 0
     836           0 :                updatd = .FALSE.
     837           0 :                task = 'RESTART_FROM_LNSRCH'
     838           0 :                CALL timer(cpu2)
     839           0 :                lnscht = lnscht + cpu2 - cpu1
     840             :                first = .FALSE.
     841             :                CYCLE
     842             :             END IF
     843        4425 :          ELSE IF (task(1:5) == 'FG_LN') THEN
     844             : !          return to the driver for calculating f and g; reenter at 666.
     845             :      CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
     846             :                         iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
     847        2377 :                             cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
     848        2377 :             RETURN
     849             :          ELSE
     850             : !          calculate and print out the quantities related to the new X.
     851        2048 :             CALL timer(cpu2)
     852        2048 :             lnscht = lnscht + cpu2 - cpu1
     853        2048 :             iter = iter + 1
     854             : 
     855             : !        Compute the infinity norm of the projected (-)gradient.
     856             : 
     857        2048 :             CALL projgr(n, lower_bound, upper_bound, nbd, x, g, g_inf_norm)
     858             : 
     859             : !        Print iteration information.
     860             : 
     861             :             CALL prn2lb(n, x, f, g, iprint, itfile, iter, nfgv, nact, &
     862        2048 :                         g_inf_norm, nseg, word, iword, iback, stp, xstep)
     863             :      CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
     864             :                         iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
     865        2048 :                             cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
     866        3905 :             RETURN
     867             :          END IF
     868             :       END IF
     869             : 
     870             : !     Test for termination.
     871             : 
     872        1857 :       IF (g_inf_norm <= pgtol) THEN
     873             : !                                terminate the algorithm.
     874           0 :          task = 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     875           0 :          CALL timer(time2)
     876           0 :          time = time2 - time1
     877             :          CALL prn3lb(n, x, f, task, iprint, info, itfile, &
     878             :                      iter, nfgv, nintol, nskip, nact, g_inf_norm, &
     879             :                      time, nseg, word, iback, stp, xstep, k, &
     880           0 :                      cachyt, sbtime, lnscht)
     881             :      CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
     882             :                         iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
     883           0 :                          cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
     884           0 :          RETURN
     885             :       END IF
     886             : 
     887        1857 :       ddum = MAX(ABS(fold), ABS(f), one)
     888        1857 :       IF ((fold - f) <= tol*ddum) THEN
     889             : !                                        terminate the algorithm.
     890           0 :          task = 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     891           0 :          IF (iback >= 10) info = -5
     892             : !           i.e., to issue a warning if iback>10 in the line search.
     893           0 :          CALL timer(time2)
     894           0 :          time = time2 - time1
     895             :          CALL prn3lb(n, x, f, task, iprint, info, itfile, &
     896             :                      iter, nfgv, nintol, nskip, nact, g_inf_norm, &
     897             :                      time, nseg, word, iback, stp, xstep, k, &
     898           0 :                      cachyt, sbtime, lnscht)
     899             :      CALL save_local(lsave, isave, dsave, x_projected, constrained, boxed, updatd, nintol, itfile, iback, nskip, head, col, itail, &
     900             :                         iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, &
     901           0 :                          cpu1, cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
     902           0 :          RETURN
     903             :       END IF
     904             : 
     905             : !     Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's.
     906        1857 :       IF (keep_space_group) THEN
     907           0 :          CALL spgr_apply_rotations_force(spgr, g)
     908           0 :          CALL spgr_apply_rotations_force(spgr, r)
     909             :       END IF
     910      880986 :       DO i = 1, n
     911      880986 :          r(i) = g(i) - r(i)
     912             :       END DO
     913        1857 :       rr = ddot(n, r, 1, r, 1)
     914        1857 :       IF (stp == one) THEN
     915        1607 :          dr = gd - gdold
     916        1607 :          ddum = -gdold
     917             :       ELSE
     918         250 :          dr = (gd - gdold)*stp
     919         250 :          CALL dscal(n, stp, d, 1)
     920         250 :          ddum = -gdold*stp
     921             :       END IF
     922             : 
     923        1857 :       IF (dr <= epsmch*ddum) THEN
     924             : !                            skip the L-BFGS update.
     925           4 :          nskip = nskip + 1
     926           4 :          updatd = .FALSE.
     927           4 :          IF (iprint >= 1) WRITE (*, 1004) dr, ddum
     928             :          first = .FALSE.
     929             :          CYCLE
     930             :       END IF
     931             : 
     932             : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     933             : !
     934             : !     Update the L-BFGS matrix.
     935             : !
     936             : !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     937             : 
     938        1853 :       updatd = .TRUE.
     939        1853 :       iupdat = iupdat + 1
     940             : 
     941             : !     Update matrices WS and WY and form the middle matrix in B.
     942             : 
     943             :       CALL matupd(n, m, ws, wy, sy, ss, d, r, itail, &
     944        1853 :                   iupdat, col, head, theta, rr, dr, stp, dtd)
     945             : 
     946             : !     Form the upper half of the pds T = theta*SS + L*D^(-1)*L';
     947             : !        Store T in the upper triangular of the array wt;
     948             : !        Cholesky factorize T to J*J' with
     949             : !           J' stored in the upper triangular of wt.
     950             : 
     951        1853 :       CALL formt(m, wt, sy, ss, col, theta, info)
     952             : 
     953        1853 :       IF (info /= 0) THEN
     954             : !          nonpositive definiteness in Cholesky factorization;
     955             : !          refresh the lbfgs memory and restart the iteration.
     956           0 :          IF (iprint >= 1) WRITE (*, 1007)
     957           0 :          info = 0
     958           0 :          col = 0
     959           0 :          head = 1
     960           0 :          theta = one
     961           0 :          iupdat = 0
     962           0 :          updatd = .FALSE.
     963             :       END IF
     964             : 
     965             : !     Now the inverse of the middle matrix in B is
     966             : 
     967             : !       [  D^(1/2)      O ] [ -D^(1/2)  D^(-1/2)*L' ]
     968             : !       [ -L*D^(-1/2)   J ] [  0        J'          ]
     969             : 
     970             :       first = .FALSE.
     971             :       END DO
     972             : 
     973             : 1001  FORMAT(//, 'ITERATION ', i5)
     974             : 1002  FORMAT &
     975             :          (/, 'At iterate', i5, 4x, 'f= ', 1p, d12.5, 4x, '|proj g|= ', 1p, d12.5)
     976             : 1003  FORMAT(2(1x, i4), 5x, '-', 5x, '-', 3x, '-', 5x, '-', 5x, '-', 8x, '-', 3x, &
     977             :              1p, 2(1x, d10.3))
     978             : 1004  FORMAT('  ys=', 1p, e10.3, '  -gs=', 1p, e10.3, ' BFGS update SKIPPED')
     979             : 1005  FORMAT(/, &
     980             :               ' Singular triangular system detected;', /, &
     981             :               '   refresh the lbfgs memory and restart the iteration.')
     982             : 1006  FORMAT(/, &
     983             :               ' Nonpositive definiteness in Cholesky factorization in formk;', /, &
     984             :               '   refresh the lbfgs memory and restart the iteration.')
     985             : 1007  FORMAT(/, &
     986             :               ' Nonpositive definiteness in Cholesky factorization in formt;', /, &
     987             :               '   refresh the lbfgs memory and restart the iteration.')
     988             : 1008  FORMAT(/, &
     989             :               ' Bad direction in the line search;', /, &
     990             :               '   refresh the lbfgs memory and restart the iteration.')
     991             : 
     992             :       RETURN
     993             : 
     994             :    END SUBROUTINE mainlb
     995             : 
     996             : ! **************************************************************************************************
     997             : !> \brief This subroutine initializes iwhere and projects the initial x to the feasible set if necessary.
     998             : !> \param n ...
     999             : !> \param lower_bound  the lower bound on x.
    1000             : !> \param upper_bound  the upper bound on x.
    1001             : !> \param nbd ...
    1002             : !> \param x ...
    1003             : !> \param iwhere  iwhere(i)=-1  if x(i) has no bounds
    1004             : !>                           3   if l(i)=u(i)
    1005             : !>                           0   otherwise.
    1006             : !>                In cauchy, iwhere is given finer gradations.
    1007             : !> \param iprint ...
    1008             : !> \param x_projected ...
    1009             : !> \param constrained ...
    1010             : !> \param boxed ...
    1011             : !> \author        NEOS, November 1994. (Latest revision June 1996.)
    1012             : !>                Optimization Technology Center.
    1013             : !>                Argonne National Laboratory and Northwestern University.
    1014             : !>                Written by
    1015             : !>                            Ciyou Zhu
    1016             : !>                in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    1017             : ! **************************************************************************************************
    1018         192 :    SUBROUTINE active(n, lower_bound, upper_bound, nbd, x, iwhere, iprint, &
    1019             :                      x_projected, constrained, boxed)
    1020             : 
    1021             :       INTEGER, INTENT(in)                                :: n
    1022             :       REAL(KIND=dp), INTENT(in)                          :: lower_bound(n), upper_bound(n)
    1023             :       INTEGER                                            :: nbd(n)
    1024             :       REAL(KIND=dp)                                      :: x(n)
    1025             :       INTEGER, INTENT(out)                               :: iwhere(n)
    1026             :       INTEGER                                            :: iprint
    1027             :       LOGICAL                                            :: x_projected, constrained, boxed
    1028             : 
    1029             :       REAL(KIND=dp), PARAMETER                           :: zero = 0.0_dp
    1030             : 
    1031             :       INTEGER                                            :: i, nbdd
    1032             : 
    1033             : !     ************
    1034             : !     Initialize nbdd, x_projected, constrained and boxed.
    1035             : 
    1036         192 :       nbdd = 0
    1037         192 :       x_projected = .FALSE.
    1038         192 :       constrained = .FALSE.
    1039         192 :       boxed = .TRUE.
    1040             : 
    1041             : !     Project the initial x to the easible set if necessary.
    1042             : 
    1043       45453 :       DO i = 1, n
    1044       45453 :          IF (nbd(i) > 0) THEN
    1045        4512 :             IF (nbd(i) <= 2 .AND. x(i) <= lower_bound(i)) THEN
    1046           0 :                IF (x(i) < lower_bound(i)) THEN
    1047           0 :                   x_projected = .TRUE.
    1048           0 :                   x(i) = lower_bound(i)
    1049             :                END IF
    1050           0 :                nbdd = nbdd + 1
    1051        4512 :             ELSE IF (nbd(i) >= 2 .AND. x(i) >= upper_bound(i)) THEN
    1052           0 :                IF (x(i) > upper_bound(i)) THEN
    1053           0 :                   x_projected = .TRUE.
    1054           0 :                   x(i) = upper_bound(i)
    1055             :                END IF
    1056           0 :                nbdd = nbdd + 1
    1057             :             END IF
    1058             :          END IF
    1059             :       END DO
    1060             : 
    1061             : !     Initialize iwhere and assign values to constrained and boxed.
    1062             : 
    1063       45453 :       DO i = 1, n
    1064       45261 :          IF (nbd(i) /= 2) boxed = .FALSE.
    1065       45453 :          IF (nbd(i) == 0) THEN
    1066             : !                                this variable is always free
    1067       40749 :             iwhere(i) = -1
    1068             : 
    1069             : !           otherwise set x(i)=mid(x(i), u(i), l(i)).
    1070             :          ELSE
    1071        4512 :             constrained = .TRUE.
    1072        4512 :             IF (nbd(i) == 2 .AND. upper_bound(i) - lower_bound(i) <= zero) THEN
    1073             : !                   this variable is always fixed
    1074           0 :                iwhere(i) = 3
    1075             :             ELSE
    1076        4512 :                iwhere(i) = 0
    1077             :             END IF
    1078             :          END IF
    1079             :       END DO
    1080             : 
    1081         192 :       IF (iprint >= 0) THEN
    1082           0 :          IF (x_projected) WRITE (*, *)                                        &
    1083           0 :      &   'The initial X is infeasible.  Restart with its projection.'
    1084           0 :          IF (.NOT. constrained) &
    1085           0 :             WRITE (*, *) 'This problem is unconstrained.'
    1086             :       END IF
    1087             : 
    1088         192 :       IF (iprint > 0) WRITE (*, 1001) nbdd
    1089             : 
    1090             : 1001  FORMAT(/, 'At X0 ', i9, ' variables are exactly at the bounds')
    1091             : 
    1092         192 :       RETURN
    1093             : 
    1094             :    END SUBROUTINE active
    1095             : 
    1096             : ! **************************************************************************************************
    1097             : !> \brief       This subroutine computes the product of the 2m x 2m middle matrix
    1098             : !>              in the compact L-BFGS formula of B and a 2m vector v;
    1099             : !>              it returns the product in p.
    1100             : !> \param m     m is the maximum number of variable metric corrections
    1101             : !>              used to define the limited memory matrix.
    1102             : !> \param sy    sy specifies the matrix S'Y.
    1103             : !> \param wt    wt specifies the upper triangular matrix J' which is
    1104             : !>              the Cholesky factor of (thetaS'S+LD^(-1)L').
    1105             : !> \param col   col specifies the number of s-vectors (or y-vectors)
    1106             : !>              stored in the compact L-BFGS formula.
    1107             : !> \param v     v specifies vector v.
    1108             : !> \param p     p is the product Mv.
    1109             : !> \param info  info = 0 for normal return,
    1110             : !>                   = nonzero for abnormal return when the system to be solved by dtrsl is singular.
    1111             : !> \author      NEOS, November 1994. (Latest revision June 1996.)
    1112             : !>              Optimization Technology Center.
    1113             : !>              Argonne National Laboratory and Northwestern University.
    1114             : !>              Written by
    1115             : !>                          Ciyou Zhu
    1116             : !>              in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    1117             : ! **************************************************************************************************
    1118          22 :    SUBROUTINE bmv(m, sy, wt, col, v, p, info)
    1119             : 
    1120             :       INTEGER                                            :: m
    1121             :       REAL(KIND=dp)                                      :: sy(m, m), wt(m, m)
    1122             :       INTEGER                                            :: col
    1123             :       REAL(KIND=dp), INTENT(in)                          :: v(2*col)
    1124             :       REAL(KIND=dp), INTENT(out)                         :: p(2*col)
    1125             :       INTEGER, INTENT(out)                               :: info
    1126             : 
    1127             :       INTEGER                                            :: i, i2, k
    1128             :       REAL(KIND=dp)                                      :: sum
    1129             : 
    1130          22 :       IF (col == 0) RETURN
    1131             : 
    1132             : !     PART I: solve [  D^(1/2)      O ] [ p1 ] = [ v1 ]
    1133             : !                   [ -L*D^(-1/2)   J ] [ p2 ]   [ v2 ].
    1134             : 
    1135             : !       solve Jp2=v2+LD^(-1)v1.
    1136          22 :       p(col + 1) = v(col + 1)
    1137          62 :       DO i = 2, col
    1138          40 :          i2 = col + i
    1139          40 :          sum = 0.0_dp
    1140         124 :          DO k = 1, i - 1
    1141         124 :             sum = sum + sy(i, k)*v(k)/sy(k, k)
    1142             :          END DO
    1143          62 :          p(i2) = v(i2) + sum
    1144             :       END DO
    1145             : !     Solve the triangular system
    1146          22 :       CALL dtrsl(wt, m, col, p(col + 1), 11, info)
    1147          22 :       IF (info /= 0) RETURN
    1148             : 
    1149             : !       solve D^(1/2)p1=v1.
    1150          84 :       DO i = 1, col
    1151          84 :          p(i) = v(i)/SQRT(sy(i, i))
    1152             :       END DO
    1153             : 
    1154             : !     PART II: solve [ -D^(1/2)   D^(-1/2)*L'  ] [ p1 ] = [ p1 ]
    1155             : !                    [  0         J'           ] [ p2 ]   [ p2 ].
    1156             : 
    1157             : !       solve J^Tp2=p2.
    1158          22 :       CALL dtrsl(wt, m, col, p(col + 1), 01, info)
    1159          22 :       IF (info /= 0) RETURN
    1160             : 
    1161             : !       compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2)
    1162             : !                 =-D^(-1/2)p1+D^(-1)L'p2.
    1163          84 :       DO i = 1, col
    1164          84 :          p(i) = -p(i)/SQRT(sy(i, i))
    1165             :       END DO
    1166          84 :       DO i = 1, col
    1167          62 :          sum = 0._dp
    1168         146 :          DO k = i + 1, col
    1169         146 :             sum = sum + sy(k, i)*p(col + k)/sy(i, i)
    1170             :          END DO
    1171          84 :          p(i) = p(i) + sum
    1172             :       END DO
    1173             : 
    1174             :       RETURN
    1175             : 
    1176             :    END SUBROUTINE bmv
    1177             : 
    1178             : ! **************************************************************************************************
    1179             : !> \brief        For given x, l, u, g (with g_inf_norm > 0), and a limited memory
    1180             : !>               BFGS matrix B defined in terms of matrices WY, WS, WT, and
    1181             : !>               scalars head, col, and theta, this subroutine computes the
    1182             : !>               generalized Cauchy point (GCP), defined as the first local
    1183             : !>               minimizer of the quadratic
    1184             : !>
    1185             : !>                    Q(x + s) = g's + 1/2 s'Bs
    1186             : !>
    1187             : !>               along the projected gradient direction P(x-tg,l,u).
    1188             : !>               The routine returns the GCP in xcp.
    1189             : !> \param n      n is the dimension of the problem.
    1190             : !> \param x      x is the starting point for the GCP computation.
    1191             : !> \param lower_bound  the lower bound on x.
    1192             : !> \param upper_bound  the upper bound on x.
    1193             : !> \param nbd    nbd represents the type of bounds imposed on the
    1194             : !>               variables, and must be specified as follows:
    1195             : !>               nbd(i)=0 if x(i) is unbounded,
    1196             : !>                      1 if x(i) has only a lower bound,
    1197             : !>                      2 if x(i) has both lower and upper bounds, and
    1198             : !>                      3 if x(i) has only an upper bound.
    1199             : !> \param g      g is the gradient of f(x).  g must be a nonzero vector.
    1200             : !> \param iorder iorder will be used to store the breakpoints in the piecewise
    1201             : !>               linear path and free variables encountered. On exit,
    1202             : !>               iorder(1),...,iorder(nleft) are indices of breakpoints
    1203             : !>                                which have not been encountered;
    1204             : !>               iorder(nleft+1),...,iorder(nbreak) are indices of
    1205             : !>                                     encountered breakpoints; and
    1206             : !>               iorder(nfree),...,iorder(n) are indices of variables which
    1207             : !>               have no bound constraits along the search direction.
    1208             : !> \param iwhere On entry iwhere indicates only the permanently fixed (iwhere=3)
    1209             : !>               or free (iwhere= -1) components of x.
    1210             : !>               On exit iwhere records the status of the current x variables.
    1211             : !>               iwhere(i)=-3  if x(i) is free and has bounds, but is not moved
    1212             : !>                          0   if x(i) is free and has bounds, and is moved
    1213             : !>                          1   if x(i) is fixed at l(i), and l(i) .ne. u(i)
    1214             : !>                          2   if x(i) is fixed at u(i), and u(i) .ne. l(i)
    1215             : !>                          3   if x(i) is always fixed, i.e.,  u(i)=x(i)=l(i)
    1216             : !>                         -1  if x(i) is always free, i.e., it has no bounds.
    1217             : !> \param t      t will be used to store the break points.
    1218             : !> \param d      d is used to store the Cauchy direction P(x-tg)-x.
    1219             : !> \param xcp    is a double precision array of dimension n used to return the GCP on exit.
    1220             : !> \param m      m is the maximum number of variable metric corrections used to define the limited memory matrix.
    1221             : !> \param wy     ws, wy, sy, and wt are double precision arrays.
    1222             : !>               On entry they store information that defines the limited memory BFGS matrix:
    1223             : !>               wy(n,m) stores Y, a set of y-vectors;
    1224             : !> \param ws     ws(n,m) stores S, a set of s-vectors;
    1225             : !> \param sy     sy(m,m) stores S'Y;
    1226             : !> \param wt     wt(m,m) stores the Cholesky factorization of (theta*S'S+LD^(-1)L').
    1227             : !> \param theta  theta is the scaling factor specifying B_0 = theta I.
    1228             : !> \param col    col is the actual number of variable metric corrections stored so far.
    1229             : !> \param head   head is the location of the first s-vector (or y-vector in S (or Y)
    1230             : !> \param p      p will be used to store the vector p = W^(T)d.
    1231             : !> \param c      c will be used to store the vector c = W^(T)(xcp-x).
    1232             : !> \param wbp    wbp will be used to store the row of W corresponding to a breakpoint.
    1233             : !> \param v      v is a double precision working array.
    1234             : !> \param nseg   On exit nseg records the number of quadratic segments explored in searching for the GCP.
    1235             : !> \param iprint iprint is an INTEGER variable that must be set by the user.
    1236             : !>               It controls the frequency and type of output generated:
    1237             : !>               iprint<0    no output is generated;
    1238             : !>               iprint=0    print only one line at the last iteration;
    1239             : !>               0<iprint<99 print also f and |proj g| every iprint iterations;
    1240             : !>               iprint=99   print details of every iteration except n-vectors;
    1241             : !>               iprint=100  print also the changes of active set and final x;
    1242             : !>               iprint>100  print details of every iteration including x and g;
    1243             : !>               When iprint > 0, the file iterate.dat will be created to summarize the iteration.
    1244             : !> \param g_inf_norm g_inf_norm is the norm of the projected gradient at x.
    1245             : !> \param info   On entry info is 0.
    1246             : !>               On exit info = 0       for normal return,
    1247             : !>                            = nonzero for abnormal return when the the system
    1248             : !>                              used in routine bmv is singular.
    1249             : !> \param epsmch ...
    1250             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    1251             : !>               Optimization Technology Center.
    1252             : !>               Argonne National Laboratory and Northwestern University.
    1253             : !>               Written by
    1254             : !>                           Ciyou Zhu
    1255             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    1256             : ! **************************************************************************************************
    1257         204 :    SUBROUTINE cauchy(n, x, lower_bound, upper_bound, nbd, g, iorder, iwhere, t, d, xcp, &
    1258         204 :                      m, wy, ws, sy, wt, theta, col, head, p, c, wbp, &
    1259         204 :                      v, nseg, iprint, g_inf_norm, info, epsmch)
    1260             :       INTEGER, INTENT(in)                                :: n
    1261             :       REAL(KIND=dp), INTENT(in)                          :: x(n), lower_bound(n), upper_bound(n)
    1262             :       INTEGER, INTENT(in)                                :: nbd(n)
    1263             :       REAL(KIND=dp), INTENT(in)                          :: g(n)
    1264             :       INTEGER                                            :: iorder(n)
    1265             :       INTEGER, INTENT(inout)                             :: iwhere(n)
    1266             :       REAL(KIND=dp)                                      :: t(n), d(n), xcp(n)
    1267             :       INTEGER, INTENT(in)                                :: m
    1268             :       REAL(KIND=dp), INTENT(in)                          :: sy(m, m), wt(m, m), theta
    1269             :       INTEGER, INTENT(in)                                :: col
    1270             :       REAL(KIND=dp), INTENT(in)                          :: ws(n, col), wy(n, col)
    1271             :       INTEGER, INTENT(in)                                :: head
    1272             :       REAL(KIND=dp)                                      :: p(2*m), c(2*m), wbp(2*m), v(2*m)
    1273             :       INTEGER                                            :: nseg, iprint
    1274             :       REAL(KIND=dp), INTENT(in)                          :: g_inf_norm
    1275             :       INTEGER, INTENT(inout)                             :: info
    1276             :       REAL(KIND=dp)                                      :: epsmch
    1277             : 
    1278             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
    1279             : 
    1280             :       INTEGER                                            :: col2, i, ibkmin, ibp, iter, j, nbreak, &
    1281             :                                                             nfree, nleft, pointr
    1282             :       LOGICAL                                            :: bnded, xlower, xupper
    1283             :       REAL(KIND=dp)                                      :: bkmin, ddot, dibp, dibp2, dt, dtm, f1, &
    1284             :                                                             f2, f2_org, neggi, tj, tj0, tl, tsum, &
    1285             :                                                             tu, wmc, wmp, wmw, zibp
    1286             : 
    1287             : !     References:
    1288             : !
    1289             : !       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
    1290             : !       memory algorithm for bound constrained optimization'',
    1291             : !       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
    1292             : !
    1293             : !       [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
    1294             : !       Subroutines for Large Scale Bound Constrained Optimization''
    1295             : !       Tech. Report, NAM-11, EECS Department, Northwestern University,
    1296             : !       1994.
    1297             : !
    1298             : !       (Postscript files of these papers are available via anonymous
    1299             : !        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
    1300             : !
    1301             : !                           *  *  *
    1302             : !     Check the status of the variables, reset iwhere(i) if necessary;
    1303             : !       compute the Cauchy direction d and the breakpoints t; initialize
    1304             : !       the derivative f1 and the vector p = W'd (for theta = 1).
    1305             : 
    1306         204 :       IF (g_inf_norm <= zero) THEN
    1307           0 :          IF (iprint >= 0) WRITE (*, *) 'Subgnorm = 0.  GCP = X.'
    1308           0 :          CALL dcopy(n, x, 1, xcp, 1)
    1309           0 :          RETURN
    1310             :       END IF
    1311         204 :       bnded = .TRUE.
    1312         204 :       nfree = n + 1
    1313         204 :       nbreak = 0
    1314         204 :       ibkmin = 0
    1315         204 :       bkmin = zero
    1316         204 :       col2 = 2*col
    1317         204 :       f1 = zero
    1318         204 :       IF (iprint >= 99) WRITE (*, 3010)
    1319             : 
    1320             : !     We set p to zero and build it up as we determine d.
    1321             : 
    1322         260 :       DO i = 1, col2
    1323         260 :          p(i) = zero
    1324             :       END DO
    1325             : 
    1326             : !     In the following loop we determine for each variable its bound
    1327             : !        status and its breakpoint, and update p accordingly.
    1328             : !        Smallest breakpoint is identified.
    1329             : 
    1330       54531 :       DO i = 1, n
    1331       54327 :          neggi = -g(i)
    1332       54327 :          IF (iwhere(i) /= 3 .AND. iwhere(i) /= -1) THEN
    1333             : !             if x(i) is not a constant and has bounds,
    1334             : !             compute the difference between x(i) and its bounds.
    1335       13590 :             IF (nbd(i) <= 2) tl = x(i) - lower_bound(i)
    1336       13590 :             IF (nbd(i) >= 2) tu = upper_bound(i) - x(i)
    1337             : 
    1338             : !           If a variable is close enough to a bound
    1339             : !             we treat it as at bound.
    1340       13590 :             xlower = nbd(i) <= 2 .AND. tl <= zero
    1341       13590 :             xupper = nbd(i) >= 2 .AND. tu <= zero
    1342             : 
    1343             : !              reset iwhere(i).
    1344       13590 :             iwhere(i) = 0
    1345       13590 :             IF (xlower) THEN
    1346           0 :                IF (neggi <= zero) iwhere(i) = 1
    1347       13590 :             ELSE IF (xupper) THEN
    1348           0 :                IF (neggi >= zero) iwhere(i) = 2
    1349             :             ELSE
    1350       13590 :                IF (ABS(neggi) <= zero) iwhere(i) = -3
    1351             :             END IF
    1352             :          END IF
    1353       54327 :          pointr = head
    1354       54531 :          IF (iwhere(i) /= 0 .AND. iwhere(i) /= -1) THEN
    1355           0 :             d(i) = zero
    1356             :          ELSE
    1357       54327 :             d(i) = neggi
    1358       54327 :             f1 = f1 - neggi*neggi
    1359             : !             calculate p := p - W'e_i* (g_i).
    1360       67995 :             DO j = 1, col
    1361       13668 :                p(j) = p(j) + wy(i, pointr)*neggi
    1362       13668 :                p(col + j) = p(col + j) + ws(i, pointr)*neggi
    1363       67995 :                pointr = MOD(pointr, m) + 1
    1364             :             END DO
    1365             :             IF (nbd(i) <= 2 .AND. nbd(i) /= 0                       &
    1366       54327 :      &                        .AND. neggi < zero) THEN
    1367             : !                                 x(i) + d(i) is bounded; compute t(i).
    1368        7107 :                nbreak = nbreak + 1
    1369        7107 :                iorder(nbreak) = i
    1370        7107 :                t(nbreak) = tl/(-neggi)
    1371        7107 :                IF (nbreak == 1 .OR. t(nbreak) < bkmin) THEN
    1372          40 :                   bkmin = t(nbreak)
    1373          40 :                   ibkmin = nbreak
    1374             :                END IF
    1375       47220 :             ELSE IF (nbd(i) >= 2 .AND. neggi > zero) THEN
    1376             : !                                 x(i) + d(i) is bounded; compute t(i).
    1377        6483 :                nbreak = nbreak + 1
    1378        6483 :                iorder(nbreak) = i
    1379        6483 :                t(nbreak) = tu/neggi
    1380        6483 :                IF (nbreak == 1 .OR. t(nbreak) < bkmin) THEN
    1381          19 :                   bkmin = t(nbreak)
    1382          19 :                   ibkmin = nbreak
    1383             :                END IF
    1384             :             ELSE
    1385             : !                x(i) + d(i) is not bounded.
    1386       40737 :                nfree = nfree - 1
    1387       40737 :                iorder(nfree) = i
    1388       40737 :                IF (ABS(neggi) > zero) bnded = .FALSE.
    1389             :             END IF
    1390             :          END IF
    1391             :       END DO
    1392             : 
    1393             : !     The indices of the nonzero components of d are now stored
    1394             : !       in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n).
    1395             : !       The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin.
    1396             : 
    1397         204 :       IF (theta /= one) THEN
    1398             : !                   complete the initialization of p for theta not= one.
    1399           9 :          CALL dscal(col, theta, p(col + 1), 1)
    1400             :       END IF
    1401             : 
    1402             : !     Initialize GCP xcp = x.
    1403             : 
    1404         204 :       CALL dcopy(n, x, 1, xcp, 1)
    1405             : 
    1406         204 :       IF (nbreak == 0 .AND. nfree == n + 1) THEN
    1407             : !                  is a zero vector, return with the initial xcp as GCP.
    1408           0 :          IF (iprint > 100) WRITE (*, 1010) (xcp(i), i=1, n)
    1409           0 :          RETURN
    1410             :       END IF
    1411             : 
    1412             : !     Initialize c = W'(xcp - x) = 0.
    1413             : 
    1414         260 :       DO j = 1, col2
    1415         260 :          c(j) = zero
    1416             :       END DO
    1417             : 
    1418             : !     Initialize derivative f2.
    1419             : 
    1420         204 :       f2 = -theta*f1
    1421         204 :       f2_org = f2
    1422         204 :       IF (col > 0) THEN
    1423           9 :          CALL bmv(m, sy, wt, col, p, v, info)
    1424           9 :          IF (info /= 0) RETURN
    1425           9 :          f2 = f2 - ddot(col2, v, 1, p, 1)
    1426             :       END IF
    1427         204 :       dtm = -f1/f2
    1428         204 :       tsum = zero
    1429         204 :       nseg = 1
    1430         204 :       IF (iprint >= 99) &
    1431           0 :          WRITE (*, *) 'There are ', nbreak, '  breakpoints '
    1432             : 
    1433         204 :       nleft = nbreak
    1434         204 :       iter = 1
    1435             : 
    1436         204 :       tj = zero
    1437             : 
    1438             : !     If there are no breakpoints, locate the GCP and return.
    1439             : 
    1440         204 :       IF (nleft == 0) THEN
    1441         189 :          IF (iprint >= 99) THEN
    1442           0 :             WRITE (*, *)
    1443           0 :             WRITE (*, *) 'GCP found in this segment'
    1444           0 :             WRITE (*, 4010) nseg, f1, f2
    1445           0 :             WRITE (*, 6010) dtm
    1446             :          END IF
    1447         189 :          IF (dtm <= zero) dtm = zero
    1448         189 :          tsum = tsum + dtm
    1449             : 
    1450             : !        Move free variables (i.e., the ones w/o breakpoints) and
    1451             : !          the variables whose breakpoints haven't been reached.
    1452             : 
    1453         189 :          CALL daxpy(n, tsum, d, 1, xcp, 1)
    1454             :       END IF
    1455             : 
    1456         208 :       DO WHILE (nleft > 0)
    1457             : 
    1458             : !     Find the next smallest breakpoint;
    1459             : !       compute dt = t(nleft) - t(nleft + 1).
    1460             : 
    1461          19 :          tj0 = tj
    1462          19 :          IF (iter == 1) THEN
    1463             : !         Since we already have the smallest breakpoint we need not do
    1464             : !         heapsort yet. Often only one breakpoint is used and the
    1465             : !         cost of heapsort is avoided.
    1466          15 :             tj = bkmin
    1467          15 :             ibp = iorder(ibkmin)
    1468             :          ELSE
    1469           4 :             IF (iter == 2) THEN
    1470             : !             Replace the already used smallest breakpoint with the
    1471             : !             breakpoint numbered nbreak > nlast, before heapsort call.
    1472           2 :                IF (ibkmin /= nbreak) THEN
    1473           2 :                   t(ibkmin) = t(nbreak)
    1474           2 :                   iorder(ibkmin) = iorder(nbreak)
    1475             :                END IF
    1476             : !        Update heap structure of breakpoints
    1477             : !           (if iter=2, initialize heap).
    1478             :             END IF
    1479           4 :             CALL hpsolb(nleft, t, iorder, iter - 2)
    1480           4 :             tj = t(nleft)
    1481           4 :             ibp = iorder(nleft)
    1482             :          END IF
    1483             : 
    1484          19 :          dt = tj - tj0
    1485             : 
    1486          19 :          IF (dt /= zero .AND. iprint >= 100) THEN
    1487           0 :             WRITE (*, 4011) nseg, f1, f2
    1488           0 :             WRITE (*, 5010) dt
    1489           0 :             WRITE (*, 6010) dtm
    1490             :          END IF
    1491             : 
    1492             : !     If a minimizer is within this interval, locate the GCP and return.
    1493             : 
    1494          19 :          IF (dtm < dt) THEN
    1495          15 :             IF (iprint >= 99) THEN
    1496           0 :                WRITE (*, *)
    1497           0 :                WRITE (*, *) 'GCP found in this segment'
    1498           0 :                WRITE (*, 4010) nseg, f1, f2
    1499           0 :                WRITE (*, 6010) dtm
    1500             :             END IF
    1501          15 :             IF (dtm <= zero) dtm = zero
    1502          15 :             tsum = tsum + dtm
    1503             : 
    1504             : !        Move free variables (i.e., the ones w/o breakpoints) and
    1505             : !          the variables whose breakpoints haven't been reached.
    1506             : 
    1507          15 :             CALL daxpy(n, tsum, d, 1, xcp, 1)
    1508          15 :             EXIT
    1509             :          END IF
    1510             : 
    1511             : !     Otherwise fix one variable and
    1512             : !       reset the corresponding component of d to zero.
    1513             : 
    1514           4 :          tsum = tsum + dt
    1515           4 :          nleft = nleft - 1
    1516           4 :          iter = iter + 1
    1517           4 :          dibp = d(ibp)
    1518           4 :          d(ibp) = zero
    1519           4 :          IF (dibp > zero) THEN
    1520           2 :             zibp = upper_bound(ibp) - x(ibp)
    1521           2 :             xcp(ibp) = upper_bound(ibp)
    1522           2 :             iwhere(ibp) = 2
    1523             :          ELSE
    1524           2 :             zibp = lower_bound(ibp) - x(ibp)
    1525           2 :             xcp(ibp) = lower_bound(ibp)
    1526           2 :             iwhere(ibp) = 1
    1527             :          END IF
    1528           4 :          IF (iprint >= 100) WRITE (*, *) 'Variable  ', ibp, '  is fixed.'
    1529           4 :          IF (nleft == 0 .AND. nbreak == n) THEN
    1530             : !                                             all n variables are fixed,
    1531             : !                                                return with xcp as GCP.
    1532           0 :             dtm = dt
    1533           0 :             EXIT
    1534             :          END IF
    1535             : 
    1536             : !     Update the derivative information.
    1537             : 
    1538           4 :          nseg = nseg + 1
    1539           4 :          dibp2 = dibp**2
    1540             : 
    1541             : !     Update f1 and f2.
    1542             : 
    1543             : !        temporarily set f1 and f2 for col=0.
    1544           4 :          f1 = f1 + dt*f2 + dibp2 - theta*dibp*zibp
    1545           4 :          f2 = f2 - theta*dibp2
    1546             : 
    1547           4 :          IF (col > 0) THEN
    1548             : !                          update c = c + dt*p.
    1549           4 :             CALL daxpy(col2, dt, p, 1, c, 1)
    1550             : 
    1551             : !           choose wbp,
    1552             : !           the row of W corresponding to the breakpoint encountered.
    1553           4 :             pointr = head
    1554          10 :             DO j = 1, col
    1555           6 :                wbp(j) = wy(ibp, pointr)
    1556           6 :                wbp(col + j) = theta*ws(ibp, pointr)
    1557          10 :                pointr = MOD(pointr, m) + 1
    1558             :             END DO
    1559             : 
    1560             : !           compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'.
    1561           4 :             CALL bmv(m, sy, wt, col, wbp, v, info)
    1562           4 :             IF (info /= 0) RETURN
    1563           4 :             wmc = ddot(col2, c, 1, v, 1)
    1564           4 :             wmp = ddot(col2, p, 1, v, 1)
    1565           4 :             wmw = ddot(col2, wbp, 1, v, 1)
    1566             : 
    1567             : !           update p = p - dibp*wbp.
    1568           4 :             CALL daxpy(col2, -dibp, wbp, 1, p, 1)
    1569             : 
    1570             : !           complete updating f1 and f2 while col > 0.
    1571           4 :             f1 = f1 + dibp*wmc
    1572           4 :             f2 = f2 + 2.0_dp*dibp*wmp - dibp2*wmw
    1573             :          END IF
    1574             : 
    1575           4 :          f2 = MAX(epsmch*f2_org, f2)
    1576         193 :          IF (nleft > 0) THEN
    1577           4 :             dtm = -f1/f2
    1578             :             CYCLE
    1579             : !                 to repeat the loop for unsearched intervals.
    1580             :          ELSE
    1581           0 :             IF (bnded) THEN
    1582           0 :                f1 = zero
    1583           0 :                f2 = zero
    1584           0 :                dtm = zero
    1585             :             ELSE
    1586           0 :                dtm = -f1/f2
    1587             :             END IF
    1588           0 :             IF (iprint >= 99) THEN
    1589           0 :                WRITE (*, *)
    1590           0 :                WRITE (*, *) 'GCP found in this segment'
    1591           0 :                WRITE (*, 4010) nseg, f1, f2
    1592           0 :                WRITE (*, 6010) dtm
    1593             :             END IF
    1594           0 :             IF (dtm <= zero) dtm = zero
    1595           0 :             tsum = tsum + dtm
    1596             : 
    1597             : !        Move free variables (i.e., the ones w/o breakpoints) and
    1598             : !          the variables whose breakpoints haven't been reached.
    1599             : 
    1600           0 :             CALL daxpy(n, tsum, d, 1, xcp, 1)
    1601           0 :             EXIT
    1602             :          END IF
    1603             :       END DO
    1604             : 
    1605             : !     Update c = c + dtm*p = W'(x^c - x)
    1606             : !       which will be used in computing r = Z'(B(x^c - x) + g).
    1607             : 
    1608         204 :       IF (col > 0) CALL daxpy(col2, dtm, p, 1, c, 1)
    1609         204 :       IF (iprint > 100) WRITE (*, 1010) (xcp(i), i=1, n)
    1610         204 :       IF (iprint >= 99) WRITE (*, 2010)
    1611             : 
    1612             : 1010  FORMAT('Cauchy X =  ', /, (4x, 1p, 6(1x, d11.4)))
    1613             : 2010  FORMAT(/, '---------------- exit CAUCHY----------------------',/)
    1614             : 3010  FORMAT(/, '---------------- CAUCHY entered-------------------')
    1615             : 4010  FORMAT('Piece    ', i3, ' --f1, f2 at start point ', 1p, 2(1x, d11.4))
    1616             : 4011  FORMAT(/, 'Piece    ', i3, ' --f1, f2 at start point ', &
    1617             :               1p, 2(1x, d11.4))
    1618             : 5010  FORMAT('Distance to the next break point =  ', 1p, d11.4)
    1619             : 6010  FORMAT('Distance to the stationary point =  ', 1p, d11.4)
    1620             : 
    1621             :       RETURN
    1622             : 
    1623             :    END SUBROUTINE cauchy
    1624             : 
    1625             : ! **************************************************************************************************
    1626             : !> \brief        This subroutine computes r=-Z'B(xcp-xk)-Z'g by using
    1627             : !>               wa(2m+1)=W'(xcp-x) from subroutine cauchy.
    1628             : !> \param n ...
    1629             : !> \param m ...
    1630             : !> \param x ...
    1631             : !> \param g ...
    1632             : !> \param ws ...
    1633             : !> \param wy ...
    1634             : !> \param sy ...
    1635             : !> \param wt ...
    1636             : !> \param z ...
    1637             : !> \param r ...
    1638             : !> \param wa ...
    1639             : !> \param index ...
    1640             : !> \param theta ...
    1641             : !> \param col ...
    1642             : !> \param head ...
    1643             : !> \param nfree ...
    1644             : !> \param constrained ...
    1645             : !> \param info ...
    1646             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    1647             : !>               Optimization Technology Center.
    1648             : !>               Argonne National Laboratory and Northwestern University.
    1649             : !>               Written by
    1650             : !>                           Ciyou Zhu
    1651             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    1652             : ! **************************************************************************************************
    1653        1853 :    SUBROUTINE cmprlb(n, m, x, g, ws, wy, sy, wt, z, r, wa, index, &
    1654             :                      theta, col, head, nfree, constrained, info)
    1655             : 
    1656             :       INTEGER, INTENT(in)                                :: n, m
    1657             :       REAL(KIND=dp), INTENT(in)                          :: x(n), g(n), ws(n, m), wy(n, m), &
    1658             :                                                             sy(m, m), wt(m, m), z(n)
    1659             :       REAL(KIND=dp), INTENT(out)                         :: r(n), wa(4*m)
    1660             :       INTEGER, INTENT(in)                                :: INDEX(n)
    1661             :       REAL(KIND=dp), INTENT(in)                          :: theta
    1662             :       INTEGER, INTENT(in)                                :: col, head, nfree
    1663             :       LOGICAL, INTENT(in)                                :: constrained
    1664             :       INTEGER                                            :: info
    1665             : 
    1666             :       INTEGER                                            :: i, j, k, pointr
    1667             :       REAL(KIND=dp)                                      :: a1, a2
    1668             : 
    1669        1853 :       IF (.NOT. constrained .AND. col > 0) THEN
    1670      871895 :          DO i = 1, n
    1671      871895 :             r(i) = -g(i)
    1672             :          END DO
    1673             :       ELSE
    1674        9059 :          DO i = 1, nfree
    1675        9050 :             k = INDEX(i)
    1676        9059 :             r(i) = -theta*(z(k) - x(k)) - g(k)
    1677             :          END DO
    1678           9 :          CALL bmv(m, sy, wt, col, wa(2*m + 1), wa(1), info)
    1679           9 :          IF (info /= 0) THEN
    1680           0 :             info = -8
    1681           0 :             RETURN
    1682             :          END IF
    1683           9 :          pointr = head
    1684          37 :          DO j = 1, col
    1685          28 :             a1 = wa(j)
    1686          28 :             a2 = theta*wa(col + j)
    1687       13690 :             DO i = 1, nfree
    1688       13662 :                k = INDEX(i)
    1689       13690 :                r(i) = r(i) + wy(k, pointr)*a1 + ws(k, pointr)*a2
    1690             :             END DO
    1691          37 :             pointr = MOD(pointr, m) + 1
    1692             :          END DO
    1693             :       END IF
    1694             : 
    1695             :       RETURN
    1696             : 
    1697             :    END SUBROUTINE cmprlb
    1698             : 
    1699             : ! **************************************************************************************************
    1700             : !> \brief       This subroutine checks the validity of the input data.
    1701             : !> \param n ...
    1702             : !> \param m ...
    1703             : !> \param factr ...
    1704             : !> \param lower_bound  the lower bound on x.
    1705             : !> \param upper_bound  the upper bound on x.
    1706             : !> \param nbd ...
    1707             : !> \param task ...
    1708             : !> \param info ...
    1709             : !> \param k ...
    1710             : !> \author      NEOS, November 1994. (Latest revision June 1996.)
    1711             : !>              Optimization Technology Center.
    1712             : !>              Argonne National Laboratory and Northwestern University.
    1713             : !>              Written by
    1714             : !>                          Ciyou Zhu
    1715             : !>              in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    1716             : ! **************************************************************************************************
    1717         192 :    SUBROUTINE errclb(n, m, factr, lower_bound, upper_bound, nbd, task, info, k)
    1718             : 
    1719             :       INTEGER, INTENT(in)                                :: n, m
    1720             :       REAL(KIND=dp), INTENT(in)                          :: factr, lower_bound(n), upper_bound(n)
    1721             :       INTEGER                                            :: nbd(n)
    1722             :       CHARACTER(LEN=60)                                  :: task
    1723             :       INTEGER                                            :: info, k
    1724             : 
    1725             :       REAL(KIND=dp), PARAMETER                           :: zero = 0.0_dp
    1726             : 
    1727             :       INTEGER                                            :: i
    1728             : 
    1729             : !     Check the input arguments for errors.
    1730             : 
    1731         192 :       IF (n <= 0) task = 'ERROR: N <= 0'
    1732         192 :       IF (m <= 0) task = 'ERROR: M <= 0'
    1733         192 :       IF (factr < zero) task = 'ERROR: FACTR < 0'
    1734             : 
    1735             : !     Check the validity of the arrays nbd(i), u(i), and l(i).
    1736             : 
    1737       45453 :       DO i = 1, n
    1738       45261 :          IF (nbd(i) < 0 .OR. nbd(i) > 3) THEN
    1739             : !                                                   return
    1740           0 :             task = 'ERROR: INVALID NBD'
    1741           0 :             info = -6
    1742           0 :             k = i
    1743             :          END IF
    1744       45453 :          IF (nbd(i) == 2) THEN
    1745        4512 :             IF (lower_bound(i) > upper_bound(i)) THEN
    1746             : !                                    return
    1747           0 :                task = 'ERROR: NO FEASIBLE SOLUTION'
    1748           0 :                info = -7
    1749           0 :                k = i
    1750             :             END IF
    1751             :          END IF
    1752             :       END DO
    1753             : 
    1754         192 :       RETURN
    1755             : 
    1756             :    END SUBROUTINE errclb
    1757             : 
    1758             : ! **************************************************************************************************
    1759             : !> \brief        This subroutine forms  the LEL^T factorization of the indefinite
    1760             : !>               matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
    1761             : !>                             [L_a -R_z           theta*S'AA'S ]
    1762             : !>               where     E = [-I  0]
    1763             : !>                             [ 0  I]
    1764             : !>               The matrix K can be shown to be equal to the matrix M^[-1]N
    1765             : !>               occurring in section 5.1 of [1], as well as to the matrix
    1766             : !>               Mbar^[-1] Nbar in section 5.3.
    1767             : !> \param n      n is the dimension of the problem.
    1768             : !> \param nsub   nsub is the number of subspace variables in free set.
    1769             : !> \param ind    ind specifies the indices of subspace variables.
    1770             : !> \param nenter nenter is the number of variables entering the free set.
    1771             : !> \param ileave indx2(ileave),...,indx2(n) are the variables leaving the free set.
    1772             : !> \param indx2  indx2(1),...,indx2(nenter) are the variables entering the free set,
    1773             : !>               while indx2(ileave),...,indx2(n) are the variables leaving the free set.
    1774             : !> \param iupdat iupdat is the total number of BFGS updates made so far.
    1775             : !> \param updatd 'updatd' is true if the L-BFGS matrix is updatd.
    1776             : !> \param wn     the upper triangle of wn stores the LEL^T factorization
    1777             : !>               of the 2*col x 2*col indefinite matrix
    1778             : !>                     [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
    1779             : !>                     [L_a -R_z           theta*S'AA'S ]
    1780             : !> \param wn1    On entry wn1 stores the lower triangular part of
    1781             : !>                     [Y' ZZ'Y   L_a'+R_z']
    1782             : !>                     [L_a+R_z   S'AA'S   ]
    1783             : !>               in the previous iteration.
    1784             : !>               On exit wn1 stores the corresponding updated matrices.
    1785             : !>               The purpose of wn1 is just to store these inner products
    1786             : !>               so they can be easily updated and inserted into wn.
    1787             : !> \param m      m is the maximum number of variable metric corrections
    1788             : !>               used to define the limited memory matrix.
    1789             : !> \param ws     ws(n,m) stores S, a set of s-vectors;
    1790             : !> \param wy     wy(n,m) stores Y, a set of y-vectors;
    1791             : !> \param sy     sy(m,m) stores S'Y;
    1792             : !> \param theta  is the scaling factor specifying B_0 = theta I;
    1793             : !> \param col    is the number of variable metric corrections stored;
    1794             : !> \param head   is the location of the 1st s- (or y-) vector in S (or Y).
    1795             : !> \param info   info =  0 for normal return;
    1796             : !>                    = -1 when the 1st Cholesky factorization failed;
    1797             : !>                    = -2 when the 2st Cholesky factorization failed.
    1798             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    1799             : !>               Optimization Technology Center.
    1800             : !>               Argonne National Laboratory and Northwestern University.
    1801             : !>               Written by
    1802             : !>                           Ciyou Zhu
    1803             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    1804             : ! **************************************************************************************************
    1805        1853 :    SUBROUTINE formk(n, nsub, ind, nenter, ileave, indx2, iupdat, &
    1806        1853 :                     updatd, wn, wn1, m, ws, wy, sy, theta, col, &
    1807             :                     head, info)
    1808             : 
    1809             :       INTEGER, INTENT(in)                                :: n, nsub, ind(n), nenter, ileave, &
    1810             :                                                             indx2(n), iupdat
    1811             :       LOGICAL                                            :: updatd
    1812             :       INTEGER, INTENT(in)                                :: m
    1813             :       REAL(KIND=dp)                                      :: wn1(2*m, 2*m)
    1814             :       REAL(KIND=dp), INTENT(out)                         :: wn(2*m, 2*m)
    1815             :       REAL(KIND=dp), INTENT(in)                          :: ws(n, m), wy(n, m), sy(m, m), theta
    1816             :       INTEGER, INTENT(in)                                :: col, head
    1817             :       INTEGER, INTENT(out)                               :: info
    1818             : 
    1819             :       REAL(KIND=dp), PARAMETER                           :: zero = 0.0_dp
    1820             : 
    1821             :       INTEGER                                            :: col2, dbegin, dend, i, ipntr, is, is1, &
    1822             :                                                             iy, jpntr, js, js1, jy, k, k1, m2, &
    1823             :                                                             pbegin, pend, upcl
    1824             :       REAL(KIND=dp)                                      :: ddot, temp1, temp2, temp3, temp4
    1825             : 
    1826             : !     References:
    1827             : !       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
    1828             : !       memory algorithm for bound constrained optimization'',
    1829             : !       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
    1830             : !
    1831             : !       [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
    1832             : !       limited memory FORTRAN code for solving bound constrained
    1833             : !       optimization problems'', Tech. Report, NAM-11, EECS Department,
    1834             : !       Northwestern University, 1994.
    1835             : !
    1836             : !       (Postscript files of these papers are available via anonymous
    1837             : !        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
    1838             : !
    1839             : !                           *  *  *
    1840             : !     Form the lower triangular part of
    1841             : !               WN1 = [Y' ZZ'Y   L_a'+R_z']
    1842             : !                     [L_a+R_z   S'AA'S   ]
    1843             : !        where L_a is the strictly lower triangular part of S'AA'Y
    1844             : !              R_z is the upper triangular part of S'ZZ'Y.
    1845             : 
    1846        1853 :       IF (updatd) THEN
    1847        1853 :          IF (iupdat > m) THEN
    1848             : !                                 shift old part of WN1.
    1849        6265 :             DO jy = 1, m - 1
    1850        5012 :                js = m + jy
    1851        5012 :                CALL dcopy(m - jy, wn1(jy + 1, jy + 1), 1, wn1(jy, jy), 1)
    1852        5012 :                CALL dcopy(m - jy, wn1(js + 1, js + 1), 1, wn1(js, js), 1)
    1853        6265 :                CALL dcopy(m - 1, wn1(m + 2, jy + 1), 1, wn1(m + 1, jy), 1)
    1854             :             END DO
    1855             :          END IF
    1856             : 
    1857             : !          put new rows in blocks (1,1), (2,1) and (2,2).
    1858        1853 :          pbegin = 1
    1859        1853 :          pend = nsub
    1860        1853 :          dbegin = nsub + 1
    1861        1853 :          dend = n
    1862        1853 :          iy = col
    1863        1853 :          is = m + col
    1864        1853 :          ipntr = head + col - 1
    1865        1853 :          IF (ipntr > m) ipntr = ipntr - m
    1866        1853 :          jpntr = head
    1867       10610 :          DO jy = 1, col
    1868        8757 :             js = m + jy
    1869        8757 :             temp1 = zero
    1870        8757 :             temp2 = zero
    1871        8757 :             temp3 = zero
    1872             : !             compute element jy of row 'col' of Y'ZZ'Y
    1873     4141386 :             DO k = pbegin, pend
    1874     4132629 :                k1 = ind(k)
    1875     4141386 :                temp1 = temp1 + wy(k1, ipntr)*wy(k1, jpntr)
    1876             :             END DO
    1877             : !             compute elements jy of row 'col' of L_a and S'AA'S
    1878        8763 :             DO k = dbegin, dend
    1879           6 :                k1 = ind(k)
    1880           6 :                temp2 = temp2 + ws(k1, ipntr)*ws(k1, jpntr)
    1881        8763 :                temp3 = temp3 + ws(k1, ipntr)*wy(k1, jpntr)
    1882             :             END DO
    1883        8757 :             wn1(iy, jy) = temp1
    1884        8757 :             wn1(is, js) = temp2
    1885        8757 :             wn1(is, jy) = temp3
    1886       10610 :             jpntr = MOD(jpntr, m) + 1
    1887             :          END DO
    1888             : 
    1889             : !          put new column in block (2,1).
    1890        1853 :          jy = col
    1891        1853 :          jpntr = head + col - 1
    1892        1853 :          IF (jpntr > m) jpntr = jpntr - m
    1893             :          ipntr = head
    1894       10610 :          DO i = 1, col
    1895        8757 :             is = m + i
    1896        8757 :             temp3 = zero
    1897             : !             compute element i of column 'col' of R_z
    1898     4141386 :             DO k = pbegin, pend
    1899     4132629 :                k1 = ind(k)
    1900     4141386 :                temp3 = temp3 + ws(k1, ipntr)*wy(k1, jpntr)
    1901             :             END DO
    1902        8757 :             ipntr = MOD(ipntr, m) + 1
    1903       10610 :             wn1(is, jy) = temp3
    1904             :          END DO
    1905        1853 :          upcl = col - 1
    1906             :       ELSE
    1907           0 :          upcl = col
    1908             :       END IF
    1909             : 
    1910             : !       modify the old parts in blocks (1,1) and (2,2) due to changes
    1911             : !       in the set of free variables.
    1912        1853 :       ipntr = head
    1913        8757 :       DO iy = 1, upcl
    1914        6904 :          is = m + iy
    1915        6904 :          jpntr = head
    1916       27447 :          DO jy = 1, iy
    1917       20543 :             js = m + jy
    1918       20543 :             temp1 = zero
    1919       20543 :             temp2 = zero
    1920       20543 :             temp3 = zero
    1921       20543 :             temp4 = zero
    1922       20549 :             DO k = 1, nenter
    1923           6 :                k1 = indx2(k)
    1924           6 :                temp1 = temp1 + wy(k1, ipntr)*wy(k1, jpntr)
    1925       20549 :                temp2 = temp2 + ws(k1, ipntr)*ws(k1, jpntr)
    1926             :             END DO
    1927       20543 :             DO k = ileave, n
    1928           0 :                k1 = indx2(k)
    1929           0 :                temp3 = temp3 + wy(k1, ipntr)*wy(k1, jpntr)
    1930       20543 :                temp4 = temp4 + ws(k1, ipntr)*ws(k1, jpntr)
    1931             :             END DO
    1932       20543 :             wn1(iy, jy) = wn1(iy, jy) + temp1 - temp3
    1933       20543 :             wn1(is, js) = wn1(is, js) - temp2 + temp4
    1934       27447 :             jpntr = MOD(jpntr, m) + 1
    1935             :          END DO
    1936        8757 :          ipntr = MOD(ipntr, m) + 1
    1937             :       END DO
    1938             : 
    1939             : !       modify the old parts in block (2,1).
    1940        1853 :       ipntr = head
    1941        8757 :       DO is = m + 1, m + upcl
    1942             :          jpntr = head
    1943       41086 :          DO jy = 1, upcl
    1944       34182 :             temp1 = zero
    1945       34182 :             temp3 = zero
    1946       34190 :             DO k = 1, nenter
    1947           8 :                k1 = indx2(k)
    1948       34190 :                temp1 = temp1 + ws(k1, ipntr)*wy(k1, jpntr)
    1949             :             END DO
    1950       34182 :             DO k = ileave, n
    1951           0 :                k1 = indx2(k)
    1952       34182 :                temp3 = temp3 + ws(k1, ipntr)*wy(k1, jpntr)
    1953             :             END DO
    1954       34182 :             IF (is <= jy + m) THEN
    1955       20543 :                wn1(is, jy) = wn1(is, jy) + temp1 - temp3
    1956             :             ELSE
    1957       13639 :                wn1(is, jy) = wn1(is, jy) - temp1 + temp3
    1958             :             END IF
    1959       41086 :             jpntr = MOD(jpntr, m) + 1
    1960             :          END DO
    1961        8757 :          ipntr = MOD(ipntr, m) + 1
    1962             :       END DO
    1963             : 
    1964             : !     Form the upper triangle of WN = [D+Y' ZZ'Y/theta   -L_a'+R_z' ]
    1965             : !                                     [-L_a +R_z        S'AA'S*theta]
    1966             : 
    1967        1853 :       m2 = 2*m
    1968       10610 :       DO iy = 1, col
    1969        8757 :          is = col + iy
    1970        8757 :          is1 = m + iy
    1971       38057 :          DO jy = 1, iy
    1972       29300 :             js = col + jy
    1973       29300 :             js1 = m + jy
    1974       29300 :             wn(jy, iy) = wn1(iy, jy)/theta
    1975       38057 :             wn(js, is) = wn1(is1, js1)*theta
    1976             :          END DO
    1977       29300 :          DO jy = 1, iy - 1
    1978       29300 :             wn(jy, is) = -wn1(is1, jy)
    1979             :          END DO
    1980       38057 :          DO jy = iy, col
    1981       38057 :             wn(jy, is) = wn1(is1, jy)
    1982             :          END DO
    1983       10610 :          wn(iy, iy) = wn(iy, iy) + sy(iy, iy)
    1984             :       END DO
    1985             : 
    1986             : !     Form the upper triangle of WN= [  LL'            L^-1(-L_a'+R_z')]
    1987             : !                                    [(-L_a +R_z)L'^-1   S'AA'S*theta  ]
    1988             : 
    1989             : !        first Cholesky factor (1,1) block of wn to get LL'
    1990             : !                          with L' stored in the upper triangle of wn.
    1991        1853 :       CALL dpofa(wn, m2, col, info)
    1992        1853 :       IF (info /= 0) THEN
    1993           0 :          info = -1
    1994           0 :          RETURN
    1995             :       END IF
    1996             : !        then form L^-1(-L_a'+R_z') in the (1,2) block.
    1997        1853 :       col2 = 2*col
    1998       10610 :       DO js = col + 1, col2
    1999       10610 :          CALL dtrsl(wn, m2, col, wn(1, js), 11, info)
    2000             :       END DO
    2001             : 
    2002             : !     Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the
    2003             : !        upper triangle of (2,2) block of wn.
    2004             : 
    2005       10610 :       DO is = col + 1, col2
    2006       39910 :          DO js = is, col2
    2007       38057 :             wn(is, js) = wn(is, js) + ddot(col, wn(1, is), 1, wn(1, js), 1)
    2008             :          END DO
    2009             :       END DO
    2010             : 
    2011             : !     Cholesky factorization of (2,2) block of wn.
    2012             : 
    2013        1853 :       CALL dpofa(wn(col + 1, col + 1), m2, col, info)
    2014        1853 :       IF (info /= 0) THEN
    2015           0 :          info = -2
    2016           0 :          RETURN
    2017             :       END IF
    2018             : 
    2019             :       RETURN
    2020             : 
    2021             :    END SUBROUTINE formk
    2022             : 
    2023             : ! **************************************************************************************************
    2024             : !> \brief       This subroutine forms the upper half of the pos. def. and symm.
    2025             : !>              T = theta*SS + L*D^(-1)*L', stores T in the upper triangle
    2026             : !>              of the array wt, and performs the Cholesky factorization of T
    2027             : !>              to produce J*J', with J' stored in the upper triangle of wt.
    2028             : !> \param m ...
    2029             : !> \param wt ...
    2030             : !> \param sy ...
    2031             : !> \param ss ...
    2032             : !> \param col ...
    2033             : !> \param theta ...
    2034             : !> \param info ...
    2035             : !> \author      NEOS, November 1994. (Latest revision June 1996.)
    2036             : !>              Optimization Technology Center.
    2037             : !>              Argonne National Laboratory and Northwestern University.
    2038             : !>              Written by
    2039             : !>                          Ciyou Zhu
    2040             : !>              in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    2041             : ! **************************************************************************************************
    2042        1853 :    SUBROUTINE formt(m, wt, sy, ss, col, theta, info)
    2043             : 
    2044             :       INTEGER                                            :: m
    2045             :       REAL(KIND=dp)                                      :: wt(m, m), sy(m, m), ss(m, m)
    2046             :       INTEGER                                            :: col
    2047             :       REAL(KIND=dp)                                      :: theta
    2048             :       INTEGER                                            :: info
    2049             : 
    2050             :       REAL(KIND=dp), PARAMETER                           :: zero = 0.0_dp
    2051             : 
    2052             :       INTEGER                                            :: i, j, k, k1
    2053             :       REAL(KIND=dp)                                      :: ddum
    2054             : 
    2055             : !     Form the upper half of  T = theta*SS + L*D^(-1)*L',
    2056             : !        store T in the upper triangle of the array wt.
    2057             : 
    2058       10610 :       DO j = 1, col
    2059       10610 :          wt(1, j) = theta*ss(1, j)
    2060             :       END DO
    2061        8757 :       DO i = 2, col
    2062       29300 :          DO j = i, col
    2063       20543 :             k1 = MIN(i, j) - 1
    2064       20543 :             ddum = zero
    2065       80325 :             DO k = 1, k1
    2066       80325 :                ddum = ddum + sy(i, k)*sy(j, k)/sy(k, k)
    2067             :             END DO
    2068       27447 :             wt(i, j) = ddum + theta*ss(i, j)
    2069             :          END DO
    2070             :       END DO
    2071             : 
    2072             : !     Cholesky factorize T to J*J' with
    2073             : !        J' stored in the upper triangle of wt.
    2074             : 
    2075        1853 :       CALL dpofa(wt, m, col, info)
    2076        1853 :       IF (info /= 0) THEN
    2077           0 :          info = -3
    2078             :       END IF
    2079             : 
    2080        1853 :       RETURN
    2081             : 
    2082             :    END SUBROUTINE formt
    2083             : 
    2084             : ! **************************************************************************************************
    2085             : !> \brief        This subroutine counts the entering and leaving variables when
    2086             : !>               iter > 0, and finds the index set of free and active variables
    2087             : !>               at the GCP.
    2088             : !> \param n ...
    2089             : !> \param nfree ...
    2090             : !> \param index  for i=1,...,nfree, index(i) are the indices of free variables
    2091             : !>               for i=nfree+1,...,n, index(i) are the indices of bound variables
    2092             : !>               On entry after the first iteration, index gives
    2093             : !>               the free variables at the previous iteration.
    2094             : !>               On exit it gives the free variables based on the determination
    2095             : !>               in cauchy using the array iwhere.
    2096             : !> \param nenter ...
    2097             : !> \param ileave ...
    2098             : !> \param indx2  On exit with iter>0, indx2 indicates which variables
    2099             : !>               have changed status since the previous iteration.
    2100             : !>               For i= 1,...,nenter, indx2(i) have changed from bound to free.
    2101             : !>               For i= ileave+1,...,n, indx2(i) have changed from free to bound.
    2102             : !> \param iwhere ...
    2103             : !> \param wrk ...
    2104             : !> \param updatd ...
    2105             : !> \param constrained     A variable indicating whether bounds are present
    2106             : !> \param iprint ...
    2107             : !> \param iter ...
    2108             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    2109             : !>               Optimization Technology Center.
    2110             : !>               Argonne National Laboratory and Northwestern University.
    2111             : !>               Written by
    2112             : !>                           Ciyou Zhu
    2113             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    2114             : ! **************************************************************************************************
    2115         204 :    SUBROUTINE freev(n, nfree, index, nenter, ileave, indx2, &
    2116         204 :                     iwhere, wrk, updatd, constrained, iprint, iter)
    2117             : 
    2118             :       INTEGER                                            :: n, nfree
    2119             :       INTEGER, INTENT(inout)                             :: INDEX(n)
    2120             :       INTEGER                                            :: nenter, ileave
    2121             :       INTEGER, INTENT(out)                               :: indx2(n)
    2122             :       INTEGER                                            :: iwhere(n)
    2123             :       LOGICAL                                            :: wrk, updatd, constrained
    2124             :       INTEGER                                            :: iprint, iter
    2125             : 
    2126             :       INTEGER                                            :: i, iact, k
    2127             : 
    2128         204 :       nenter = 0
    2129         204 :       ileave = n + 1
    2130         204 :       IF (iter > 0 .AND. constrained) THEN
    2131             : !                           count the entering and leaving variables.
    2132        9087 :          DO i = 1, nfree
    2133        9074 :             k = INDEX(i)
    2134             : 
    2135             : !            write(*,*) ' k  = index(i) ', k
    2136             : !            write(*,*) ' index = ', i
    2137             : 
    2138        9087 :             IF (iwhere(k) > 0) THEN
    2139           2 :                ileave = ileave - 1
    2140           2 :                indx2(ileave) = k
    2141           2 :                IF (iprint >= 100) WRITE (*, *)                         &
    2142           0 :      &             'Variable ', k, ' leaves the set of free variables'
    2143             :             END IF
    2144             :          END DO
    2145          17 :          DO i = 1 + nfree, n
    2146           4 :             k = INDEX(i)
    2147          17 :             IF (iwhere(k) <= 0) THEN
    2148           2 :                nenter = nenter + 1
    2149           2 :                indx2(nenter) = k
    2150           2 :                IF (iprint >= 100) WRITE (*, *)                         &
    2151           0 :      &             'Variable ', k, ' enters the set of free variables'
    2152             :             END IF
    2153             :          END DO
    2154          13 :          IF (iprint >= 99) WRITE (*, *) &
    2155           0 :             n + 1 - ileave, ' variables leave; ', nenter, ' variables enter'
    2156             :       END IF
    2157         204 :       wrk = (ileave < n + 1) .OR. (nenter > 0) .OR. updatd
    2158             : 
    2159             : !     Find the index set of free and active variables at the GCP.
    2160             : 
    2161         204 :       nfree = 0
    2162         204 :       iact = n + 1
    2163       54531 :       DO i = 1, n
    2164       54531 :          IF (iwhere(i) <= 0) THEN
    2165       54323 :             nfree = nfree + 1
    2166       54323 :             INDEX(nfree) = i
    2167             :          ELSE
    2168           4 :             iact = iact - 1
    2169           4 :             INDEX(iact) = i
    2170             :          END IF
    2171             :       END DO
    2172         204 :       IF (iprint >= 99) WRITE (*, *) &
    2173           0 :          nfree, ' variables are free at GCP ', iter + 1
    2174             : 
    2175         204 :       RETURN
    2176             : 
    2177             :    END SUBROUTINE freev
    2178             : 
    2179             : ! **************************************************************************************************
    2180             : !> \brief        This subroutine sorts out the least element of t, and puts the
    2181             : !>               remaining elements of t in a heap.
    2182             : !> \param n      n is the dimension of the arrays t and iorder.
    2183             : !> \param t      On entry t stores the elements to be sorted,
    2184             : !>               On exit t(n) stores the least elements of t, and t(1) to t(n-1)
    2185             : !>               stores the remaining elements in the form of a heap.
    2186             : !> \param iorder On entry iorder(i) is the index of t(i).
    2187             : !>               On exit iorder(i) is still the index of t(i), but iorder may be
    2188             : !>               permuted in accordance with t.
    2189             : !> \param iheap  iheap should be set as follows:
    2190             : !>               iheap .eq. 0 if t(1) to t(n) is not in the form of a heap,
    2191             : !>               iheap .ne. 0 if otherwise.
    2192             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    2193             : !>               Optimization Technology Center.
    2194             : !>               Argonne National Laboratory and Northwestern University.
    2195             : !>               Written by
    2196             : !>                           Ciyou Zhu
    2197             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    2198             : ! **************************************************************************************************
    2199           4 :    SUBROUTINE hpsolb(n, t, iorder, iheap)
    2200             :       INTEGER, INTENT(in)                                :: n
    2201             :       REAL(KIND=dp), INTENT(inout)                       :: t(n)
    2202             :       INTEGER, INTENT(inout)                             :: iorder(n)
    2203             :       INTEGER, INTENT(in)                                :: iheap
    2204             : 
    2205             :       INTEGER                                            :: i, indxin, indxou, j, k
    2206             :       REAL(KIND=dp)                                      :: ddum, out
    2207             : 
    2208             : !
    2209             : !     References:
    2210             : !       Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT.
    2211             : !
    2212             : !                           *  *  *
    2213             : 
    2214           4 :       IF (iheap == 0) THEN
    2215             : 
    2216             : !        Rearrange the elements t(1) to t(n) to form a heap.
    2217             : 
    2218          10 :          DO k = 2, n
    2219           8 :             ddum = t(k)
    2220           8 :             indxin = iorder(k)
    2221             : 
    2222             : !           Add ddum to the heap.
    2223           8 :             i = k
    2224          12 :             DO WHILE (i > 1)
    2225          11 :                j = i/2
    2226          12 :                IF (ddum < t(j)) THEN
    2227           4 :                   t(i) = t(j)
    2228           4 :                   iorder(i) = iorder(j)
    2229           4 :                   i = j
    2230             :                ELSE
    2231             :                   EXIT
    2232             :                END IF
    2233             :             END DO
    2234           8 :             t(i) = ddum
    2235          10 :             iorder(i) = indxin
    2236             :          END DO
    2237             :       END IF
    2238             : 
    2239             : !     Assign to 'out' the value of t(1), the least member of the heap,
    2240             : !        and rearrange the remaining members to form a heap as
    2241             : !        elements 1 to n-1 of t.
    2242             : 
    2243           4 :       IF (n > 1) THEN
    2244           4 :          i = 1
    2245           4 :          out = t(1)
    2246           4 :          indxou = iorder(1)
    2247           4 :          ddum = t(n)
    2248           4 :          indxin = iorder(n)
    2249             : 
    2250             : !        Restore the heap
    2251           4 :          j = 2*i
    2252           8 :          DO WHILE (j <= n - 1)
    2253           6 :             IF (t(j + 1) < t(j)) j = j + 1
    2254           6 :             IF (t(j) < ddum) THEN
    2255           4 :                t(i) = t(j)
    2256           4 :                iorder(i) = iorder(j)
    2257           4 :                i = j
    2258             :             ELSE
    2259             :                EXIT
    2260             :             END IF
    2261           6 :             j = 2*i
    2262             :          END DO
    2263           4 :          t(i) = ddum
    2264           4 :          iorder(i) = indxin
    2265             : 
    2266             : !     Put the least member in t(n).
    2267             : 
    2268           4 :          t(n) = out
    2269           4 :          iorder(n) = indxou
    2270             :       END IF
    2271             : 
    2272           4 :       RETURN
    2273             : 
    2274             :    END SUBROUTINE hpsolb
    2275             : 
    2276             : ! **************************************************************************************************
    2277             : !> \brief        This subroutine calls subroutine dcsrch from the Minpack2 library
    2278             : !>               to perform the line search.  Subroutine dscrch is safeguarded so
    2279             : !>               that all trial points lie within the feasible region.
    2280             : !> \param n ...
    2281             : !> \param lower_bound  the lower bound on x.
    2282             : !> \param upper_bound  the upper bound on x.
    2283             : !> \param nbd ...
    2284             : !> \param x ...
    2285             : !> \param f ...
    2286             : !> \param fold ...
    2287             : !> \param gd ...
    2288             : !> \param gdold ...
    2289             : !> \param g ...
    2290             : !> \param d ...
    2291             : !> \param r ...
    2292             : !> \param t ...
    2293             : !> \param z ...
    2294             : !> \param stp ...
    2295             : !> \param dnorm ...
    2296             : !> \param dtd ...
    2297             : !> \param xstep ...
    2298             : !> \param step_max ...
    2299             : !> \param iter ...
    2300             : !> \param ifun ...
    2301             : !> \param iback ...
    2302             : !> \param nfgv ...
    2303             : !> \param info ...
    2304             : !> \param task ...
    2305             : !> \param boxed ...
    2306             : !> \param constrained ...
    2307             : !> \param csave ...
    2308             : !> \param isave ...
    2309             : !> \param dsave ...
    2310             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    2311             : !>               Optimization Technology Center.
    2312             : !>               Argonne National Laboratory and Northwestern University.
    2313             : !>               Written by
    2314             : !>                           Ciyou Zhu
    2315             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    2316             : ! **************************************************************************************************
    2317        4425 :    SUBROUTINE lnsrlb(n, lower_bound, upper_bound, nbd, x, f, fold, gd, gdold, g, d, r, t, &
    2318        4425 :                      z, stp, dnorm, dtd, xstep, step_max, iter, ifun, &
    2319             :                      iback, nfgv, info, task, boxed, constrained, csave, &
    2320             :                      isave, dsave)
    2321             : 
    2322             :       INTEGER, INTENT(in)                                :: n
    2323             :       REAL(KIND=dp), INTENT(in)                          :: lower_bound(n), upper_bound(n)
    2324             :       INTEGER                                            :: nbd(n)
    2325             :       REAL(KIND=dp)                                      :: x(n), f, fold, gd, gdold, g(n), d(n), &
    2326             :                                                             r(n), t(n), z(n), stp, dnorm, dtd, &
    2327             :                                                             xstep, step_max
    2328             :       INTEGER                                            :: iter, ifun, iback, nfgv, info
    2329             :       CHARACTER(LEN=60)                                  :: task
    2330             :       LOGICAL                                            :: boxed, constrained
    2331             :       CHARACTER(LEN=60)                                  :: csave
    2332             :       INTEGER                                            :: isave(2)
    2333             :       REAL(KIND=dp)                                      :: dsave(13)
    2334             : 
    2335             :       REAL(KIND=dp), PARAMETER                           :: big = 1.0E10_dp, ftol = 1.0E-3_dp, &
    2336             :                                                             gtol = 0.9_dp, one = 1.0_dp, &
    2337             :                                                             xtol = 0.1_dp, zero = 0.0_dp
    2338             : 
    2339             :       INTEGER                                            :: i
    2340             :       REAL(KIND=dp)                                      :: a1, a2, ddot
    2341             : 
    2342        4425 :       IF (.NOT. (task(1:5) == 'FG_LN')) THEN
    2343             : 
    2344        2048 :          dtd = ddot(n, d, 1, d, 1)
    2345        2048 :          dnorm = SQRT(dtd)
    2346             : 
    2347             : !     Determine the maximum step length.
    2348             : 
    2349        2048 :          step_max = big
    2350        2048 :          IF (constrained) THEN
    2351          15 :             IF (iter == 0) THEN
    2352           2 :                step_max = one
    2353             :             ELSE
    2354        9091 :                DO i = 1, n
    2355        9078 :                   a1 = d(i)
    2356        9091 :                   IF (nbd(i) /= 0) THEN
    2357        9078 :                      IF (a1 < zero .AND. nbd(i) <= 2) THEN
    2358        1834 :                         a2 = lower_bound(i) - x(i)
    2359        1834 :                         IF (a2 >= zero) THEN
    2360           0 :                            step_max = zero
    2361        1834 :                         ELSE IF (a1*step_max < a2) THEN
    2362          11 :                            step_max = a2/a1
    2363             :                         END IF
    2364        7244 :                      ELSE IF (a1 > zero .AND. nbd(i) >= 2) THEN
    2365        1410 :                         a2 = upper_bound(i) - x(i)
    2366        1410 :                         IF (a2 <= zero) THEN
    2367           0 :                            step_max = zero
    2368        1410 :                         ELSE IF (a1*step_max > a2) THEN
    2369          13 :                            step_max = a2/a1
    2370             :                         END IF
    2371             :                      END IF
    2372             :                   END IF
    2373             :                END DO
    2374             :             END IF
    2375             :          END IF
    2376             : 
    2377        2048 :          IF (iter == 0 .AND. .NOT. boxed) THEN
    2378         189 :             stp = MIN(one/dnorm, step_max)
    2379             :          ELSE
    2380        1859 :             stp = one
    2381             :          END IF
    2382             : 
    2383        2048 :          CALL dcopy(n, x, 1, t, 1)
    2384        2048 :          CALL dcopy(n, g, 1, r, 1)
    2385        2048 :          fold = f
    2386        2048 :          ifun = 0
    2387        2048 :          iback = 0
    2388        2048 :          csave = 'START'
    2389             :       END IF
    2390        4425 :       gd = ddot(n, g, 1, d, 1)
    2391        4425 :       IF (ifun == 0) THEN
    2392        2048 :          gdold = gd
    2393        2048 :          IF (gd >= zero) THEN
    2394             : !                               the directional derivative >=0.
    2395             : !                               Line search is impossible.
    2396           0 :             WRITE (*, *) ' ascent direction in projection gd = ', gd
    2397           0 :             info = -4
    2398           0 :             RETURN
    2399             :          END IF
    2400             :       END IF
    2401             : 
    2402        4425 :       CALL dcsrch(f, gd, stp, ftol, gtol, xtol, zero, step_max, csave, isave, dsave)
    2403             : 
    2404        4425 :       xstep = stp*dnorm
    2405        4425 :       IF (csave(1:4) /= 'CONV' .AND. csave(1:4) /= 'WARN') THEN
    2406        2377 :          task = 'FG_LNSRCH'
    2407        2377 :          ifun = ifun + 1
    2408        2377 :          nfgv = nfgv + 1
    2409        2377 :          iback = ifun - 1
    2410        2377 :          IF (stp == one) THEN
    2411        1859 :             CALL dcopy(n, z, 1, x, 1)
    2412             :          ELSE
    2413      144746 :             DO i = 1, n
    2414      144746 :                x(i) = stp*d(i) + t(i)
    2415             :             END DO
    2416             :          END IF
    2417             :       ELSE
    2418        2048 :          task = 'NEW_X'
    2419             :       END IF
    2420             : 
    2421             :       RETURN
    2422             : 
    2423             :    END SUBROUTINE lnsrlb
    2424             : 
    2425             : ! **************************************************************************************************
    2426             : !> \brief        This subroutine updates matrices WS and WY, and forms the middle matrix in B.
    2427             : !> \param n ...
    2428             : !> \param m ...
    2429             : !> \param ws ...
    2430             : !> \param wy ...
    2431             : !> \param sy ...
    2432             : !> \param ss ...
    2433             : !> \param d ...
    2434             : !> \param r ...
    2435             : !> \param itail ...
    2436             : !> \param iupdat ...
    2437             : !> \param col ...
    2438             : !> \param head ...
    2439             : !> \param theta ...
    2440             : !> \param rr ...
    2441             : !> \param dr ...
    2442             : !> \param stp ...
    2443             : !> \param dtd ...
    2444             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    2445             : !>               Optimization Technology Center.
    2446             : !>               Argonne National Laboratory and Northwestern University.
    2447             : !>               Written by
    2448             : !>                           Ciyou Zhu
    2449             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    2450             : ! **************************************************************************************************
    2451        1853 :    SUBROUTINE matupd(n, m, ws, wy, sy, ss, d, r, itail, &
    2452             :                      iupdat, col, head, theta, rr, dr, stp, dtd)
    2453             : 
    2454             :       INTEGER                                            :: n, m
    2455             :       REAL(KIND=dp)                                      :: ws(n, m), wy(n, m), sy(m, m), ss(m, m), &
    2456             :                                                             d(n), r(n)
    2457             :       INTEGER                                            :: itail, iupdat, col, head
    2458             :       REAL(KIND=dp)                                      :: theta, rr, dr, stp, dtd
    2459             : 
    2460             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp
    2461             : 
    2462             :       INTEGER                                            :: j, pointr
    2463             :       REAL(KIND=dp)                                      :: ddot
    2464             : 
    2465             : !     ************
    2466             : !     Set pointers for matrices WS and WY.
    2467             : 
    2468        1853 :       IF (iupdat <= m) THEN
    2469         600 :          col = iupdat
    2470         600 :          itail = MOD(head + iupdat - 2, m) + 1
    2471             :       ELSE
    2472        1253 :          itail = MOD(itail, m) + 1
    2473        1253 :          head = MOD(head, m) + 1
    2474             :       END IF
    2475             : 
    2476             : !     Update matrices WS and WY.
    2477             : 
    2478        1853 :       CALL dcopy(n, d, 1, ws(1, itail), 1)
    2479        1853 :       CALL dcopy(n, r, 1, wy(1, itail), 1)
    2480             : 
    2481             : !     Set theta=yy/ys.
    2482             : 
    2483        1853 :       theta = rr/dr
    2484             : 
    2485             : !     Form the middle matrix in B.
    2486             : 
    2487             : !        update the upper triangle of SS,
    2488             : !                                         and the lower triangle of SY:
    2489        1853 :       IF (iupdat > m) THEN
    2490             : !                              move old information
    2491        6265 :          DO j = 1, col - 1
    2492        5012 :             CALL dcopy(j, ss(2, j + 1), 1, ss(1, j), 1)
    2493        6265 :             CALL dcopy(col - j, sy(j + 1, j + 1), 1, sy(j, j), 1)
    2494             :          END DO
    2495             :       END IF
    2496             : !        add new information: the last row of SY
    2497             : !                                             and the last column of SS:
    2498        1853 :       pointr = head
    2499        8757 :       DO j = 1, col - 1
    2500        6904 :          sy(col, j) = ddot(n, d, 1, wy(1, pointr), 1)
    2501        6904 :          ss(j, col) = ddot(n, ws(1, pointr), 1, d, 1)
    2502        8757 :          pointr = MOD(pointr, m) + 1
    2503             :       END DO
    2504        1853 :       IF (stp == one) THEN
    2505        1606 :          ss(col, col) = dtd
    2506             :       ELSE
    2507         247 :          ss(col, col) = stp*stp*dtd
    2508             :       END IF
    2509        1853 :       sy(col, col) = dr
    2510             : 
    2511        1853 :       RETURN
    2512             : 
    2513             :    END SUBROUTINE matupd
    2514             : 
    2515             : ! **************************************************************************************************
    2516             : !> \brief        This subroutine prints the input data, initial point, upper and
    2517             : !>               lower bounds of each variable, machine precision, as well as
    2518             : !>               the headings of the output.
    2519             : !>
    2520             : !> \param n ...
    2521             : !> \param m ...
    2522             : !> \param lower_bound  the lower bound on x.
    2523             : !> \param upper_bound  the upper bound on x.
    2524             : !> \param x ...
    2525             : !> \param iprint ...
    2526             : !> \param itfile ...
    2527             : !> \param epsmch ...
    2528             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    2529             : !>               Optimization Technology Center.
    2530             : !>               Argonne National Laboratory and Northwestern University.
    2531             : !>               Written by
    2532             : !>                           Ciyou Zhu
    2533             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    2534             : ! **************************************************************************************************
    2535         192 :    SUBROUTINE prn1lb(n, m, lower_bound, upper_bound, x, iprint, itfile, epsmch)
    2536             : 
    2537             :       INTEGER, INTENT(in)                                :: n, m
    2538             :       REAL(KIND=dp), INTENT(in)                          :: lower_bound(n), upper_bound(n), x(n)
    2539             :       INTEGER                                            :: iprint, itfile
    2540             :       REAL(KIND=dp)                                      :: epsmch
    2541             : 
    2542             :       INTEGER                                            :: i
    2543             : 
    2544         192 :       IF (iprint >= 0) THEN
    2545           0 :          WRITE (*, 7001) epsmch
    2546           0 :          WRITE (*, *) 'N = ', n, '    M = ', m
    2547           0 :          IF (iprint >= 1) THEN
    2548           0 :             WRITE (itfile, 2001) epsmch
    2549           0 :             WRITE (itfile, *) 'N = ', n, '    M = ', m
    2550           0 :             WRITE (itfile, 9001)
    2551           0 :             IF (iprint > 100) THEN
    2552           0 :                WRITE (*, 1004) 'L =', (lower_bound(i), i=1, n)
    2553           0 :                WRITE (*, 1004) 'X0 =', (x(i), i=1, n)
    2554           0 :                WRITE (*, 1004) 'U =', (upper_bound(i), i=1, n)
    2555             :             END IF
    2556             :          END IF
    2557             :       END IF
    2558             : 
    2559             : 1004  FORMAT(/, a4, 1p, 6(1x, d11.4), /, (4x, 1p, 6(1x, d11.4)))
    2560             : 2001  FORMAT('RUNNING THE L-BFGS-B CODE', /, /, &
    2561             :              'it    = iteration number', /, &
    2562             :              'nf    = number of function evaluations', /, &
    2563             :              'nseg  = number of segments explored during the Cauchy search', /, &
    2564             :              'nact  = number of active bounds at the generalized Cauchy point' &
    2565             :              , /, &
    2566             :              'sub   = manner in which the subspace minimization terminated:' &
    2567             :              , /, '        con = converged, bnd = a bound was reached', /, &
    2568             :              'itls  = number of iterations performed in the line search', /, &
    2569             :              'stepl = step length used', /, &
    2570             :              'tstep = norm of the displacement (total step)', /, &
    2571             :              'projg = norm of the projected gradient', /, &
    2572             :              'f     = function value', /, /, &
    2573             :              '           * * *', /, /, &
    2574             :              'Machine precision =', 1p, d10.3)
    2575             : 7001  FORMAT('RUNNING THE L-BFGS-B CODE', /, /, &
    2576             :              '           * * *', /, /, &
    2577             :              'Machine precision =', 1p, d10.3)
    2578             : 9001  FORMAT(/, 3x, 'it', 3x, 'nf', 2x, 'nseg', 2x, 'nact', 2x, 'sub', 2x, 'itls', &
    2579             :               2x, 'stepl', 4x, 'tstep', 5x, 'projg', 8x, 'f')
    2580             : 
    2581         192 :       RETURN
    2582             : 
    2583             :    END SUBROUTINE prn1lb
    2584             : 
    2585             : ! **************************************************************************************************
    2586             : !> \brief        This subroutine prints out new information after a successful line search.
    2587             : !> \param n ...
    2588             : !> \param x ...
    2589             : !> \param f ...
    2590             : !> \param g ...
    2591             : !> \param iprint ...
    2592             : !> \param itfile ...
    2593             : !> \param iter ...
    2594             : !> \param nfgv ...
    2595             : !> \param nact ...
    2596             : !> \param g_inf_norm ...
    2597             : !> \param nseg ...
    2598             : !> \param word ...
    2599             : !> \param iword ...
    2600             : !> \param iback ...
    2601             : !> \param stp ...
    2602             : !> \param xstep ...
    2603             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    2604             : !>               Optimization Technology Center.
    2605             : !>               Argonne National Laboratory and Northwestern University.
    2606             : !>               Written by
    2607             : !>                           Ciyou Zhu
    2608             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    2609             : ! **************************************************************************************************
    2610        2048 :    SUBROUTINE prn2lb(n, x, f, g, iprint, itfile, iter, nfgv, nact, &
    2611             :                      g_inf_norm, nseg, word, iword, iback, stp, xstep)
    2612             : 
    2613             :       INTEGER, INTENT(in)                                :: n
    2614             :       REAL(KIND=dp), INTENT(in)                          :: x(n), f, g(n)
    2615             :       INTEGER, INTENT(in)                                :: iprint, itfile, iter, nfgv, nact
    2616             :       REAL(KIND=dp), INTENT(in)                          :: g_inf_norm
    2617             :       INTEGER, INTENT(in)                                :: nseg
    2618             :       CHARACTER(LEN=3)                                   :: word
    2619             :       INTEGER                                            :: iword, iback
    2620             :       REAL(KIND=dp)                                      :: stp, xstep
    2621             : 
    2622             :       INTEGER                                            :: i, imod
    2623             : 
    2624             : !           'word' records the status of subspace solutions.
    2625             : 
    2626        2048 :       IF (iword == 0) THEN
    2627             : !                            the subspace minimization converged.
    2628        1853 :          word = 'con'
    2629         195 :       ELSE IF (iword == 1) THEN
    2630             : !                          the subspace minimization stopped at a bound.
    2631           0 :          word = 'bnd'
    2632         195 :       ELSE IF (iword == 5) THEN
    2633             : !                             the truncated Newton step has been used.
    2634           0 :          word = 'TNT'
    2635             :       ELSE
    2636         195 :          word = '---'
    2637             :       END IF
    2638        2048 :       IF (iprint >= 99) THEN
    2639           0 :          WRITE (*, *) 'LINE SEARCH', iback, ' times; norm of step = ', xstep
    2640           0 :          WRITE (*, 2001) iter, f, g_inf_norm
    2641           0 :          IF (iprint > 100) THEN
    2642           0 :             WRITE (*, 1004) 'X =', (x(i), i=1, n)
    2643           0 :             WRITE (*, 1004) 'G =', (g(i), i=1, n)
    2644             :          END IF
    2645        2048 :       ELSE IF (iprint > 0) THEN
    2646           0 :          imod = MOD(iter, iprint)
    2647           0 :          IF (imod == 0) WRITE (*, 2001) iter, f, g_inf_norm
    2648             :       END IF
    2649        2048 :       IF (iprint >= 1) WRITE (itfile, 3001) &
    2650           0 :          iter, nfgv, nseg, nact, word, iback, stp, xstep, g_inf_norm, f
    2651             : 
    2652             : 1004  FORMAT(/, a4, 1p, 6(1x, d11.4), /, (4x, 1p, 6(1x, d11.4)))
    2653             : 2001  FORMAT &
    2654             :          (/, 'At iterate', i5, 4x, 'f= ', 1p, d12.5, 4x, '|proj g|= ', 1p, d12.5)
    2655             : 3001  FORMAT(2(1x, i4), 2(1x, i5), 2x, a3, 1x, i4, 1p, 2(2x, d7.1), 1p, 2(1x, d10.3))
    2656             : 
    2657        2048 :       RETURN
    2658             : 
    2659             :    END SUBROUTINE prn2lb
    2660             : 
    2661             : ! **************************************************************************************************
    2662             : !> \brief        This subroutine prints out information when either a built-in
    2663             : !>               convergence test is satisfied or when an error message is
    2664             : !>               generated.
    2665             : !> \param n ...
    2666             : !> \param x ...
    2667             : !> \param f ...
    2668             : !> \param task ...
    2669             : !> \param iprint ...
    2670             : !> \param info ...
    2671             : !> \param itfile ...
    2672             : !> \param iter ...
    2673             : !> \param nfgv ...
    2674             : !> \param nintol ...
    2675             : !> \param nskip ...
    2676             : !> \param nact ...
    2677             : !> \param g_inf_norm ...
    2678             : !> \param time ...
    2679             : !> \param nseg ...
    2680             : !> \param word ...
    2681             : !> \param iback ...
    2682             : !> \param stp ...
    2683             : !> \param xstep ...
    2684             : !> \param k ...
    2685             : !> \param cachyt ...
    2686             : !> \param sbtime ...
    2687             : !> \param lnscht ...
    2688             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    2689             : !>               Optimization Technology Center.
    2690             : !>               Argonne National Laboratory and Northwestern University.
    2691             : !>               Written by
    2692             : !>                           Ciyou Zhu
    2693             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    2694             : ! **************************************************************************************************
    2695           0 :    SUBROUTINE prn3lb(n, x, f, task, iprint, info, itfile, &
    2696             :                      iter, nfgv, nintol, nskip, nact, g_inf_norm, &
    2697             :                      time, nseg, word, iback, stp, xstep, k, &
    2698             :                      cachyt, sbtime, lnscht)
    2699             : 
    2700             :       INTEGER, INTENT(in)                                :: n
    2701             :       REAL(KIND=dp), INTENT(in)                          :: x(n), f
    2702             :       CHARACTER(LEN=60), INTENT(in)                      :: task
    2703             :       INTEGER, INTENT(in)                                :: iprint, info, itfile, iter, nfgv, &
    2704             :                                                             nintol, nskip, nact
    2705             :       REAL(KIND=dp), INTENT(in)                          :: g_inf_norm, time
    2706             :       INTEGER, INTENT(in)                                :: nseg
    2707             :       CHARACTER(LEN=3)                                   :: word
    2708             :       INTEGER                                            :: iback
    2709             :       REAL(KIND=dp)                                      :: stp, xstep
    2710             :       INTEGER                                            :: k
    2711             :       REAL(KIND=dp)                                      :: cachyt, sbtime, lnscht
    2712             : 
    2713             :       INTEGER                                            :: i
    2714             : 
    2715           0 :       IF (iprint >= 0 .AND. .NOT. (task(1:5) == 'ERROR')) THEN
    2716           0 :          WRITE (*, 3003)
    2717           0 :          WRITE (*, 3004)
    2718           0 :          WRITE (*, 3005) n, iter, nfgv, nintol, nskip, nact, g_inf_norm, f
    2719           0 :          IF (iprint >= 100) THEN
    2720           0 :             WRITE (*, 1004) 'X =', (x(i), i=1, n)
    2721             :          END IF
    2722           0 :          IF (iprint >= 1) WRITE (*, *) ' F =', f
    2723             :       END IF
    2724           0 :       IF (iprint >= 0) THEN
    2725           0 :          WRITE (*, 3009) task
    2726           0 :          IF (info /= 0) THEN
    2727           0 :             IF (info == -1) WRITE (*, 9011)
    2728           0 :             IF (info == -2) WRITE (*, 9012)
    2729           0 :             IF (info == -3) WRITE (*, 9013)
    2730           0 :             IF (info == -4) WRITE (*, 9014)
    2731           0 :             IF (info == -5) WRITE (*, 9015)
    2732           0 :             IF (info == -6) WRITE (*, *) ' Input nbd(', k, ') is invalid.'
    2733           0 :             IF (info == -7) &
    2734           0 :                WRITE (*, *) ' l(', k, ') > u(', k, ').  No feasible solution.'
    2735           0 :             IF (info == -8) WRITE (*, 9018)
    2736           0 :             IF (info == -9) WRITE (*, 9019)
    2737             :          END IF
    2738           0 :          IF (iprint >= 1) WRITE (*, 3007) cachyt, sbtime, lnscht
    2739           0 :          WRITE (*, 3008) time
    2740           0 :          IF (iprint >= 1) THEN
    2741           0 :             IF (info == -4 .OR. info == -9) THEN
    2742             :                WRITE (itfile, 3002) &
    2743           0 :                   iter, nfgv, nseg, nact, word, iback, stp, xstep
    2744             :             END IF
    2745           0 :             WRITE (itfile, 3009) task
    2746           0 :             IF (info /= 0) THEN
    2747           0 :                IF (info == -1) WRITE (itfile, 9011)
    2748           0 :                IF (info == -2) WRITE (itfile, 9012)
    2749           0 :                IF (info == -3) WRITE (itfile, 9013)
    2750           0 :                IF (info == -4) WRITE (itfile, 9014)
    2751           0 :                IF (info == -5) WRITE (itfile, 9015)
    2752           0 :                IF (info == -8) WRITE (itfile, 9018)
    2753           0 :                IF (info == -9) WRITE (itfile, 9019)
    2754             :             END IF
    2755           0 :             WRITE (itfile, 3008) time
    2756             :          END IF
    2757             :       END IF
    2758             : 
    2759             : 1004  FORMAT(/, a4, 1p, 6(1x, d11.4), /, (4x, 1p, 6(1x, d11.4)))
    2760             : 3002  FORMAT(2(1x, i4), 2(1x, i5), 2x, a3, 1x, i4, 1p, 2(2x, d7.1), 6x, '-', 10x, '-')
    2761             : 3003  FORMAT(/, &
    2762             :               '           * * *', /, /, &
    2763             :               'Tit   = total number of iterations', /, &
    2764             :               'Tnf   = total number of function evaluations', /, &
    2765             :               'Tnint = total number of segments explored during', &
    2766             :               ' Cauchy searches', /, &
    2767             :               'Skip  = number of BFGS updates skipped', /, &
    2768             :               'Nact  = number of active bounds at final generalized', &
    2769             :               ' Cauchy point', /, &
    2770             :               'Projg = norm of the final projected gradient', /, &
    2771             :               'F     = final function value', /, /, &
    2772             :               '           * * *')
    2773             : 3004  FORMAT(/, 3x, 'N', 4x, 'Tit', 5x, 'Tnf', 2x, 'Tnint', 2x, &
    2774             :               'Skip', 2x, 'Nact', 5x, 'Projg', 8x, 'F')
    2775             : 3005  FORMAT(i5, 2(1x, i6), (1x, i6), (2x, i4), (1x, i5), 1p, 2(2x, d10.3))
    2776             : 3007  FORMAT(/, ' Cauchy                time', 1p, e10.3, ' seconds.', / &
    2777             :               ' Subspace minimization time', 1p, e10.3, ' seconds.', / &
    2778             :               ' Line search           time', 1p, e10.3, ' seconds.')
    2779             : 3008  FORMAT(/, ' Total User time', 1p, e10.3, ' seconds.',/)
    2780             : 3009  FORMAT(/, a60)
    2781             : 9011  FORMAT(/, &
    2782             :               ' Matrix in 1st Cholesky factorization in formk is not Pos. Def.')
    2783             : 9012  FORMAT(/, &
    2784             :               ' Matrix in 2st Cholesky factorization in formk is not Pos. Def.')
    2785             : 9013  FORMAT(/, &
    2786             :               ' Matrix in the Cholesky factorization in formt is not Pos. Def.')
    2787             : 9014  FORMAT(/, &
    2788             :               ' Derivative >= 0, backtracking line search impossible.', /, &
    2789             :               '   Previous x, f and g restored.', /, &
    2790             :               ' Possible causes: 1 error in function or gradient evaluation;', /, &
    2791             :               '                  2 rounding errors dominate computation.')
    2792             : 9015  FORMAT(/, &
    2793             :               ' Warning:  more than 10 function and gradient', /, &
    2794             :               '   evaluations in the last line search.  Termination', /, &
    2795             :               '   may possibly be caused by a bad search direction.')
    2796             : 9018  FORMAT(/, ' The triangular system is singular.')
    2797             : 9019  FORMAT(/, &
    2798             :               ' Line search cannot locate an adequate point after 20 function', /, &
    2799             :               '  and gradient evaluations.  Previous x, f and g restored.', /, &
    2800             :               ' Possible causes: 1 error in function or gradient evaluation;', /, &
    2801             :               '                  2 rounding error dominate computation.')
    2802             : 
    2803           0 :       RETURN
    2804             : 
    2805             :    END SUBROUTINE prn3lb
    2806             : 
    2807             : ! **************************************************************************************************
    2808             : !> \brief        This subroutine computes the infinity norm of the projected  gradient.
    2809             : !> \param n ...
    2810             : !> \param lower_bound  the lower bound on x.
    2811             : !> \param upper_bound  the upper bound on x.
    2812             : !> \param nbd ...
    2813             : !> \param x ...
    2814             : !> \param g ...
    2815             : !> \param g_inf_norm ...
    2816             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    2817             : !>               Optimization Technology Center.
    2818             : !>               Argonne National Laboratory and Northwestern University.
    2819             : !>               Written by
    2820             : !>                           Ciyou Zhu
    2821             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    2822             : ! **************************************************************************************************
    2823        2239 :    SUBROUTINE projgr(n, lower_bound, upper_bound, nbd, x, g, g_inf_norm)
    2824             : 
    2825             :       INTEGER, INTENT(in)                                :: n
    2826             :       REAL(KIND=dp), INTENT(in)                          :: lower_bound(n), upper_bound(n)
    2827             :       INTEGER, INTENT(in)                                :: nbd(n)
    2828             :       REAL(KIND=dp), INTENT(in)                          :: x(n), g(n)
    2829             :       REAL(KIND=dp)                                      :: g_inf_norm
    2830             : 
    2831             :       REAL(KIND=dp), PARAMETER                           :: zero = 0.0_dp
    2832             : 
    2833             :       INTEGER                                            :: i
    2834             :       REAL(KIND=dp)                                      :: gi
    2835             : 
    2836        2239 :       g_inf_norm = zero
    2837      971866 :       DO i = 1, n
    2838      969627 :          gi = g(i)
    2839      969627 :          IF (nbd(i) /= 0) THEN
    2840       18102 :             IF (gi < zero) THEN
    2841        8812 :                IF (nbd(i) >= 2) gi = MAX((x(i) - upper_bound(i)), gi)
    2842             :             ELSE
    2843        9290 :                IF (nbd(i) <= 2) gi = MIN((x(i) - lower_bound(i)), gi)
    2844             :             END IF
    2845             :          END IF
    2846      971866 :          g_inf_norm = MAX(g_inf_norm, ABS(gi))
    2847             :       END DO
    2848             : 
    2849        2239 :       RETURN
    2850             : 
    2851             :    END SUBROUTINE projgr
    2852             : 
    2853             : ! **************************************************************************************************
    2854             : !> \brief        This routine contains the major changes in the updated version.
    2855             : !>               The changes are described in the accompanying paper
    2856             : !>
    2857             : !>               Jose Luis Morales, Jorge Nocedal
    2858             : !>               "Remark On Algorithm 788: L-BFGS-B: Fortran Subroutines for Large
    2859             : !>               Bound Constrained Optimization". Decemmber 27, 2010.
    2860             : !>
    2861             : !>               J.L. Morales  Departamento de Matematicas,
    2862             : !>                             Instituto Tecnologico Autonomo de Mexico
    2863             : !>                             Mexico D.F.
    2864             : !>
    2865             : !>               J, Nocedal    Department of Electrical Engineering and
    2866             : !>                            Computer Science.
    2867             : !>                             Northwestern University. Evanston, IL. USA
    2868             : !>
    2869             : !>                             January 17, 2011
    2870             : !>
    2871             : !>      *****************************************************************
    2872             : !>
    2873             : !>               Given xcp, l, u, r, an index set that specifies
    2874             : !>               the active set at xcp, and an l-BFGS matrix B
    2875             : !>               (in terms of WY, WS, SY, WT, head, col, and theta),
    2876             : !>               this subroutine computes an approximate solution
    2877             : !>               of the subspace problem
    2878             : !>
    2879             : !>               (P)   min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp)
    2880             : !>
    2881             : !>               subject to l<=x<=u
    2882             : !>                       x_i=xcp_i for all i in A(xcp)
    2883             : !>
    2884             : !>               along the subspace unconstrained Newton direction
    2885             : !>
    2886             : !>               d = -(Z'BZ)^(-1) r.
    2887             : !>
    2888             : !>               The formula for the Newton direction, given the L-BFGS matrix
    2889             : !>               and the Sherman-Morrison formula, is
    2890             : !>
    2891             : !>               d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r.
    2892             : !>
    2893             : !>               where
    2894             : !>                 K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
    2895             : !>                     [L_a -R_z           theta*S'AA'S ]
    2896             : !>
    2897             : !>               Note that this procedure for computing d differs
    2898             : !>               from that described in [1]. One can show that the matrix K is
    2899             : !>               equal to the matrix M^[-1]N in that paper.
    2900             : !> \param n      n is the dimension of the problem.
    2901             : !> \param m      m is the maximum number of variable metric corrections
    2902             : !>               used to define the limited memory matrix.
    2903             : !> \param nsub   nsub is the number of free variables.
    2904             : !> \param ind    ind specifies the coordinate indices of free variables.
    2905             : !> \param lower_bound  the lower bound on x.
    2906             : !> \param upper_bound  the upper bound on x.
    2907             : !> \param nbd    nbd represents the type of bounds imposed on the
    2908             : !>               variables, and must be specified as follows:
    2909             : !>               nbd(i)=0 if x(i) is unbounded,
    2910             : !>                      1 if x(i) has only a lower bound,
    2911             : !>                      2 if x(i) has both lower and upper bounds, and
    2912             : !>                      3 if x(i) has only an upper bound.
    2913             : !> \param x      x is a double precision array of dimension n.
    2914             : !>               On entry x specifies the Cauchy point xcp.
    2915             : !>               On exit x(i) is the minimizer of Q over the subspace of free variables.
    2916             : !> \param d      On entry d is the reduced gradient of Q at xcp.
    2917             : !>               On exit d is the Newton direction of Q.
    2918             : !> \param xp     xp is a double precision array of dimension n.
    2919             : !>               used to safeguard the projected Newton direction
    2920             : !> \param ws     ws and wy are double precision arrays;
    2921             : !>               On entry they store the information defining the limited memory BFGS matrix:
    2922             : !>               ws(n,m) stores S, a set of s-vectors;
    2923             : !> \param wy     wy(n,m) stores Y, a set of y-vectors;
    2924             : !> \param theta  theta is the scaling factor specifying B_0 = theta I;
    2925             : !> \param xx     xx holds the current iterate
    2926             : !> \param gg     gg holds the gradient at the current iterate
    2927             : !> \param col    is the number of variable metric corrections stored;
    2928             : !> \param head   head is the location of the 1st s- (or y-) vector in S (or Y).
    2929             : !> \param iword  iword specifies the status of the subspace solution.
    2930             : !>               iword = 0 if the solution is in the box,
    2931             : !>                       1 if some bound is encountered.
    2932             : !> \param wv     wv is a working array
    2933             : !> \param wn     the upper triangle of wn stores the LEL^T factorization
    2934             : !>               of the indefinite matrix
    2935             : !>
    2936             : !>               K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ]
    2937             : !>                   [L_a -R_z           theta*S'AA'S ]
    2938             : !>               where E = [-I  0]
    2939             : !>                         [ 0  I]
    2940             : !> \param iprint iprint is an INTEGER variable that must be set by the user.
    2941             : !>               It controls the frequency and type of output generated:
    2942             : !>               iprint<0    no output is generated;
    2943             : !>               iprint=0    print only one line at the last iteration;
    2944             : !>               0<iprint<99 print also f and |proj g| every iprint iterations;
    2945             : !>               iprint=99   print details of every iteration except n-vectors;
    2946             : !>               iprint=100  print also the changes of active set and final x;
    2947             : !>               iprint>100  print details of every iteration including x and g;
    2948             : !>               When iprint > 0, the file iterate.dat will be created to summarize the iteration.
    2949             : !> \param info   info = 0       for normal return,
    2950             : !>                    = nonzero for abnormal return when the matrix K is ill-conditioned.
    2951             : !> \author       NEOS, November 1994. (Latest revision June 1996.)
    2952             : !>               Optimization Technology Center.
    2953             : !>               Argonne National Laboratory and Northwestern University.
    2954             : !>               Written by
    2955             : !>                           Ciyou Zhu
    2956             : !>               in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
    2957             : ! **************************************************************************************************
    2958        1853 :    SUBROUTINE subsm(n, m, nsub, ind, lower_bound, upper_bound, nbd, x, d, xp, ws, wy, &
    2959        1853 :                     theta, xx, gg, &
    2960        1853 :                     col, head, iword, wv, wn, iprint, info)
    2961             :       INTEGER, INTENT(in)                                :: n, m, nsub, ind(nsub)
    2962             :       REAL(KIND=dp), INTENT(in)                          :: lower_bound(n), upper_bound(n)
    2963             :       INTEGER, INTENT(in)                                :: nbd(n)
    2964             :       REAL(KIND=dp), INTENT(inout)                       :: x(n), d(n)
    2965             :       REAL(KIND=dp)                                      :: xp(n)
    2966             :       REAL(KIND=dp), INTENT(in)                          :: ws(n, m), wy(n, m), theta, xx(n), gg(n)
    2967             :       INTEGER, INTENT(in)                                :: col, head
    2968             :       INTEGER, INTENT(out)                               :: iword
    2969             :       REAL(KIND=dp)                                      :: wv(2*m)
    2970             :       REAL(KIND=dp), INTENT(in)                          :: wn(2*m, 2*m)
    2971             :       INTEGER                                            :: iprint
    2972             :       INTEGER, INTENT(out)                               :: info
    2973             : 
    2974             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
    2975             : 
    2976             :       INTEGER                                            :: col2, i, ibd, j, js, jy, k, m2, pointr
    2977             :       REAL(KIND=dp)                                      :: alpha, dd_p, dk, temp1, temp2, xk
    2978             : 
    2979             : !     References:
    2980             : !
    2981             : !       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
    2982             : !       memory algorithm for bound constrained optimization'',
    2983             : !       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
    2984             : !
    2985             : !
    2986             : !
    2987             : !                           *  *  *
    2988             : !
    2989             : 
    2990        1853 :       IF (nsub <= 0) RETURN
    2991        1853 :       IF (iprint >= 99) WRITE (*, 1001)
    2992             : 
    2993             : !     Compute wv = W'Zd.
    2994             : 
    2995        1853 :       pointr = head
    2996       10610 :       DO i = 1, col
    2997             :          temp1 = zero
    2998             :          temp2 = zero
    2999     4141386 :          DO j = 1, nsub
    3000     4132629 :             k = ind(j)
    3001     4132629 :             temp1 = temp1 + wy(k, pointr)*d(j)
    3002     4141386 :             temp2 = temp2 + ws(k, pointr)*d(j)
    3003             :          END DO
    3004        8757 :          wv(i) = temp1
    3005        8757 :          wv(col + i) = theta*temp2
    3006       10610 :          pointr = MOD(pointr, m) + 1
    3007             :       END DO
    3008             : 
    3009             : !     Compute wv:=K^(-1)wv.
    3010             : 
    3011        1853 :       m2 = 2*m
    3012        1853 :       col2 = 2*col
    3013        1853 :       CALL dtrsl(wn, m2, col2, wv, 11, info)
    3014        1853 :       IF (info /= 0) RETURN
    3015       10610 :       DO i = 1, col
    3016       10610 :          wv(i) = -wv(i)
    3017             :       END DO
    3018        1853 :       CALL dtrsl(wn, m2, col2, wv, 01, info)
    3019        1853 :       IF (info /= 0) RETURN
    3020             : 
    3021             : !     Compute d = (1/theta)d + (1/theta**2)Z'W wv.
    3022             : 
    3023             :       pointr = head
    3024       10610 :       DO jy = 1, col
    3025        8757 :          js = col + jy
    3026     4141386 :          DO i = 1, nsub
    3027     4132629 :             k = ind(i)
    3028             :             d(i) = d(i) + wy(k, pointr)*wv(jy)/theta                     &
    3029     4141386 :      &                  + ws(k, pointr)*wv(js)
    3030             :          END DO
    3031       10610 :          pointr = MOD(pointr, m) + 1
    3032             :       END DO
    3033             : 
    3034        1853 :       CALL dscal(nsub, one/theta, d, 1)
    3035             : !
    3036             : !-----------------------------------------------------------------
    3037             : !     Let us try the projection, d is the Newton direction
    3038             : 
    3039        1853 :       iword = 0
    3040             : 
    3041        1853 :       CALL dcopy(n, x, 1, xp, 1)
    3042             : !
    3043      880954 :       DO i = 1, nsub
    3044      879101 :          k = ind(i)
    3045      879101 :          dk = d(i)
    3046      879101 :          xk = x(k)
    3047      880954 :          IF (nbd(k) /= 0) THEN
    3048             : !
    3049             :             ! lower bounds only
    3050        9050 :             IF (nbd(k) .EQ. 1) THEN
    3051           0 :                x(k) = MAX(lower_bound(k), xk + dk)
    3052           0 :                IF (x(k) .EQ. lower_bound(k)) iword = 1
    3053             :             ELSE
    3054             : !
    3055             :                ! upper and lower bounds
    3056        9050 :                IF (nbd(k) .EQ. 2) THEN
    3057        9050 :                   xk = MAX(lower_bound(k), xk + dk)
    3058        9050 :                   x(k) = MIN(upper_bound(k), xk)
    3059        9050 :                   IF (x(k) .EQ. lower_bound(k) .OR. x(k) .EQ. upper_bound(k)) iword = 1
    3060             :                ELSE
    3061             : !
    3062             :                   ! upper bounds only
    3063           0 :                   IF (nbd(k) .EQ. 3) THEN
    3064           0 :                      x(k) = MIN(upper_bound(k), xk + dk)
    3065           0 :                      IF (x(k) .EQ. upper_bound(k)) iword = 1
    3066             :                   END IF
    3067             :                END IF
    3068             :             END IF
    3069             : !
    3070             :             ! free variables
    3071             :          ELSE
    3072      870051 :             x(k) = xk + dk
    3073             :          END IF
    3074             :       END DO
    3075             : !
    3076        1853 :       IF (.NOT. (iword .EQ. 0)) THEN
    3077             : !
    3078             : !     check sign of the directional derivative
    3079             : !
    3080             :          dd_p = zero
    3081           0 :          DO i = 1, n
    3082           0 :             dd_p = dd_p + (x(i) - xx(i))*gg(i)
    3083             :          END DO
    3084           0 :          IF (dd_p .GT. zero) THEN
    3085           0 :             CALL dcopy(n, xp, 1, x, 1)
    3086           0 :             IF (iprint > 0) WRITE (*, *) ' Positive dir derivative in projection '
    3087           0 :             IF (iprint > 0) WRITE (*, *) ' Using the backtracking step '
    3088           0 :             alpha = one
    3089           0 :             temp1 = alpha
    3090           0 :             ibd = 0
    3091           0 :             DO i = 1, nsub
    3092           0 :                k = ind(i)
    3093           0 :                dk = d(i)
    3094           0 :                IF (nbd(k) /= 0) THEN
    3095           0 :                   IF (dk < zero .AND. nbd(k) <= 2) THEN
    3096           0 :                      temp2 = lower_bound(k) - x(k)
    3097           0 :                      IF (temp2 >= zero) THEN
    3098             :                         temp1 = zero
    3099           0 :                      ELSE IF (dk*alpha < temp2) THEN
    3100           0 :                         temp1 = temp2/dk
    3101             :                      END IF
    3102           0 :                   ELSE IF (dk > zero .AND. nbd(k) >= 2) THEN
    3103           0 :                      temp2 = upper_bound(k) - x(k)
    3104           0 :                      IF (temp2 <= zero) THEN
    3105             :                         temp1 = zero
    3106           0 :                      ELSE IF (dk*alpha > temp2) THEN
    3107           0 :                         temp1 = temp2/dk
    3108             :                      END IF
    3109             :                   END IF
    3110           0 :                   IF (temp1 < alpha) THEN
    3111           0 :                      alpha = temp1
    3112           0 :                      ibd = i
    3113             :                   END IF
    3114             :                END IF
    3115             :             END DO
    3116             : 
    3117           0 :             IF (alpha < one) THEN
    3118           0 :                dk = d(ibd)
    3119           0 :                k = ind(ibd)
    3120           0 :                IF (dk > zero) THEN
    3121           0 :                   x(k) = upper_bound(k)
    3122           0 :                   d(ibd) = zero
    3123           0 :                ELSE IF (dk < zero) THEN
    3124           0 :                   x(k) = lower_bound(k)
    3125           0 :                   d(ibd) = zero
    3126             :                END IF
    3127             :             END IF
    3128           0 :             DO i = 1, nsub
    3129           0 :                k = ind(i)
    3130           0 :                x(k) = x(k) + alpha*d(i)
    3131             :             END DO
    3132             :          END IF
    3133             :       END IF
    3134             : 
    3135        1853 :       IF (iprint >= 99) WRITE (*, 1004)
    3136             : 
    3137             : 1001  FORMAT(/, '----------------SUBSM entered-----------------',/)
    3138             : 1004  FORMAT(/, '----------------exit SUBSM --------------------',/)
    3139             : 
    3140             :       RETURN
    3141             : 
    3142             :    END SUBROUTINE subsm
    3143             : 
    3144             : ! **************************************************************************************************
    3145             : !> \brief         This subroutine finds a step that satisfies a sufficient
    3146             : !>                decrease condition and a curvature condition.
    3147             : !>
    3148             : !>                Each call of the subroutine updates an interval with
    3149             : !>                endpoints stx and sty. The interval is initially chosen
    3150             : !>                so that it contains a minimizer of the modified function
    3151             : !>
    3152             : !>                      psi(stp) = f(stp) - f(0) - ftol*stp*f'(0).
    3153             : !>
    3154             : !>                If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
    3155             : !>                interval is chosen so that it contains a minimizer of f.
    3156             : !>
    3157             : !>                The algorithm is designed to find a step that satisfies
    3158             : !>                the sufficient decrease condition
    3159             : !>
    3160             : !>                f(stp) <= f(0) + ftol*stp*f'(0),
    3161             : !>
    3162             : !>                and the curvature condition
    3163             : !>
    3164             : !>                abs(f'(stp)) <= gtol*abs(f'(0)).
    3165             : !>
    3166             : !>                If ftol is less than gtol and if, for example, the function
    3167             : !>                is bounded below, then there is always a step which satisfies
    3168             : !>                both conditions.
    3169             : !>
    3170             : !>                If no step can be found that satisfies both conditions, then
    3171             : !>                the algorithm stops with a warning. In this case stp only
    3172             : !>                satisfies the sufficient decrease condition.
    3173             : !>
    3174             : !>                A typical invocation of dcsrch has the following outline:
    3175             : !>
    3176             : !>                task = 'START'
    3177             : !>                DO WHILE (.TRUE.)
    3178             : !>                   call dcsrch( ... )
    3179             : !>                   if (task .eq. 'FG') then
    3180             : !>                      Evaluate the function and the gradient at stp
    3181             : !>                   else
    3182             : !>                      exit
    3183             : !>                   end if
    3184             : !>                END DO
    3185             : !> \param f       On initial entry f is the value of the function at 0.
    3186             : !>                On subsequent entries f is the value of the
    3187             : !>                function at stp.
    3188             : !>                On exit f is the value of the function at stp.
    3189             : !> \param g       On initial entry g is the derivative of the function at 0.
    3190             : !>                On subsequent entries g is the derivative of the
    3191             : !>                function at stp.
    3192             : !>                On exit g is the derivative of the function at stp.
    3193             : !> \param stp     On entry stp is the current estimate of a satisfactory
    3194             : !>                step. On initial entry, a positive initial estimate
    3195             : !>                must be provided.
    3196             : !>                On exit stp is the current estimate of a satisfactory step
    3197             : !>                if task = 'FG'. If task = 'CONV' then stp satisfies
    3198             : !>                the sufficient decrease and curvature condition.
    3199             : !> \param ftol    ftol specifies a nonnegative tolerance for the
    3200             : !>                sufficient decrease condition.
    3201             : !> \param gtol    gtol specifies a nonnegative tolerance for the
    3202             : !>                curvature condition.
    3203             : !> \param xtol    xtol specifies a nonnegative relative tolerance
    3204             : !>                for an acceptable step. The subroutine exits with a
    3205             : !>                warning if the relative difference between sty and stx
    3206             : !>                is less than xtol.
    3207             : !> \param stpmin  stpmin is a nonnegative lower bound for the step.
    3208             : !> \param stpmax  stpmax is a nonnegative upper bound for the step.
    3209             : !> \param task    task is a character variable of length at least 60.
    3210             : !>                On initial entry task must be set to 'START'.
    3211             : !>                On exit task indicates the required action:
    3212             : !>
    3213             : !>                If task(1:2) = 'FG' then evaluate the function and
    3214             : !>                derivative at stp and call dcsrch again.
    3215             : !>
    3216             : !>                If task(1:4) = 'CONV' then the search is successful.
    3217             : !>
    3218             : !>                If task(1:4) = 'WARN' then the subroutine is not able
    3219             : !>                to satisfy the convergence conditions. The exit value of
    3220             : !>                stp contains the best point found during the search.
    3221             : !>
    3222             : !>                If task(1:5) = 'ERROR' then there is an error in the
    3223             : !>                input arguments.
    3224             : !>
    3225             : !>                On exit with convergence, a warning or an error, the
    3226             : !>                variable task contains additional information.
    3227             : !> \param isave   is work array
    3228             : !> \param dsave   is a work array
    3229             : ! **************************************************************************************************
    3230        4425 :    SUBROUTINE dcsrch(f, g, stp, ftol, gtol, xtol, stpmin, stpmax, &
    3231             :                      task, isave, dsave)
    3232             :       REAL(KIND=dp)                                      :: f, g
    3233             :       REAL(KIND=dp), INTENT(inout)                       :: stp
    3234             :       REAL(KIND=dp)                                      :: ftol, gtol, xtol, stpmin, stpmax
    3235             :       CHARACTER(LEN=*)                                   :: task
    3236             :       INTEGER                                            :: isave(2)
    3237             :       REAL(KIND=dp)                                      :: dsave(13)
    3238             : 
    3239             :       REAL(KIND=dp), PARAMETER                           :: p5 = 0.5_dp, p66 = 0.66_dp, &
    3240             :                                                             xtrapl = 1.1_dp, xtrapu = 4.0_dp, &
    3241             :                                                             zero = 0.0_dp
    3242             : 
    3243             :       INTEGER                                            :: stage
    3244             :       LOGICAL                                            :: brackt
    3245             :       REAL(KIND=dp)                                      :: finit, fm, ftest, fx, fxm, fy, fym, &
    3246             :                                                             ginit, gm, gtest, gx, gxm, gy, gym, &
    3247             :                                                             stmax, stmin, stx, sty, width, width1
    3248             : 
    3249             : !
    3250             : !     NOTE: The user must no alter work arrays between calls.
    3251             : !
    3252             : !
    3253             : !     MINPACK-1 Project. June 1983.
    3254             : !     Argonne National Laboratory.
    3255             : !     Jorge J. More' and David J. Thuente.
    3256             : !
    3257             : !     MINPACK-2 Project. October 1993.
    3258             : !     Argonne National Laboratory and University of Minnesota.
    3259             : !     Brett M. Averick, Richard G. Carter, and Jorge J. More'.
    3260             : !
    3261             : !     **********
    3262             : !     Initialization block.
    3263             : 
    3264        4425 :       IF (task(1:5) == 'START') THEN
    3265             : 
    3266             : !        Check the input arguments for errors.
    3267             : 
    3268        2048 :          IF (stp < stpmin) task = 'ERROR: STP < STPMIN'
    3269        2048 :          IF (stp > stpmax) task = 'ERROR: STP > STPMAX'
    3270        2048 :          IF (g >= zero) task = 'ERROR: INITIAL G >= ZERO'
    3271        2048 :          IF (ftol < zero) task = 'ERROR: FTOL < ZERO'
    3272        2048 :          IF (gtol < zero) task = 'ERROR: GTOL < ZERO'
    3273        2048 :          IF (xtol < zero) task = 'ERROR: XTOL < ZERO'
    3274        2048 :          IF (stpmin < zero) task = 'ERROR: STPMIN < ZERO'
    3275        2048 :          IF (stpmax < stpmin) task = 'ERROR: STPMAX < STPMIN'
    3276             : 
    3277             : !        Exit if there are errors on input.
    3278             : 
    3279        2048 :          IF (task(1:5) == 'ERROR') RETURN
    3280             : 
    3281             : !        Initialize local variables.
    3282             : 
    3283        2048 :          brackt = .FALSE.
    3284        2048 :          stage = 1
    3285        2048 :          finit = f
    3286        2048 :          ginit = g
    3287        2048 :          gtest = ftol*ginit
    3288        2048 :          width = stpmax - stpmin
    3289        2048 :          width1 = width/p5
    3290             : 
    3291             : !        The variables stx, fx, gx contain the values of the step,
    3292             : !        function, and derivative at the best step.
    3293             : !        The variables sty, fy, gy contain the value of the step,
    3294             : !        function, and derivative at sty.
    3295             : !        The variables stp, f, g contain the values of the step,
    3296             : !        function, and derivative at stp.
    3297             : 
    3298        2048 :          stx = zero
    3299        2048 :          fx = finit
    3300        2048 :          gx = ginit
    3301        2048 :          sty = zero
    3302        2048 :          fy = finit
    3303        2048 :          gy = ginit
    3304        2048 :          stmin = zero
    3305        2048 :          stmax = stp + xtrapu*stp
    3306        2048 :          task = 'FG'
    3307             : 
    3308             :       ELSE
    3309             : 
    3310             : !        Restore local variables.
    3311             : 
    3312        2377 :          IF (isave(1) == 1) THEN
    3313         315 :             brackt = .TRUE.
    3314             :          ELSE
    3315        2062 :             brackt = .FALSE.
    3316             :          END IF
    3317        2377 :          stage = isave(2)
    3318        2377 :          ginit = dsave(1)
    3319        2377 :          gtest = dsave(2)
    3320        2377 :          gx = dsave(3)
    3321        2377 :          gy = dsave(4)
    3322        2377 :          finit = dsave(5)
    3323        2377 :          fx = dsave(6)
    3324        2377 :          fy = dsave(7)
    3325        2377 :          stx = dsave(8)
    3326        2377 :          sty = dsave(9)
    3327        2377 :          stmin = dsave(10)
    3328        2377 :          stmax = dsave(11)
    3329        2377 :          width = dsave(12)
    3330        2377 :          width1 = dsave(13)
    3331             : 
    3332             : !        If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the
    3333             : !        algorithm enters the second stage.
    3334             : 
    3335        2377 :          ftest = finit + stp*gtest
    3336        2377 :          IF (stage == 1 .AND. f <= ftest .AND. g >= zero) &
    3337         550 :             stage = 2
    3338             : 
    3339             : !        Test for warnings.
    3340             : 
    3341        2377 :          IF (brackt .AND. (stp <= stmin .OR. stp >= stmax)) &
    3342           2 :             task = 'WARNING: ROUNDING ERRORS PREVENT PROGRESS'
    3343        2377 :          IF (brackt .AND. stmax - stmin <= xtol*stmax) &
    3344           2 :             task = 'WARNING: XTOL TEST SATISFIED'
    3345        2377 :          IF (stp == stpmax .AND. f <= ftest .AND. g <= gtest) &
    3346           8 :             task = 'WARNING: STP = STPMAX'
    3347        2377 :          IF (stp == stpmin .AND. (f > ftest .OR. g >= gtest)) &
    3348           0 :             task = 'WARNING: STP = STPMIN'
    3349             : 
    3350             : !        Test for convergence.
    3351             : 
    3352        2377 :          IF (f <= ftest .AND. ABS(g) <= gtol*(-ginit)) &
    3353        2041 :             task = 'CONVERGENCE'
    3354             : 
    3355             : !        Test for termination.
    3356             : 
    3357        2377 :          IF (.NOT. (task(1:4) == 'WARN' .OR. task(1:4) == 'CONV')) THEN
    3358             : 
    3359             : !        A modified function is used to predict the step during the
    3360             : !        first stage if a lower function value has been obtained but
    3361             : !        the decrease is not sufficient.
    3362             : 
    3363         329 :             IF (stage == 1 .AND. f <= fx .AND. f > ftest) THEN
    3364             : 
    3365             : !             Define the modified function and derivative values.
    3366             : 
    3367           0 :                fm = f - stp*gtest
    3368           0 :                fxm = fx - stx*gtest
    3369           0 :                fym = fy - sty*gtest
    3370           0 :                gm = g - gtest
    3371           0 :                gxm = gx - gtest
    3372           0 :                gym = gy - gtest
    3373             : 
    3374             : !             Call dcstep to update stx, sty, and to compute the new step.
    3375             : 
    3376             :                CALL dcstep(stx, fxm, gxm, sty, fym, gym, stp, fm, gm, &
    3377           0 :                            brackt, stmin, stmax)
    3378             : 
    3379             : !             Reset the function and derivative values for f.
    3380             : 
    3381           0 :                fx = fxm + stx*gtest
    3382           0 :                fy = fym + sty*gtest
    3383           0 :                gx = gxm + gtest
    3384           0 :                gy = gym + gtest
    3385             : 
    3386             :             ELSE
    3387             : 
    3388             : !             Call dcstep to update stx, sty, and to compute the new step.
    3389             : 
    3390             :                CALL dcstep(stx, fx, gx, sty, fy, gy, stp, f, g, &
    3391         329 :                            brackt, stmin, stmax)
    3392             : 
    3393             :             END IF
    3394             : 
    3395             : !          Decide if a bisection step is needed.
    3396             : 
    3397         329 :             IF (brackt) THEN
    3398         315 :                IF (ABS(sty - stx) >= p66*width1) stp = stx + p5*(sty - stx)
    3399             :                width1 = width
    3400             :                width = ABS(sty - stx)
    3401             :             END IF
    3402             : 
    3403             : !          Set the minimum and maximum steps allowed for stp.
    3404             : 
    3405         329 :             IF (brackt) THEN
    3406         315 :                stmin = MIN(stx, sty)
    3407         315 :                stmax = MAX(stx, sty)
    3408             :             ELSE
    3409          14 :                stmin = stp + xtrapl*(stp - stx)
    3410          14 :                stmax = stp + xtrapu*(stp - stx)
    3411             :             END IF
    3412             : 
    3413             : !          Force the step to be within the bounds stpmax and stpmin.
    3414             : 
    3415         329 :             stp = MAX(stp, stpmin)
    3416         329 :             stp = MIN(stp, stpmax)
    3417             : 
    3418             : !          If further progress is not possible, let stp be the best
    3419             : !          point obtained during the search.
    3420             : 
    3421             :             IF (brackt .AND. (stp <= stmin .OR. stp >= stmax) &
    3422         329 :                 .OR. (brackt .AND. stmax - stmin <= xtol*stmax)) stp = stx
    3423             : 
    3424             : !          Obtain another function and derivative.
    3425             : 
    3426         329 :             task = 'FG'
    3427             : 
    3428             :          END IF
    3429             :       END IF
    3430             : 
    3431             : !     Save local variables.
    3432             : 
    3433        4425 :       IF (brackt) THEN
    3434         613 :          isave(1) = 1
    3435             :       ELSE
    3436        3812 :          isave(1) = 0
    3437             :       END IF
    3438        4425 :       isave(2) = stage
    3439        4425 :       dsave(1) = ginit
    3440        4425 :       dsave(2) = gtest
    3441        4425 :       dsave(3) = gx
    3442        4425 :       dsave(4) = gy
    3443        4425 :       dsave(5) = finit
    3444        4425 :       dsave(6) = fx
    3445        4425 :       dsave(7) = fy
    3446        4425 :       dsave(8) = stx
    3447        4425 :       dsave(9) = sty
    3448        4425 :       dsave(10) = stmin
    3449        4425 :       dsave(11) = stmax
    3450        4425 :       dsave(12) = width
    3451        4425 :       dsave(13) = width1
    3452             : 
    3453        4425 :       RETURN
    3454        4425 :    END SUBROUTINE dcsrch
    3455             : 
    3456             : ! **************************************************************************************************
    3457             : !> \brief          This subroutine computes a safeguarded step for a search
    3458             : !>                 procedure and updates an interval that contains a step that
    3459             : !>                 satisfies a sufficient decrease and a curvature condition.
    3460             : !>
    3461             : !>                 The parameter stx contains the step with the least function
    3462             : !>                 value. If brackt is set to .true. then a minimizer has
    3463             : !>                 been bracketed in an interval with endpoints stx and sty.
    3464             : !>                 The parameter stp contains the current step.
    3465             : !>                 The subroutine assumes that if brackt is set to .true. then
    3466             : !>
    3467             : !>                    min(stx,sty) < stp < max(stx,sty),
    3468             : !>
    3469             : !>                 and that the derivative at stx is negative in the direction
    3470             : !>                 of the step.
    3471             : !> \param stx      On entry stx is the best step obtained so far and is an
    3472             : !>                 endpoint of the interval that contains the minimizer.
    3473             : !>                 On exit stx is the updated best step.
    3474             : !> \param fx       fx is the function at stx.
    3475             : !> \param dx       On entry dx is the derivative of the function at
    3476             : !>                 stx. The derivative must be negative in the direction of
    3477             : !>                 the step, that is, dx and stp - stx must have opposite
    3478             : !>                 signs.
    3479             : !>                 On exit dx is the derivative of the function at stx.
    3480             : !> \param sty      On entry sty is the second endpoint of the interval that
    3481             : !>                 contains the minimizer.
    3482             : !>                 On exit sty is the updated endpoint of the interval that
    3483             : !>                 contains the minimizer.
    3484             : !> \param fy       fy is the function at sty.
    3485             : !> \param dy       On entry dy is the derivative of the function at sty.
    3486             : !>                 On exit dy is the derivative of the function at the exit sty.
    3487             : !> \param stp      On entry stp is the current step. If brackt is set to .true.
    3488             : !>                 then on input stp must be between stx and sty.
    3489             : !>                 On exit stp is a new trial step.
    3490             : !> \param fp       fp is the function at stp
    3491             : !> \param dp_loc   dp_loc is the the derivative of the function at stp.
    3492             : !> \param brackt   On entry brackt specifies if a minimizer has been bracketed.
    3493             : !>                 Initially brackt must be set to .false.
    3494             : !>                 On exit brackt specifies if a minimizer has been bracketed.
    3495             : !>                 When a minimizer is bracketed brackt is set to .true.
    3496             : !> \param stpmin   stpmin is a lower bound for the step.
    3497             : !> \param stpmax   stpmax is an upper bound for the step.
    3498             : ! **************************************************************************************************
    3499         329 :    SUBROUTINE dcstep(stx, fx, dx, sty, fy, dy, stp, fp, dp_loc, brackt, &
    3500             :                      stpmin, stpmax)
    3501             :       REAL(KIND=dp), INTENT(inout)                       :: stx, fx, dx, sty, fy, dy, stp
    3502             :       REAL(KIND=dp), INTENT(in)                          :: fp, dp_loc
    3503             :       LOGICAL, INTENT(inout)                             :: brackt
    3504             :       REAL(KIND=dp), INTENT(in)                          :: stpmin, stpmax
    3505             : 
    3506             :       REAL(KIND=dp), PARAMETER                           :: p66 = 0.66_dp, three = 3.0_dp, &
    3507             :                                                             two = 2.0_dp, zero = 0.0_dp
    3508             : 
    3509             :       REAL(KIND=dp)                                      :: gamma, p, q, r, s, sgnd, stpc, stpf, &
    3510             :                                                             stpq, theta
    3511             : 
    3512             : !
    3513             : !     MINPACK-1 Project. June 1983
    3514             : !     Argonne National Laboratory.
    3515             : !     Jorge J. More' and David J. Thuente.
    3516             : !
    3517             : !     MINPACK-2 Project. October 1993.
    3518             : !     Argonne National Laboratory and University of Minnesota.
    3519             : !     Brett M. Averick and Jorge J. More'.
    3520             : !
    3521             : !     **********
    3522             : 
    3523         329 :       sgnd = dp_loc*SIGN(1.0_dp, dx)
    3524             : 
    3525             : !     First case: A higher function value. The minimum is bracketed.
    3526             : !     If the cubic step is closer to stx than the quadratic step, the
    3527             : !     cubic step is taken, otherwise the average of the cubic and
    3528             : !     quadratic steps is taken.
    3529             : 
    3530         329 :       IF (fp > fx) THEN
    3531         291 :          theta = three*(fx - fp)/(stp - stx) + dx + dp_loc
    3532         291 :          s = MAX(ABS(theta), ABS(dx), ABS(dp_loc))
    3533         291 :          gamma = s*SQRT((theta/s)**2 - (dx/s)*(dp_loc/s))
    3534         291 :          IF (stp < stx) gamma = -gamma
    3535         291 :          p = (gamma - dx) + theta
    3536         291 :          q = ((gamma - dx) + gamma) + dp_loc
    3537         291 :          r = p/q
    3538         291 :          stpc = stx + r*(stp - stx)
    3539             :          stpq = stx + ((dx/((fx - fp)/(stp - stx) + dx))/two)*          &
    3540         291 :      &                                                       (stp - stx)
    3541         291 :          IF (ABS(stpc - stx) < ABS(stpq - stx)) THEN
    3542             :             stpf = stpc
    3543             :          ELSE
    3544         123 :             stpf = stpc + (stpq - stpc)/two
    3545             :          END IF
    3546         291 :          brackt = .TRUE.
    3547             : 
    3548             : !     Second case: A lower function value and derivatives of opposite
    3549             : !     sign. The minimum is bracketed. If the cubic step is farther from
    3550             : !     stp than the secant step, the cubic step is taken, otherwise the
    3551             : !     secant step is taken.
    3552             : 
    3553          38 :       ELSE IF (sgnd < zero) THEN
    3554          17 :          theta = three*(fx - fp)/(stp - stx) + dx + dp_loc
    3555          17 :          s = MAX(ABS(theta), ABS(dx), ABS(dp_loc))
    3556          17 :          gamma = s*SQRT((theta/s)**2 - (dx/s)*(dp_loc/s))
    3557          17 :          IF (stp > stx) gamma = -gamma
    3558          17 :          p = (gamma - dp_loc) + theta
    3559          17 :          q = ((gamma - dp_loc) + gamma) + dx
    3560          17 :          r = p/q
    3561          17 :          stpc = stp + r*(stx - stp)
    3562          17 :          stpq = stp + (dp_loc/(dp_loc - dx))*(stx - stp)
    3563          17 :          IF (ABS(stpc - stp) > ABS(stpq - stp)) THEN
    3564             :             stpf = stpc
    3565             :          ELSE
    3566          14 :             stpf = stpq
    3567             :          END IF
    3568          17 :          brackt = .TRUE.
    3569             : 
    3570             : !     Third case: A lower function value, derivatives of the same sign,
    3571             : !     and the magnitude of the derivative decreases.
    3572             : 
    3573          21 :       ELSE IF (ABS(dp_loc) < ABS(dx)) THEN
    3574             : 
    3575             : !        The cubic step is computed only if the cubic tends to infinity
    3576             : !        in the direction of the step or if the minimum of the cubic
    3577             : !        is beyond stp. Otherwise the cubic step is defined to be the
    3578             : !        secant step.
    3579             : 
    3580          10 :          theta = three*(fx - fp)/(stp - stx) + dx + dp_loc
    3581          10 :          s = MAX(ABS(theta), ABS(dx), ABS(dp_loc))
    3582             : 
    3583             : !        The case gamma = 0 only arises if the cubic does not tend
    3584             : !        to infinity in the direction of the step.
    3585             : 
    3586          10 :          gamma = s*SQRT(MAX(zero, (theta/s)**2 - (dx/s)*(dp_loc/s)))
    3587          10 :          IF (stp > stx) gamma = -gamma
    3588          10 :          p = (gamma - dp_loc) + theta
    3589          10 :          q = (gamma + (dx - dp_loc)) + gamma
    3590          10 :          r = p/q
    3591          10 :          IF (r < zero .AND. gamma /= zero) THEN
    3592          10 :             stpc = stp + r*(stx - stp)
    3593           0 :          ELSE IF (stp > stx) THEN
    3594           0 :             stpc = stpmax
    3595             :          ELSE
    3596           0 :             stpc = stpmin
    3597             :          END IF
    3598          10 :          stpq = stp + (dp_loc/(dp_loc - dx))*(stx - stp)
    3599             : 
    3600          10 :          IF (brackt) THEN
    3601             : 
    3602             : !           A minimizer has been bracketed. If the cubic step is
    3603             : !           closer to stp than the secant step, the cubic step is
    3604             : !           taken, otherwise the secant step is taken.
    3605             : 
    3606           5 :             IF (ABS(stpc - stp) < ABS(stpq - stp)) THEN
    3607             :                stpf = stpc
    3608             :             ELSE
    3609           0 :                stpf = stpq
    3610             :             END IF
    3611           5 :             IF (stp > stx) THEN
    3612           5 :                stpf = MIN(stp + p66*(sty - stp), stpf)
    3613             :             ELSE
    3614           0 :                stpf = MAX(stp + p66*(sty - stp), stpf)
    3615             :             END IF
    3616             :          ELSE
    3617             : 
    3618             : !           A minimizer has not been bracketed. If the cubic step is
    3619             : !           farther from stp than the secant step, the cubic step is
    3620             : !           taken, otherwise the secant step is taken.
    3621             : 
    3622           5 :             IF (ABS(stpc - stp) > ABS(stpq - stp)) THEN
    3623             :                stpf = stpc
    3624             :             ELSE
    3625           4 :                stpf = stpq
    3626             :             END IF
    3627           5 :             stpf = MIN(stpmax, stpf)
    3628           5 :             stpf = MAX(stpmin, stpf)
    3629             :          END IF
    3630             : 
    3631             : !     Fourth case: A lower function value, derivatives of the same sign,
    3632             : !     and the magnitude of the derivative does not decrease. If the
    3633             : !     minimum is not bracketed, the step is either stpmin or stpmax,
    3634             : !     otherwise the cubic step is taken.
    3635             : 
    3636             :       ELSE
    3637          11 :          IF (brackt) THEN
    3638           2 :             theta = three*(fp - fy)/(sty - stp) + dy + dp_loc
    3639           2 :             s = MAX(ABS(theta), ABS(dy), ABS(dp_loc))
    3640           2 :             gamma = s*SQRT((theta/s)**2 - (dy/s)*(dp_loc/s))
    3641           2 :             IF (stp > sty) gamma = -gamma
    3642           2 :             p = (gamma - dp_loc) + theta
    3643           2 :             q = ((gamma - dp_loc) + gamma) + dy
    3644           2 :             r = p/q
    3645           2 :             stpc = stp + r*(sty - stp)
    3646           2 :             stpf = stpc
    3647           9 :          ELSE IF (stp > stx) THEN
    3648           9 :             stpf = stpmax
    3649             :          ELSE
    3650           0 :             stpf = stpmin
    3651             :          END IF
    3652             :       END IF
    3653             : 
    3654             : !     Update the interval which contains a minimizer.
    3655             : 
    3656         329 :       IF (fp > fx) THEN
    3657         291 :          sty = stp
    3658         291 :          fy = fp
    3659         291 :          dy = dp_loc
    3660             :       ELSE
    3661          38 :          IF (sgnd < zero) THEN
    3662          17 :             sty = stx
    3663          17 :             fy = fx
    3664          17 :             dy = dx
    3665             :          END IF
    3666          38 :          stx = stp
    3667          38 :          fx = fp
    3668          38 :          dx = dp_loc
    3669             :       END IF
    3670             : 
    3671             : !     Compute the new step.
    3672             : 
    3673         329 :       stp = stpf
    3674             : 
    3675         329 :       RETURN
    3676             :    END SUBROUTINE dcstep
    3677             : 
    3678             : !MK LINPACK
    3679             : 
    3680             : ! **************************************************************************************************
    3681             : !> \brief         factors a double precision symmetric positive definite
    3682             : !>                matrix.
    3683             : !>
    3684             : !>                dpofa is usually called by dpoco, but it can be called
    3685             : !>                directly with a saving in time if  rcond  is not needed.
    3686             : !>                (time for dpoco) = (1 + 18/n)*(time for dpofa) .
    3687             : !> \param a       the symmetric matrix to be factored.  only the
    3688             : !>                diagonal and upper triangle are used.
    3689             : !>                on return
    3690             : !>                an upper triangular matrix  r  so that  a = trans(r)*r
    3691             : !>                where  trans(r)  is the transpose.
    3692             : !>                the strict lower triangle is unaltered.
    3693             : !>                if  info .ne. 0 , the factorization is not complete.
    3694             : !> \param lda     the leading dimension of the array  a .
    3695             : !> \param n       the order of the matrix  a .
    3696             : !> \param info    = 0  for normal return.
    3697             : !>                = k  signals an error condition.  the leading minor
    3698             : !>                     of order  k  is not positive definite.
    3699             : ! **************************************************************************************************
    3700        5559 :    SUBROUTINE dpofa(a, lda, n, info)
    3701             :       INTEGER, INTENT(in)                                :: lda
    3702             :       REAL(KIND=dp)                                      :: a(lda, *)
    3703             :       INTEGER, INTENT(in)                                :: n
    3704             :       INTEGER                                            :: info
    3705             : 
    3706             :       INTEGER                                            :: j, jm1, k
    3707             :       REAL(KIND=dp)                                      :: ddot, s, t
    3708             : 
    3709             : !
    3710             : !     linpack.  this version dated 08/14/78 .
    3711             : !     cleve moler, university of new mexico, argonne national lab.
    3712             : !
    3713             : !     begin block with ...exits to 40
    3714             : !
    3715             : !
    3716             : 
    3717       31830 :       DO j = 1, n
    3718       26271 :          info = j
    3719       26271 :          s = 0.0_dp
    3720       26271 :          jm1 = j - 1
    3721       26271 :          IF (.NOT. (jm1 < 1)) THEN
    3722       82341 :             DO k = 1, jm1
    3723       61629 :                t = a(k, j) - ddot(k - 1, a(1, k), 1, a(1, j), 1)
    3724       61629 :                t = t/a(k, k)
    3725       61629 :                a(k, j) = t
    3726       82341 :                s = s + t*t
    3727             :             END DO
    3728             :          END IF
    3729       26271 :          s = a(j, j) - s
    3730             : !     ......exit
    3731       26271 :          IF (s <= 0.0_dp) EXIT
    3732       26271 :          a(j, j) = SQRT(s)
    3733       31830 :          info = 0
    3734             :       END DO
    3735        5559 :       RETURN
    3736             :    END SUBROUTINE dpofa
    3737             : 
    3738             : ! **************************************************************************************************
    3739             : !> \brief           dtrsl solves systems of the form
    3740             : !>
    3741             : !>                  t * x = b
    3742             : !>                  or
    3743             : !>                  trans(t) * x = b
    3744             : !>
    3745             : !>                  where t is a triangular matrix of order n. here trans(t)
    3746             : !>                  denotes the transpose of the matrix t.
    3747             : !> \param t         t contains the matrix of the system. the zero
    3748             : !>                  elements of the matrix are not referenced, and
    3749             : !>                  the corresponding elements of the array can be
    3750             : !>                  used to store other information.
    3751             : !> \param ldt       ldt is the leading dimension of the array t.
    3752             : !> \param n         n is the order of the system.
    3753             : !> \param b         contains the right hand side of the system.
    3754             : !>                  on return
    3755             : !>                  b contains the solution, if info .eq. 0.
    3756             : !>                  otherwise b is unaltered.
    3757             : !> \param job       job specifies what kind of system is to be solved.
    3758             : !>                   if job is
    3759             : !>                       00   solve t*x=b, t lower triangular,
    3760             : !>                       01   solve t*x=b, t upper triangular,
    3761             : !>                       10   solve trans(t)*x=b, t lower triangular,
    3762             : !>                       11   solve trans(t)*x=b, t upper triangular.
    3763             : !> \param info      on return
    3764             : !>                  info contains zero if the system is nonsingular.
    3765             : !>                  otherwise info contains the index of
    3766             : !>                  the first zero diagonal element of t.
    3767             : ! **************************************************************************************************
    3768       12507 :    SUBROUTINE dtrsl(t, ldt, n, b, job, info)
    3769             :       INTEGER, INTENT(in)                                :: ldt
    3770             :       REAL(KIND=dp), INTENT(in)                          :: t(ldt, *)
    3771             :       INTEGER, INTENT(in)                                :: n
    3772             :       REAL(KIND=dp), INTENT(inout)                       :: b(*)
    3773             :       INTEGER, INTENT(in)                                :: job
    3774             :       INTEGER, INTENT(out)                               :: info
    3775             : 
    3776             :       INTEGER                                            :: CASE, j, jj
    3777             :       REAL(KIND=dp)                                      :: ddot, temp
    3778             : 
    3779             : !     linpack. this version dated 08/14/78 .
    3780             : !     g. w. stewart, university of maryland, argonne national lab.
    3781             : !
    3782             : !     begin block permitting ...exits to 150
    3783             : !
    3784             : !        check for zero diagonal elements.
    3785             : !
    3786             : 
    3787       97502 :       DO info = 1, n
    3788             : !     ......exit
    3789       97502 :          IF (t(info, info) == 0.0_dp) RETURN
    3790             :       END DO
    3791       12507 :       info = 0
    3792             : !
    3793             : !        determine the task and go to it.
    3794             : !
    3795       12507 :       CASE = 1
    3796       12507 :       IF (MOD(job, 10) /= 0) CASE = 2
    3797       12507 :       IF (MOD(job, 100)/10 /= 0) CASE = CASE + 2
    3798             : 
    3799           0 :       SELECT CASE (CASE)
    3800             :       CASE (1)
    3801             : !
    3802             : !        solve t*x=b for t lower triangular
    3803             : !
    3804           0 :          b(1) = b(1)/t(1, 1)
    3805           0 :          IF (n > 1) THEN
    3806           0 :             DO j = 2, n
    3807           0 :                temp = -b(j - 1)
    3808           0 :                CALL daxpy(n - j + 1, temp, t(j, j - 1), 1, b(j), 1)
    3809           0 :                b(j) = b(j)/t(j, j)
    3810             :             END DO
    3811             :          END IF
    3812             :       CASE (2)
    3813             : !
    3814             : !        solve t*x=b for t upper triangular.
    3815             : !
    3816        1875 :          b(n) = b(n)/t(n, n)
    3817        1875 :          IF (n > 1) THEN
    3818       17570 :             DO jj = 2, n
    3819       15701 :                j = n - jj + 1
    3820       15701 :                temp = -b(j + 1)
    3821       15701 :                CALL daxpy(j, temp, t(1, j + 1), 1, b(1), 1)
    3822       17570 :                b(j) = b(j)/t(j, j)
    3823             :             END DO
    3824             :          END IF
    3825             :       CASE (3)
    3826             : !
    3827             : !        solve trans(t)*x=b for t lower triangular.
    3828             : !
    3829           0 :          b(n) = b(n)/t(n, n)
    3830           0 :          IF (n > 1) THEN
    3831           0 :             DO jj = 2, n
    3832           0 :                j = n - jj + 1
    3833           0 :                b(j) = b(j) - ddot(jj - 1, t(j + 1, j), 1, b(j + 1), 1)
    3834           0 :                b(j) = b(j)/t(j, j)
    3835             :             END DO
    3836             :          END IF
    3837             :       CASE (4)
    3838             : !
    3839             : !        solve trans(t)*x=b for t upper triangular.
    3840             : !
    3841       10632 :          b(1) = b(1)/t(1, 1)
    3842       10632 :          IF (.NOT. (n < 2)) THEN
    3843       67281 :             DO j = 2, n
    3844       56787 :                b(j) = b(j) - ddot(j - 1, t(1, j), 1, b(1), 1)
    3845       67281 :                b(j) = b(j)/t(j, j)
    3846             :             END DO
    3847             :          END IF
    3848             :       CASE DEFAULT
    3849       12507 :          CPABORT("unexpected case")
    3850             :       END SELECT
    3851             : 
    3852             :       RETURN
    3853             :    END SUBROUTINE dtrsl
    3854             : 
    3855             : !MK Timer
    3856             : 
    3857             : ! **************************************************************************************************
    3858             : !> \brief This routine computes cpu time in double precision; it makes use o
    3859             : !>        the intrinsic f90 cpu_time therefore a conversion type is
    3860             : !>        needed.
    3861             : !> \param ttime ...
    3862             : ! **************************************************************************************************
    3863        8402 :    SUBROUTINE timer(ttime)
    3864             :       REAL(KIND=dp)                                      :: ttime
    3865             : 
    3866             : !
    3867             : !     REAL temp
    3868             : !
    3869             : !           J.L Morales  Departamento de Matematicas,
    3870             : !                        Instituto Tecnologico Autonomo de Mexico
    3871             : !                        Mexico D.F.
    3872             : !
    3873             : !           J.L Nocedal  Department of Electrical Engineering and
    3874             : !                        Computer Science.
    3875             : !                        Northwestern University. Evanston, IL. USA
    3876             : !
    3877             : !                        January 21, 2011
    3878             : !
    3879             : !MK      temp = sngl(ttime)
    3880             : !MK      CALL cpu_time(temp)
    3881             : !MK      ttime = REAL(temp, KIND=dp)
    3882             : 
    3883        8402 :       ttime = m_walltime()
    3884             : 
    3885        8402 :    END SUBROUTINE timer
    3886             : 
    3887             : ! **************************************************************************************************
    3888             : !> \brief  Saves the lcoal variables, long term this should be replaces by a lbfgs type
    3889             : !> \param lsave    lsave is a working array
    3890             : !>                 On exit with 'task' = NEW_X, the following information is available:
    3891             : !>                 If lsave(1) = .true.  then  the initial X has been replaced by
    3892             : !>                               its projection in the feasible set
    3893             : !>                 If lsave(2) = .true.  then  the problem is constrained;
    3894             : !>                 If lsave(3) = .true.  then  each variable has upper and lower bounds;
    3895             : !> \param isave    isave is a working array
    3896             : !>                 On exit with 'task' = NEW_X, the following information is available:
    3897             : !>                 isave(22) = the total number of intervals explored in the
    3898             : !>                         search of Cauchy points;
    3899             : !>                 isave(26) = the total number of skipped BFGS updates before the current iteration;
    3900             : !>                 isave(30) = the number of current iteration;
    3901             : !>                 isave(31) = the total number of BFGS updates prior the current iteration;
    3902             : !>                 isave(33) = the number of intervals explored in the search of
    3903             : !>                             Cauchy point in the current iteration;
    3904             : !>                 isave(34) = the total number of function and gradient evaluations;
    3905             : !>                 isave(36) = the number of function value or gradient
    3906             : !>                             evaluations in the current iteration;
    3907             : !>                 if isave(37) = 0  then the subspace argmin is within the box;
    3908             : !>                 if isave(37) = 1  then the subspace argmin is beyond the box;
    3909             : !>                 isave(38) = the number of free variables in the current iteration;
    3910             : !>                 isave(39) = the number of active constraints in the current iteration;
    3911             : !>                 n + 1 - isave(40) = the number of variables leaving the set of
    3912             : !>                                     active constraints in the current iteration;
    3913             : !>                 isave(41) = the number of variables entering the set of active
    3914             : !>                             constraints in the current iteration.
    3915             : !> \param dsave    dsave is a working array of dimension 29.
    3916             : !>                 On exit with 'task' = NEW_X, the following information is available:
    3917             : !>                 dsave(1) = current 'theta' in the BFGS matrix;
    3918             : !>                 dsave(2) = f(x) in the previous iteration;
    3919             : !>                 dsave(3) = factr*epsmch;
    3920             : !>                 dsave(4) = 2-norm of the line search direction vector;
    3921             : !>                 dsave(5) = the machine precision epsmch generated by the code;
    3922             : !>                 dsave(7) = the accumulated time spent on searching for Cauchy points;
    3923             : !>                 dsave(8) = the accumulated time spent on subspace minimization;
    3924             : !>                 dsave(9) = the accumulated time spent on line search;
    3925             : !>                 dsave(11) = the slope of the line search function at the current point of line search;
    3926             : !>                 dsave(12) = the maximum relative step length imposed in line search;
    3927             : !>                 dsave(13) = the infinity norm of the projected gradient;
    3928             : !>                 dsave(14) = the relative step length in the line search;
    3929             : !>                 dsave(15) = the slope of the line search function at the starting point of the line search;
    3930             : !>                 dsave(16) = the square of the 2-norm of the line search direction vector.
    3931             : !> \param x_projected ...
    3932             : !> \param constrained ...
    3933             : !> \param boxed ...
    3934             : !> \param updatd ...
    3935             : !> \param nintol ...
    3936             : !> \param itfile ...
    3937             : !> \param iback ...
    3938             : !> \param nskip ...
    3939             : !> \param head ...
    3940             : !> \param col ...
    3941             : !> \param itail ...
    3942             : !> \param iter ...
    3943             : !> \param iupdat ...
    3944             : !> \param nseg ...
    3945             : !> \param nfgv ...
    3946             : !> \param info ...
    3947             : !> \param ifun ...
    3948             : !> \param iword ...
    3949             : !> \param nfree ...
    3950             : !> \param nact ...
    3951             : !> \param ileave ...
    3952             : !> \param nenter ...
    3953             : !> \param theta ...
    3954             : !> \param fold ...
    3955             : !> \param tol ...
    3956             : !> \param dnorm ...
    3957             : !> \param epsmch ...
    3958             : !> \param cpu1 ...
    3959             : !> \param cachyt ...
    3960             : !> \param sbtime ...
    3961             : !> \param lnscht ...
    3962             : !> \param time1 ...
    3963             : !> \param gd ...
    3964             : !> \param step_max ...
    3965             : !> \param g_inf_norm ...
    3966             : !> \param stp ...
    3967             : !> \param gdold ...
    3968             : !> \param dtd ...
    3969             : !> \author Samuel Andermatt (01.15)
    3970             : ! **************************************************************************************************
    3971             : 
    3972        4617 :    SUBROUTINE save_local(lsave,isave,dsave,x_projected,constrained,boxed,updatd,nintol,itfile,iback,nskip,head,col,itail,&
    3973             :                   iter, iupdat, nseg, nfgv, info, ifun, iword, nfree, nact, ileave, nenter, theta, fold, tol, dnorm, epsmch, cpu1, &
    3974             :                          cachyt, sbtime, lnscht, time1, gd, step_max, g_inf_norm, stp, gdold, dtd)
    3975             :       LOGICAL, INTENT(out)                               :: lsave(4)
    3976             :       INTEGER, INTENT(out)                               :: isave(23)
    3977             :       REAL(KIND=dp), INTENT(out)                         :: dsave(29)
    3978             :       LOGICAL, INTENT(in)                                :: x_projected, constrained, boxed, updatd
    3979             :       INTEGER, INTENT(in)                                :: nintol, itfile, iback, nskip, head, col, &
    3980             :                                                             itail, iter, iupdat, nseg, nfgv, info, &
    3981             :                                                             ifun, iword, nfree, nact, ileave, &
    3982             :                                                             nenter
    3983             :       REAL(KIND=dp), INTENT(in)                          :: theta, fold, tol, dnorm, epsmch, cpu1, &
    3984             :                                                             cachyt, sbtime, lnscht, time1, gd, &
    3985             :                                                             step_max, g_inf_norm, stp, gdold, dtd
    3986             : 
    3987        4617 :       lsave(1) = x_projected
    3988        4617 :       lsave(2) = constrained
    3989        4617 :       lsave(3) = boxed
    3990        4617 :       lsave(4) = updatd
    3991             : 
    3992        4617 :       isave(1) = nintol
    3993        4617 :       isave(3) = itfile
    3994        4617 :       isave(4) = iback
    3995        4617 :       isave(5) = nskip
    3996        4617 :       isave(6) = head
    3997        4617 :       isave(7) = col
    3998        4617 :       isave(8) = itail
    3999        4617 :       isave(9) = iter
    4000        4617 :       isave(10) = iupdat
    4001        4617 :       isave(12) = nseg
    4002        4617 :       isave(13) = nfgv
    4003        4617 :       isave(14) = info
    4004        4617 :       isave(15) = ifun
    4005        4617 :       isave(16) = iword
    4006        4617 :       isave(17) = nfree
    4007        4617 :       isave(18) = nact
    4008        4617 :       isave(19) = ileave
    4009        4617 :       isave(20) = nenter
    4010             : 
    4011        4617 :       dsave(1) = theta
    4012        4617 :       dsave(2) = fold
    4013        4617 :       dsave(3) = tol
    4014        4617 :       dsave(4) = dnorm
    4015        4617 :       dsave(5) = epsmch
    4016        4617 :       dsave(6) = cpu1
    4017        4617 :       dsave(7) = cachyt
    4018        4617 :       dsave(8) = sbtime
    4019        4617 :       dsave(9) = lnscht
    4020        4617 :       dsave(10) = time1
    4021        4617 :       dsave(11) = gd
    4022        4617 :       dsave(12) = step_max
    4023        4617 :       dsave(13) = g_inf_norm
    4024        4617 :       dsave(14) = stp
    4025        4617 :       dsave(15) = gdold
    4026        4617 :       dsave(16) = dtd
    4027             : 
    4028        4617 :    END SUBROUTINE save_local
    4029             : 
    4030             : END MODULE cp_lbfgs

Generated by: LCOV version 1.15