MODULE tron

! ************************************************************************* !
!                                                                           !
!            COPYRIGHT NOTIFICATION                                         !
!                                                                           !
! This program discloses material protectable under copyright laws of       !
! the United States. Permission to copy and modify this software and its    !
! documentation for internal research use is hereby granted, provided       !
! that this notice is retained thereon and on all copies or modifications.  !
! The University of Chicago makes no representations as to the suitability  !
! and operability of this software for any purpose.                         !
! It is provided "as is" without express or implied warranty.               !
!                                                                           !
! Use of this software for commercial purposes is expressly prohibited      !
! without contacting                                                        !
!                                                                           !
!    Jorge J. More'                                                         !
!    Mathematics and Computer Science Division                              !
!    Argonne National Laboratory                                            !
!    9700 S. Cass Ave.                                                      !
!    Argonne, Illinois 60439-4844                                           !
!    e-mail: more@mcs.anl.gov                                               !
!                                                                           !
! Argonne National Laboratory with facilities in the states of              !
! Illinois and Idaho, is owned by The United States Government, and         !
! operated by the University of Chicago under provision of a contract       !
! with the Department of Energy.                                            !
!                                                                           !
! ************************************************************************* !
!                                                                           !
!            ADDITIONAL INFORMATION                                         !
!                                                                           !
! Chih-Jen Lin and Jorge J. More',                                          !
! Newton's method for large bound-constrained optimization problems,        !
! Argonne National Laboratory,                                              !
! Mathematics and Computer Science Division,                                !
! Preprint ANL/MCS-P724-0898,                                               !
! August 1998 (Revised March 1999).                                         !
!                                                                           !
! http://www.mcs.anl.gov/~more/papers/tron.ps.gz                            !
!                                                                           !
! The Fortran 77 code can be downloaded from:                               !
!                                                                           !
! http://www-unix.mcs.anl.gov/~more/tron                                    !
!                                                                           !
! ************************************************************************* !
!                                                                           !
! Last modification (Fortran 77 version): April 29, 1999                    !


IMPLICIT NONE
INTEGER, PARAMETER                     :: dp = SELECTED_REAL_KIND(14, 60)
REAL (dp), PARAMETER, PRIVATE          :: zero = 0.0_dp, one = 1.0_dp
REAL (dp), ALLOCATABLE, SAVE, PRIVATE  :: xc(:), s(:), dsave(:), wa(:)
INTEGER, ALLOCATABLE, SAVE, PRIVATE    :: indfree(:), isave(:)

CONTAINS


SUBROUTINE dtron(n, x, xl, xu, f, g, a, adiag, acol_ptr, arow_ind,  &
                 frtol, fatol, fmin, cgtol, itermax, delta, task,  &
                 b, bdiag, bcol_ptr, brow_ind,   &
                 l, ldiag, lcol_ptr, lrow_ind, iterscg)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 11:18:32
! Latest revision - 5 July 1999

INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: x(:)
REAL (dp), INTENT(IN)      :: xl(:)
REAL (dp), INTENT(IN)      :: xu(:)
REAL (dp), INTENT(IN OUT)  :: f
REAL (dp), INTENT(IN)      :: g(:)
REAL (dp), INTENT(IN)      :: a(:)
REAL (dp), INTENT(IN)      :: adiag(:)
INTEGER, INTENT(IN)        :: acol_ptr(:)    ! acol_ptr(n+1)
INTEGER, INTENT(IN)        :: arow_ind(:)
REAL (dp), INTENT(IN)      :: frtol
REAL (dp), INTENT(IN)      :: fatol
REAL (dp), INTENT(IN)      :: fmin
REAL (dp), INTENT(IN)      :: cgtol
INTEGER, INTENT(IN)        :: itermax
REAL (dp), INTENT(IN OUT)  :: delta
CHARACTER (LEN=*), INTENT(IN OUT) :: task
REAL (dp), INTENT(OUT)     :: b(:)
REAL (dp), INTENT(OUT)     :: bdiag(:)
INTEGER, INTENT(OUT)       :: bcol_ptr(:)    ! bcol_ptr(n+1)
INTEGER, INTENT(OUT)       :: brow_ind(:)
REAL (dp), INTENT(OUT)     :: l(:)
REAL (dp), INTENT(OUT)     :: ldiag(:)
INTEGER, INTENT(OUT)       :: lcol_ptr(:)    ! lcol_ptr(n+1)
INTEGER, INTENT(OUT)       :: lrow_ind(:)
INTEGER, INTENT(OUT)       :: iterscg

