!    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.

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

!  **********

INTEGER   :: i

DO i = 1, n
  x(i) = MAX(xl(i), MIN(x(i), xu(i)))
END DO

RETURN
END SUBROUTINE dmid



SUBROUTINE dprsrch(n, x, xl, xu, a, diag, col_ptr, row_ind, g, w)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 11:18:32

INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: x(:)
REAL (dp), INTENT(IN)      :: xl(:)
REAL (dp), INTENT(IN)      :: xu(:)
REAL (dp), INTENT(IN)      :: a(:)
REAL (dp), INTENT(IN)      :: diag(:)
INTEGER, INTENT(IN)        :: col_ptr(:)    ! col_ptr(n+1)
INTEGER, INTENT(IN)        :: row_ind(:)
REAL (dp), INTENT(IN)      :: g(:)
REAL (dp), INTENT(IN OUT)  :: w(:)

!  **********

!  Subroutine dprsrch

!  This subroutine uses a projected search to compute a step
!  that satisfies a sufficient decrease condition for the quadratic

!        q(s) = 0.5*s'*A*s + g'*s,

!  where A is a symmetric matrix in compressed column storage,
!  and g is a vector. Given the parameter alpha, the step is

!        s[alpha] = P[x + alpha*w] - x,

!  where w is the search direction and P the projection onto the
!  n-dimensional interval [xl,xu]. The final step s = s[alpha]
!  satisfies the sufficient decrease condition

!        q(s) <= mu_0*(g'*s),

!  where mu_0 is a constant in (0,1).

!  The search direction w must be a descent direction for the
!  quadratic q at x such that the quadratic is decreasing
!  in the ray  x + alpha*w for 0 <= alpha <= 1.

!  The subroutine statement is

!    subroutine dprsrch(n, x, xl, xu, a, diag, col_ptr, row_ind, g, w)

!  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 set to the final point P[x + alpha*w].

!    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.

!    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.

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

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

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

!    g is a REAL (dp) array of dimension n.
!      On entry g specifies the vector g.
!      On exit g is unchanged.

!    w is a double prevision array of dimension n.
!      On entry w specifies the search direction.
!      On exit w is the step s[alpha].

!  Subprograms called

!    MINPACK-2  ......  dbreakpt, dgpstep, dmid, dssyax

!    Level 1 BLAS  ...  daxpy, dcopy, ddot

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

!  **********

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

!     Constant that defines sufficient decrease.

REAL (dp), PARAMETER :: mu0=0.01_dp

!     Interpolation factor.

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

REAL (dp) :: wa1(n), wa2(n)
LOGICAL   :: search
INTEGER   :: nbrpt, nsteps
REAL (dp) :: alpha, brptmin, brptmax, gts, q

! REAL (dp) :: ddot
! EXTERNAL daxpy, dcopy, ddot
! EXTERNAL dbreakpt, dgpstep, dmid, dssyax

!     Set the initial alpha = 1 because the quadratic function is
!     decreasing in the ray x + alpha*w for 0 <= alpha <= 1.

alpha = one
nsteps = 0

!     Find the smallest break-point on the ray x + alpha*w.

CALL dbreakpt(n, x, xl, xu, w, nbrpt, brptmin, brptmax)

!     Reduce alpha until the sufficient decrease condition is
!     satisfied or x + alpha*w is feasible.

search = .true.
DO WHILE (search .AND. alpha > brptmin)
  
!        Calculate P[x + alpha*w] - x and check the sufficient
!        decrease condition.
  
  nsteps = nsteps + 1
  CALL dgpstep(n, x, xl, xu, alpha, w, wa1)
  CALL dssyax(n, a, diag, col_ptr, row_ind, wa1, wa2)
  gts = DOT_PRODUCT( g(1:n), wa1(1:n) )
  q = p5*DOT_PRODUCT( wa1(1:n), wa2(1:n) ) + gts
  IF (q <= mu0*gts) THEN
    search = .false.
  ELSE
    
!           This is a crude interpolation procedure that
!           will be replaced in future versions of the code.
    
    alpha = interpf*alpha
    
  END IF
END DO

!     Force at least one more constraint to be added to the active
!     set if alpha < brptmin and the full step is not successful.
!     There is sufficient decrease because the quadratic function
!     is decreasing in the ray x + alpha*w for 0 <= alpha <= 1.

IF (alpha < one .AND. alpha < brptmin) alpha = brptmin

!     Compute the final iterate and step.

CALL dgpstep(n, x, xl, xu, alpha, w, wa1)
x(1:n) = x(1:n) + alpha * w(1:n)
CALL dmid(n, x, xl, xu)
w(1:n) = wa1(1:n)

RETURN

END SUBROUTINE dprsrch



SUBROUTINE dtrpcg(n, a, adiag, acol_ptr, arow_ind, g, delta,  &
                  l, ldiag, lcol_ptr, lrow_ind,   &
                  tol, stol, itermax, w, iters, info)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 11:18:33

INTEGER, INTENT(IN)     :: n
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)   :: g(:)
REAL (dp), INTENT(IN)   :: delta
REAL (dp), INTENT(IN)   :: l(:)
REAL (dp), INTENT(IN)   :: ldiag(:)
INTEGER, INTENT(IN)     :: lcol_ptr(:)    ! lcol_ptr(n+1)
INTEGER, INTENT(IN)     :: lrow_ind(:)
REAL (dp), INTENT(IN)   :: tol
REAL (dp), INTENT(IN)   :: stol
INTEGER, INTENT(IN)     :: itermax
REAL (dp), INTENT(OUT)  :: w(:)
INTEGER, INTENT(OUT)    :: iters
INTEGER, INTENT(OUT)    :: info

