!        Compute q by solving L*q = A*z and save L*q for
!        use in updating the residual t.
  
  CALL dssyax(n, a, adiag, acol_ptr, arow_ind, z, q)
  z(1:n) = q(1:n)
  CALL dstrsol(n, l, ldiag, lcol_ptr, lrow_ind, q, 'N')
  
!        Compute alpha and determine sigma such that the trust region
!        constraint || w + sigma*p || = delta is satisfied.
  
  ptq = DOT_PRODUCT( p(1:n), q(1:n) )
  IF (ptq > zero) THEN
    alpha = rho/ptq
  ELSE
    alpha = zero
  END IF
  CALL dtrqsol(n, w, p, delta, sigma)
  
!        Exit if there is negative curvature or if the
!        iterates exit the trust region.
  
  IF (ptq <= zero .OR. alpha >= sigma) THEN
    w(1:n) = w(1:n) + sigma * p(1:n)
    IF (ptq <= zero) THEN
      info = 3
    ELSE
      info = 4
    END IF
    
    RETURN
    
  END IF
  
!        Update w and the residuals r and t.
!        Note that t = L*r.
  
  w(1:n) = w(1:n) + alpha * p(1:n)
  r(1:n) = r(1:n) - alpha * q(1:n)
  t(1:n) = t(1:n) - alpha * z(1:n)
  
!        Exit if the residual convergence test is satisfied.
  
  rtr = SUM( r(1:n)**2 )
  rnorm = SQRT(rtr)
  tnorm = dnrm2(n, t, 1)
  
  IF (tnorm <= tol) THEN
    info = 1
    RETURN
  END IF
  
  IF (rnorm <= stol) THEN
    info = 2
    RETURN
  END IF
  
!        Compute p = r + beta*p and update rho.
  
  beta = rtr/rho
  p(1:n) = beta * p(1:n)
  p(1:n) = p(1:n) + r(1:n)
  rho = rtr
  
END DO

info = 5

RETURN

END SUBROUTINE dtrpcg



SUBROUTINE dtrqsol(n, x, p, delta, sigma)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 11:18:32

INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(IN)   :: x(:)
REAL (dp), INTENT(IN)   :: p(:)
REAL (dp), INTENT(IN)   :: delta
REAL (dp), INTENT(OUT)  :: sigma

!  **********

!  Subroutine dtrqsol

!  This subroutine computes the largest (non-negative) solution
!  of the quadratic trust region equation

!        ||x + sigma*p|| = delta.

!  The code is only guaranteed to produce a non-negative solution
!  if ||x|| <= delta, and p != 0. If the trust region equation has
!  no solution, sigma = 0.

!  The subroutine statement ix

!    dtrqsol(n, x, p, delta, sigma)

!  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 must contain the vector x.
!      On exit x is unchanged.

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

!    delta is a REAL (dp) variable.
!      On entry delta specifies the scalar delta.
!      On exit delta is unchanged.

!    sigma is a REAL (dp) variable.
!      On entry sigma need not be specified.
!      On exit sigma contains the non-negative solution.

!  Subprograms called

!    Level 1 BLAS  ...  ddot

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

!  **********


REAL (dp) :: dsq, ptp, ptx, rad, xtx

! REAL (dp) :: ddot

ptx = DOT_PRODUCT( p(1:n), x(1:n) )
ptp = SUM( p(1:n)**2 )
xtx = SUM( x(1:n)**2 )
dsq = delta**2

!     Guard against abnormal cases.

rad = ptx**2 + ptp*(dsq - xtx)
rad = SQRT(MAX(rad, zero))

IF (ptx > zero) THEN
  sigma = (dsq - xtx)/(ptx + rad)
ELSE IF (rad > zero) THEN
  sigma = (rad - ptx)/ptp
ELSE
  sigma = zero
END IF

RETURN

END SUBROUTINE dtrqsol



SUBROUTINE dssyax(n, a, adiag, jptr, indr, x, y)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 14:03:11

INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(IN)   :: a(:)
REAL (dp), INTENT(IN)   :: adiag(:)
INTEGER, INTENT(IN)     :: jptr(:)    ! jptr(n+1)
INTEGER, INTENT(IN)     :: indr(:)
REAL (dp), INTENT(IN)   :: x(:)
REAL (dp), INTENT(OUT)  :: y(:)

!  **********

!  Subroutine dssyax

!  This subroutine computes the matrix-vector product y = A*x,
!  where A is a symmetric matrix with the strict lower triangular
!  part in compressed column storage.

!  The subroutine statement is

!    subroutine dssyax(n, a, adiag, jptr, indr, x, y)

!  where

!    n is an integer variable.
!      On entry n is the order of A.
!      On exit n is unchanged.

!    a is a REAL (dp) array of dimension *.
!      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.

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

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

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

!    y is a REAL (dp) array of dimension n.
!      On entry y need not be specified.
!      On exit y contains the product A*x.

!  MINPACK-2 Project. May 1998.
!  Argonne National Laboratory.

!  **********


INTEGER   :: i, j, k
REAL (dp) :: sum

y(1:n) = adiag(1:n)*x(1:n)

DO j = 1, n
  sum = zero
  DO i = jptr(j), jptr(j+1) - 1
    k = indr(i)
    sum = sum + a(i)*x(k)
    y(k) = y(k) + a(i)*x(j)
  END DO
  y(j) = y(j) + sum
END DO

RETURN

END SUBROUTINE dssyax