!  *********

!  Subroutine dtron

!  This subroutine implements a trust region Newton method for the
!  solution of large bound-constrained optimization problems

!        min { f(x) : xl <= x <= xu }

!  where the Hessian matrix is sparse. The user must evaluate the
!  function, gradient, and the Hessian matrix.

!  This subroutine uses reverse communication.
!  The user must choose an initial approximation x to the minimizer,
!  and make an initial call with task set to 'START'.
!  On exit task indicates the required action.

!  A typical invocation has the following outline:

!  Compute a starting vector x.
!  Compute the sparsity pattern of the Hessian matrix and
!  store in compressed column storage in (acol_ptr,arow_ind).

!  task = 'START'
!  do while (search)

!     if (task .eq. 'F' .or. task .eq. 'START') then
!        Evaluate the function at x and store in f.
!     end if
!     if (task .eq. 'GH' .or. task .eq. 'START') then
!        Evaluate the gradient at x and store in g.
!        Evaluate the Hessian at x and store in compressed
!        column storage in (a, adiag, acol_ptr, arow_ind)
!     end if

!     call dtron(n, x, xl, xu, f, g, a, adiag, acol_ptr, arow_ind,
!                frtol, fatol, fmin, cgtol, itermax, delta, task,
!                b, bdiag, bcol_ptr, brow_ind,
!                l, ldiag, lcol_ptr, lrow_ind, iterscg)

!     if (task(1:4) .eq. 'CONV') search = .false.

!   end do

!  NOTE: The user must not alter work arrays between calls.

!  The subroutine statement is

!    subroutine dtron(n, x, xl, xu, f, g, a, adiag, acol_ptr, arow_ind,
!                     frtol, fatol, fmin, cgtol, itermax, delta, task,
!                     b, bdiag, bcol_ptr, brow_ind,
!                     l, ldiag, lcol_ptr, lrow_ind, isave)

!  where

!    n is an integer variable.
!      On entry n is the number of variables.
!      On exit n is unchanged.

!    x is a REAL (dp) array of dimension n.
!      On entry x specifies the vector x.
!      On exit x is the final minimizer.

!    xl is a REAL (dp) array of dimension n.
!      On entry xl is the vector of lower bounds.
!      On exit xl is unchanged.

!    xu is a REAL (dp) array of dimension n.
!      On entry xu is the vector of upper bounds.
!      On exit xu is unchanged.

!    f is a REAL (dp) variable.
!      On entry f must contain the function at x.
!      On exit f is unchanged.   This is NOT true.

!    g is a REAL (dp) array of dimension n.
!      On entry g must contain the gradient at x.
!      On exit g is unchanged.

!    a is a REAL (dp) array of dimension nnz.
!      On entry a must contain the strict lower triangular part
!         of A in compressed column storage.
!      On exit a is unchanged.

!    adiag is a REAL (dp) array of dimension n.
!      On entry adiag must contain the diagonal elements of A.
!      On exit adiag is unchanged.

!    acol_ptr is an integer array of dimension n + 1.
!      On entry acol_ptr must contain pointers to the columns of A.
!         The nonzeros in column j of A must be in positions
!         acol_ptr(j), ... , acol_ptr(j+1) - 1.
!      On exit acol_ptr is unchanged.

!    arow_ind is an integer array of dimension nnz.
!      On entry arow_ind must contain row indices for the strict
!         lower triangular part of A in compressed column storage.
!      On exit arow_ind is unchanged.