!  *********

!  Subroutine dtrpcg

!  Given a sparse symmetric matrix A in compressed column storage,
!  this subroutine uses a preconditioned conjugate gradient method
!  to find an approximate minimizer of the trust region subproblem

!        min { q(s) : || L'*s || <= delta }.

!  where q is the quadratic

!        q(s) = 0.5*s'*A*s + g'*s,

!  A is a symmetric matrix in compressed column storage, L is a lower
!  triangular matrix in compressed column storage, and g is a vector.

!  This subroutine generates the conjugate gradient iterates for
!  the equivalent problem

!        min { Q(w) : || w || <= delta }.

!  where Q is the quadratic defined by

!        Q(w) = q(s),      w = L'*s.

!  Termination occurs if the conjugate gradient iterates leave
!  the trust region, a negative curvature direction is generated,
!  or one of the following two convergence tests is satisfied.

!  Convergence in the original variables:

!        || grad q(s) || <= tol

!  Convergence in the scaled variables:

!        || grad Q(w) || <= stol

!  Note that if w = L'*s, then L*grad Q(w) = grad q(s).

!  The subroutine statement is

!    subroutine dtrcg(n, a, adiag, acol_ptr, arow_ind, g, delta,
!                    l, ldiag, lcol_ptr, lrow_ind,
!                    tol, stol, itermax, w, iters, info)

!  where

!    n is an integer variable.
!      On entry n is the number of variables.
!      On exit n 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.

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

!    delta is a REAL (dp) variable.
!      On entry delta is the trust region size.
!      On exit delta is unchanged.

! The descriptions below for l, ldiag, lcol_ptr & lrow_ind appear to be wrong.
! All of these quantities must be INPUT; they are not changed.

!    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.

!    tol is a REAL (dp) variable.
!      On entry tol specifies the convergence test
!         in the un-scaled variables.
!      On exit tol is unchanged

!    stol is a REAL (dp) variable.
!      On entry stol specifies the convergence test
!         in the scaled variables.
!      On exit stol 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.

!    w is a REAL (dp) array of dimension n.
!      On entry w need not be specified.
!      On exit w contains the final conjugate gradient iterate.

!    iters is an integer variable.
!      On entry iters need not be specified.
!      On exit iters is set to the number of conjugate
!         gradient iterations.

!    info is an integer variable.
!      On entry info need not be specified.
!      On exit info is set as follows:

!          info = 1  Convergence in the original variables.
!                    || grad q(s) || <= tol

!          info = 2  Convergence in the scaled variables.
!                    || grad Q(w) || <= stol

!          info = 3  Negative curvature direction generated.
!                    In this case || w || = delta and a direction
!                    of negative curvature w can be recovered by
!                    solving L'*w = p.

!          info = 4  Conjugate gradient iterates exit the
!                    trust region. In this case || w || = delta.

!          info = 5  Failure to converge within itermax iterations.

!  Subprograms called

!    MINPACK-2  ......  dtrqsol, dstrsol, dssyax

!    Level 1 BLAS  ...  daxpy, dcopy, ddot

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

!  **********

REAL (dp) :: p(n), q(n), r(n), t(n), z(n)
REAL (dp) :: alpha, beta, ptq, rho, rtr, sigma
REAL (dp) :: rnorm, rnorm0, tnorm

! REAL (dp) :: ddot
! EXTERNAL dtrqsol, dstrsol, dssyax
! EXTERNAL daxpy, dcopy, ddot

!     Initialize the iterate w and the residual r.

w(1:n) = zero

!     Initialize the residual t of grad q to -g.
!     Initialize the residual r of grad Q by solving L*r = -g.
!     Note that t = L*r.

t(1:n) = g(1:n)
t(1:n) = - t(1:n)
r(1:n) = t(1:n)
CALL dstrsol(n, l, ldiag, lcol_ptr, lrow_ind, r, 'N')

!     Initialize the direction p.

p(1:n) = r(1:n)

!     Initialize rho and the norms of r and t.

rho = SUM( r(1:n)**2 )
rnorm0 = SQRT(rho)

!     Exit if g = 0.

IF (rnorm0 == zero) THEN
  info = 1
  RETURN
END IF

iters = 0
DO iters = 1, itermax
  
!        Compute z by solving L'*z = p.
  
  z(1:n) = p(1:n)
  CALL dstrsol(n, l, ldiag, lcol_ptr, lrow_ind, z, 'T')