dtron.f90

牛顿优化算法源fortran代码

      
!              Unsuccessful iterate.
      
      task = 'NFGH'
      x(1:n) = xc(1:n)
      f = fc
      
    END IF
    
!           Test for convergence.
    
    IF (f < fmin) task = 'WARNING: F < FMIN'
    IF (ABS(actred) <= fatol .AND. prered <= fatol) task =  &
        'CONVERGENCE: FATOL TEST SATISFIED'
    IF (ABS(actred) <= frtol*ABS(f) .AND. prered <= frtol*ABS(f)) task =  &
        'CONVERGENCE: FRTOL TEST SATISFIED'
    
  END IF
  
!        Determine what needs to be done on the next call.
  
  IF (task == 'F') THEN
    work = 'EVALUATE'
  ELSE
    work = 'COMPUTE'
  END IF
  
!        Continue the search if the step is not successful.
  
  IF (task /= 'NFGH') search = .false.
  
END DO

!     Save local variables.

IF (work == 'COMPUTE') THEN
  isave(1) = 1
ELSE IF (work == 'EVALUATE') THEN
  isave(1) = 2
END IF

isave(2) = iter
isave(3) = iterscg
dsave(1) = fc
dsave(2) = alphac

RETURN
END SUBROUTINE dtron



SUBROUTINE dcauchy(n, x, xl, xu, a, diag, col_ptr, row_ind, g, delta,  &
                   alpha, s)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 11:18:33

INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN)      :: 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)      :: delta
REAL (dp), INTENT(IN OUT)  :: alpha
REAL (dp), INTENT(OUT)     :: s(:)


!  **********

!  Subroutine dcauchy

!  This subroutine computes a Cauchy step that satisfies a trust
!  region constraint and a sufficient decrease condition.

!  The Cauchy step is computed for the quadratic

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

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

!        s[alpha] = P[x - alpha*g] - x,

!  with P the projection onto the n-dimensional interval [xl,xu].
!  The Cauchy step satisfies the trust region constraint and the
!  sufficient decrease condition

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

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

!  The subroutine statement is

!    subroutine dcauchy(n, x, xl, xu, a, diag, col_ptr, row_ind, g, delta,
!                       alpha, s)

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

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

!    alpha is a REAL (dp) variable.
!      On entry alpha is the current estimate of the step.
!      On exit alpha defines the Cauchy step s[alpha].

!    s is a REAL (dp) array of dimension n.
!      On entry s need not be specified.
!      On exit s is the Cauchy step s[alpha].

!  Subprograms called

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

!    Level 1 BLAS  ...  ddot, dnrm2

!  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 and extrapolation factors.

REAL (dp), PARAMETER :: interpf=0.1_dp, extrapf=10.0_dp

LOGICAL   :: search, interp
INTEGER   :: nbrpt, nsteps
REAL (dp) :: alphas, brptmax, brptmin, gts, q
REAL (dp) :: wa(n)

! REAL (dp) :: dnrm2, ddot
! EXTERNAL dnrm2, ddot

!     Find the minimal and maximal break-point on x - alpha*g.

wa(1:n) = g(1:n)
wa(1:n) = - wa(1:n)
CALL dbreakpt(n, x, xl, xu, g, nbrpt, brptmin, brptmax)

!     Evaluate the initial alpha and decide if the algorithm
!     must interpolate or extrapolate.

CALL dgpstep(n, x, xl, xu, -alpha, g, s)
IF (dnrm2(n, s, 1) > delta) THEN
  interp = .true.
ELSE
  CALL dssyax(n, a, diag, col_ptr, row_ind, s, wa)
  gts = DOT_PRODUCT( g(1:n), s(1:n) )
  q = p5*DOT_PRODUCT( s(1:n), wa(1:n) ) + gts
  interp = (q >= mu0*gts)
END IF

nsteps = 1

!     Either interpolate or extrapolate to find a successful step.

IF (interp) THEN
  
!        Reduce alpha until a successful step is found.
  
  search = .true.
  DO WHILE (search)
    
!           This is a crude interpolation procedure that
!           will be replaced in future versions of the code.
    
    nsteps = nsteps + 1
    alpha = interpf*alpha
    CALL dgpstep(n, x, xl, xu, -alpha, g, s)
    IF (dnrm2(n, s, 1) <= delta) THEN
      CALL dssyax(n, a, diag, col_ptr, row_ind, s, wa)
      gts = DOT_PRODUCT( g(1:n), s(1:n) )
      q = p5*DOT_PRODUCT( s(1:n), wa(1:n) ) + gts
      search = (q > mu0*gts)
    END IF
  END DO
  