!    frtol is a REAL (dp) variable.
!      On entry frtol specifies the relative error desired in the function.
!         Convergence occurs if the estimate of the relative error between f(x)
!         and f(xsol), where xsol is a local minimizer, is less than frtol.
!      On exit frtol is unchanged.

!    fatol is a REAL (dp) variable.
!      On entry fatol specifies the absolute error desired in the function.
!         Convergence occurs if the estimate of the absolute error between f(x)
!         and f(xsol), where xsol is a local minimizer, is less than fatol.
!      On exit fatol is unchanged.

!    fmin is a REAL (dp) variable.
!      On entry fmin specifies a lower bound for the function.
!         The subroutine exits with a warning if f < fmin.
!      On exit fmin is unchanged.

!    cgtol is a REAL (dp) variable.
!      On entry cgtol specifies the convergence criteria for
!         the conjugate gradient method.
!      On exit cgtol is unchanged.

!    itermax is an integer variable.
!      On entry itermax specifies the limit on the number of
!         conjugate gradient iterations.
!      On exit itermax is unchanged.

!    delta is a REAL (dp) variable.
!      On entry delta is the trust region bound.
!      On exit delta is unchanged.  This is NOT true.

!    task is a character variable of length at least 60.
!      On initial entry task must be set to 'START'.
!      On exit task indicates the required action:

!         If task(1:1) = 'F' then evaluate the function at x.

!         If task(1:2) = 'GH' then evaluate the gradient and the
!         Hessian matrix at x.

!         If task(1:4) = 'CONV' then the search is successful.

!         If task(1:4) = 'WARN' then the subroutine is not able
!         to satisfy the convergence conditions. The exit value
!         of x contains the best approximation found.

!    bdiag is a REAL (dp) array of dimension n.
!      On entry bdiag need not be specified.
!      On exit bdiag contains the diagonal elements of B.

!    bcol_ptr is an integer array of dimension n + 1.
!      On entry bcol_ptr need not be specified
!      On exit bcol_ptr contains pointers to the columns of B.
!         The nonzeros in column j of B are in the
!         bcol_ptr(j), ... , bcol_ptr(j+1) - 1 positions of b.

!    brow_ind is an integer array of dimension nnz.
!      On entry brow_ind need not be specified.
!      On exit brow_ind contains row indices for the strict lower
!         triangular part of B in compressed column storage.

!    l is a REAL (dp) array of dimension nnz + n*p.
!      On entry l need not be specified.
!      On exit l contains the strict lower triangular part
!         of L in compressed column storage.

!    ldiag is a REAL (dp) array of dimension n.
!      On entry ldiag need not be specified.
!      On exit ldiag contains the diagonal elements of L.

!    lcol_ptr is an integer array of dimension n + 1.
!      On entry lcol_ptr need not be specified.
!      On exit lcol_ptr contains pointers to the columns of L.
!         The nonzeros in column j of L are in the
!         lcol_ptr(j), ... , lcol_ptr(j+1) - 1 positions of l.

!    lrow_ind is an integer array of dimension nnz + n*p.
!      On entry lrow_ind need not be specified.
!      On exit lrow_ind contains row indices for the strict lower
!         triangular part of L in compressed column storage.

!  Subprograms called

!    MINPACK-2  ......  dcauchy, dspcg, dssyax

!    Level 1 BLAS  ...  dcopy

!  MINPACK-2 Project. May 1999.
!  Argonne National Laboratory.
!  Chih-Jen Lin and Jorge J. More'.

!  **********

REAL (dp), PARAMETER :: p5=0.5_dp

!     Parameters for updating the iterates.

REAL (dp), PARAMETER :: eta0=1D-4, eta1=0.25_dp, eta2=0.75_dp

!     Parameters for updating the trust region size delta.

REAL (dp), PARAMETER :: sigma1=0.25_dp, sigma2=0.5_dp, sigma3=4.0_dp

LOGICAL   :: search
INTEGER   :: info, iter, iters
REAL (dp) :: alphac, fc, prered, actred, snorm, gs, alpha
CHARACTER (LEN=60) :: work