FUNCTION dnrm2 ( n, x, incx) RESULT(fn_val)

!  Euclidean norm of the n-vector stored in x() with storage increment incx .
!  if n <= 0 return with result = 0.
!  if n >= 1 then incx must be >= 1

!  c.l.lawson, 1978 jan 08
!  modified to correct failure to update ix, 1/25/92.
!  modified 3/93 to return if incx <= 0.
!  This version by Alan.Miller @ vic.cmis.csiro.au
!  Latest revision - 22 January 1999

!  four phase method using two built-in constants that are
!  hopefully applicable to all machines.
!      cutlo = maximum of  SQRT(u/eps)  over all known machines.
!      cuthi = minimum of  SQRT(v)      over all known machines.
!  where
!      eps = smallest no. such that eps + 1. > 1.
!      u   = smallest positive no.   (underflow limit)
!      v   = largest  no.            (overflow  limit)

!  brief outline of algorithm..

!  phase 1    scans zero components.
!  move to phase 2 when a component is nonzero and <= cutlo
!  move to phase 3 when a component is > cutlo
!  move to phase 4 when a component is >= cuthi/m
!  where m = n for x() real and m = 2*n for complex.

INTEGER, INTENT(IN)   :: n, incx
REAL (dp), INTENT(IN) :: x(:)
REAL (dp)             :: fn_val

! Local variables
INTEGER     :: i, ix, j, next
REAL (dp)   :: cuthi, cutlo, hitest, sum, xmax

IF(n <= 0 .OR. incx <= 0) THEN
  fn_val = zero
  RETURN
END IF

! Set machine-dependent constants

cutlo = SQRT( TINY(one) / EPSILON(one) )
cuthi = SQRT( HUGE(one) )

next = 1
sum = zero
i = 1
ix = 1
!                                                 begin main loop
20 SELECT CASE (next)
  CASE (1)
     IF( ABS(x(i)) > cutlo) GO TO 85
     next = 2
     xmax = zero
     GO TO 20

  CASE (2)
!                   phase 1.  sum is zero

     IF( x(i) == zero) GO TO 200
     IF( ABS(x(i)) > cutlo) GO TO 85

!                                prepare for phase 2.   x(i) is very small.
     next = 3
     GO TO 105

  CASE (3)
!                   phase 2.  sum is small.
!                             scale to avoid destructive underflow.

     IF( ABS(x(i)) > cutlo ) THEN
!                  prepare for phase 3.

       sum = (sum * xmax) * xmax
       GO TO 85
     END IF

  CASE (4)
     GO TO 110
END SELECT

! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                     common code for phases 2 and 4.
!                     in phase 4 sum is large.  scale to avoid overflow.

110 IF( ABS(x(i)) <= xmax ) GO TO 115
sum = one + sum * (xmax / x(i))**2
xmax = ABS(x(i))
GO TO 200

! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!                   phase 3.  sum is mid-range.  no scaling.

!     for real or d.p. set hitest = cuthi/n
!     for complex      set hitest = cuthi/(2*n)

85 hitest = cuthi / REAL( n, dp )

DO j = ix, n
  IF(ABS(x(i)) >= hitest) GO TO 100
  sum = sum + x(i)**2
  i = i + incx
END DO
fn_val = SQRT( sum )
RETURN

! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!                                prepare for phase 4.
!                                ABS(x(i)) is very large
100 ix = j
next = 4
sum = (sum / x(i)) / x(i)
!                                Set xmax; large if next = 4, small if next = 3
105 xmax = ABS(x(i))

115 sum = sum + (x(i)/xmax)**2

200 ix = ix + 1
i = i + incx
IF( ix <= n ) GO TO 20

!              end of main loop.

!              compute square root and adjust for scaling.

fn_val = xmax * SQRT(sum)

RETURN
END FUNCTION dnrm2



SUBROUTINE dicfs(n, nnz, a, adiag, acol_ptr, arow_ind,  &
                 l, ldiag, lcol_ptr, lrow_ind,  &
                 p, alpha)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 14:02:36

INTEGER, INTENT(IN)        :: n
INTEGER, INTENT(IN)        :: nnz
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(OUT)     :: l(:)          ! l(nnz+n*p)
REAL (dp), INTENT(OUT)     :: ldiag(:)
INTEGER, INTENT(OUT)       :: lcol_ptr(:)   ! lcol_ptr(n+1)
INTEGER, INTENT(OUT)       :: lrow_ind(:)   ! lrow_ind(nnz+n*p)
INTEGER, INTENT(IN)        :: p
REAL (dp), INTENT(IN OUT)  :: alpha

!  *********

!  Subroutine dicfs

!  Given a symmetric matrix A in compressed column storage, this
!  subroutine computes an incomplete Cholesky factor of A + alpha*D,
!  where alpha is a shift and D is the diagonal matrix with entries
!  set to the l2 norms of the columns of A.

!  The subroutine statement is

!    subroutine dicfs(n, nnz, a, adiag, acol_ptr, arow_ind,
!                     l, ldiag, lcol_ptr, lrow_ind,
!                     p, alpha)

!  where

!    n is an integer variable.
!      On entry n is the order of A.
!      On exit n is unchanged.

!    nnz is an integer variable.
!      On entry nnz is the number of nonzeros in the strict lower
!         triangular part of A.
!      On exit nnz is unchanged.

!    a is a REAL (dp) array of dimension nnz.
!      On entry a must