ELSE
  
!        Increase alpha until a successful step is found.
  
  search = .true.
  DO WHILE (search .AND. alpha <= brptmax)
    
!           This is a crude extrapolation procedure that
!           will be replaced in future versions of the code.
    
    nsteps = nsteps + 1
    alphas = alpha
    alpha = extrapf*alpha
    CALL dgpstep(n, x, xl, xu, -alpha, g, s)
    IF (dnrm2(n, s, 1) <= delta) THEN
      CALL dssyax(n, a, diag, col_ptr, row_ind, s, wa)
      gts = DOT_PRODUCT( g(1:n), s(1:n) )
      q = p5*DOT_PRODUCT( s(1:n), wa(1:n) ) + gts
      search = (q < mu0*gts)
    ELSE
      search = .false.
    END IF
  END DO
  
!        Recover the last successful step.
  
  alpha = alphas
  CALL dgpstep(n, x, xl, xu, -alpha, g, s)
END IF

RETURN

END SUBROUTINE dcauchy



SUBROUTINE dspcg(n, x, xl, xu, a, adiag, acol_ptr, arow_ind, g, delta,  &
                 rtol, s, nv, itermax, iters, info,  &
                 b, bdiag, bcol_ptr, brow_ind,  &
                 l, ldiag, lcol_ptr, lrow_ind)
 
! 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)      :: 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)      :: rtol
REAL (dp), INTENT(IN OUT)  :: s(:)
INTEGER, INTENT(IN)        :: nv
INTEGER, INTENT(IN)        :: itermax
INTEGER, INTENT(OUT)       :: iters
INTEGER, INTENT(OUT)       :: info
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(:)


!  *********

!  Subroutine dspcg

!  This subroutine generates a sequence of approximate minimizers
!  for the subproblem

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

!  The quadratic is defined by

!        q(x[0]+s) = 0.5*s'*A*s + g'*s,

!  where x[0] is a base point provided by the user, A is a symmetric
!  matrix in compressed column storage, and g is a vector.

!  At each stage we have an approximate minimizer x[k], and generate a direction
!  p[k] by using a preconditioned conjugate gradient method on the subproblem

!        min { q(x[k]+p) : || L'*p || <= delta, s(fixed) = 0 },

!  where fixed is the set of variables fixed at x[k], delta is the trust region
!  bound, and L is an incomplete Cholesky factorization of the submatrix

!        B = A(free:free),

!  where free is the set of free variables at x[k]. Given p[k],
!  the next minimizer x[k+1] is generated by a projected search.

!  The starting point for this subroutine is x[1] = x[0] + s, where
!  x[0] is a base point and s is the Cauchy step.

!  The subroutine converges when the step s satisfies

!        || (g + A*s)[free] || <= rtol*|| g[free] ||

!  In this case the final x is an approximate minimizer in the
!  face defined by the free variables.

!  The subroutine terminates when the trust region bound does
!  not allow further progress, that is, || L'*p[k] || = delta.
!  In this case the final x satisfies q(x) < q(x[k]).

!  The subroutine statement is

!    subroutine dspcg(n, x, xl, xu, a, adiag, acol_ptr, arow_ind, g, delta,
!                     rtol, s, nv, itermax, iters, info,
!                     b, bdiag, bcol_ptr, brow_ind,
!                     l, ldiag, lcol_ptr, lrow_ind)

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

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

!    rtol is a REAL (dp) variable.
!      On entry rtol specifies the accuracy of the final minimizer.
!      On exit rtol is unchanged.

!    s is a REAL (dp) array of dimension n.
!      On entry s is the Cauchy step.
!      On exit s contain the final step.

!    nv is an integer variable.
!      On entry nv specifies the amount of memory available for the
!         incomplete Cholesky factorization.
!      On exit nv 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.

!    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. The final step s satisfies
!                    || (g + A*s)[free] || <= rtol*|| g[free] ||,
!                    and the final x is an approximate minimizer
!                    in the face defined by the free variables.

!          info = 2  Termination. The trust region bound does
!                    not allow further progress.

!          info = 3  Failure to converge within itermax iterations.

!    b is a REAL (dp) array of dimension nnz.