! REAL (dp) :: ddot, dnrm2
! EXTERNAL dcauchy, dspcg, dssyax
! EXTERNAL dcopy, ddot, dnrm2

!     Initialization section.

IF (task(1:5) == 'START') THEN
  
!        Initialize local variables.
  
  iter = 1
  iterscg = 0
  alphac = one
  work = 'COMPUTE'
  IF (ALLOCATED(xc)) DEALLOCATE( xc, s, indfree, dsave, wa, isave )
  ALLOCATE( xc(n), s(n), indfree(n), dsave(n), wa(n), isave(n) )
  
ELSE
  
!        Restore local variables.
  
  IF (isave(1) == 1) THEN
    work = 'COMPUTE'
  ELSE IF (isave(1) == 2) THEN
    work = 'EVALUATE'
  END IF
  iter = isave(2)
  iterscg = isave(3)
  fc = dsave(1)
  alphac = dsave(2)
END IF

!     Search for a lower function value.

search = .true.
DO WHILE (search)
  
!        Compute a step and evaluate the function at the trial point.
  
  IF (work == 'COMPUTE') THEN
    
!           Save the best function value and the best x.
    
    fc = f
    xc(1:n) = x(1:n)
    
!           Compute the Cauchy step and store in s.
    
    CALL dcauchy(n, x, xl, xu, a, adiag, acol_ptr, arow_ind, g, delta,  &
                 alphac, s)
    
!           Compute the projected Newton step.
    
    CALL dspcg(n, x, xl, xu, a, adiag, acol_ptr, arow_ind, g, delta,  &
               cgtol, s, 5, itermax, iters, info, b, bdiag, bcol_ptr, brow_ind,  &
               l, ldiag, lcol_ptr, lrow_ind)
    
    iterscg = iterscg + iters
    task = 'F'
  END IF
  
!        Evaluate the step and determine if the step is successful.
  
  IF (work == 'EVALUATE') THEN
    
!           Compute the predicted reduction.
    
    CALL dssyax(n, a, adiag, acol_ptr, arow_ind, s, wa)
    prered = DOT_PRODUCT( s(1:n), p5*wa(1:n) - g(1:n) )
    
!           Compute the actual reduction.
    
    actred =  fc - f
    
!           On the first iteration, adjust the initial step bound.
    
    snorm = dnrm2(n, s, 1)
    IF (iter == 1)  delta = MIN(delta, snorm)
    
!           Compute prediction alpha*snorm of the step.
    
    gs = DOT_PRODUCT( g(1:n), s(1:n) )
    IF (f-fc-gs <= zero) THEN
      alpha = sigma3
    ELSE
      alpha = MAX(sigma1, -p5*(gs/(f-fc-gs)))
    END IF
    
!           Update the trust region bound according to the ratio
!           of actual to predicted reduction.
    
    IF (actred < eta0*prered) THEN
      
!              Reduce delta.  Step is not successful.
      
      delta = MIN(MAX(alpha, sigma1)*snorm, sigma2*delta)
      
    ELSE IF (actred < eta1*prered) THEN
      
!              Reduce delta.  Step is not sufficiently successful.
      
      delta = MAX(sigma1*delta, MIN(alpha*snorm, sigma2*delta))
      
    ELSE IF (actred < eta2*prered) THEN
      
!              The ratio of actual to predicted reduction is in
!              the interval (eta1, eta2).  We are allowed to either
!              increase or decrease delta.
      
      delta = MAX(sigma1*delta, MIN(alpha*snorm, sigma3*delta))
      
    ELSE
      
!              The ratio of actual to predicted reduction exceeds eta2.
!              Do not decrease delta.
      
      delta = MAX(delta, MIN(alpha*snorm, sigma3*delta))
      
    END IF
    
!           Update the iterate.
    
    IF (actred > eta0*prered) THEN
      
!              Successful iterate.
      
      task = 'GH'
      iter = iter + 1
      
    ELSE