dtron.f90

牛顿优化算法源fortran代码

!      On entry b need not be specified.
!      On exit b contains the strict lower triangular part
!         of B in compressed column storage.

!    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*nv.
!      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*nv.
!      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  ......  dmid, dprsrch, dtrpcg,
!                       dssyax, dstrsol

!    Level 1 BLAS  ...  daxpy, dnrm2

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

!  **********

INTEGER   :: indfree(n), iwa(3*n)
REAL (dp) :: gfree(n), w(n), wa(3*n)
INTEGER   :: infotr, ip, j, jfree, itertr, nfaces, nfree, nnz
REAL (dp) :: alpha, gfnorm, gfnormf
REAL (dp) :: tol, stol

! REAL (dp) :: dnrm2
! EXTERNAL dmid, dprsrch, dtrpcg, dssyax, dstrsol
! EXTERNAL daxpy, dnrm2

!     Compute A*(x[1] - x[0]) and store in w.

CALL dssyax(n, a, adiag, acol_ptr, arow_ind, s, w)

!     Compute the Cauchy point.

x(1:n) = x(1:n) + s(1:n)
CALL dmid(n, x, xl, xu)

!     Start the main iteration loop.
!     There are at most n iterations because at each iteration
!     at least one variable becomes active.

iters = 0
DO nfaces = 1, n
  
!        Determine the free variables at the current minimizer.
!        The indices of the free variables are stored in the first
!        n free positions of the array indfree.
!        The array iwa is used to detect free variables by setting
!        iwa(i) = 1 if the ith variable is free, otherwise iwa(i) = 0.
  
  nfree  = 0
  DO j = 1, n
    IF (xl(j) < x(j) .AND. x(j) < xu(j)) THEN
      nfree = nfree + 1
      indfree(nfree) = j
      iwa(j) = nfree
    ELSE
      iwa(j) = 0
    END IF
  END DO
  
!        Exit if there are no free constraints.
  
  IF (nfree == 0) THEN
    info = 1
    RETURN
  END IF
  
!        Obtain the submatrix of A for the free variables.
!        Recall that iwa allows the detection of free variables.
  
  bcol_ptr(1) = 1
  nnz = 0
  DO j = 1, nfree
    jfree = indfree(j)
    bdiag(j) = adiag(jfree)
    DO ip = acol_ptr(jfree), acol_ptr(jfree+1)-1
      IF (iwa(arow_ind(ip)) > 0) THEN
        nnz = nnz + 1
        brow_ind(nnz) = iwa(arow_ind(ip))
        b(nnz) = a(ip)
      END IF
    END DO
    bcol_ptr(j+1) = nnz + 1
  END DO
  
!        Compute the incomplete Cholesky factorization.
!        The variable nnz is used to specify an upper bound on the
!        length of b, so we need to make sure that nnz >= 1.
  
  alpha = zero
  nnz = MAX(nnz,1)
  CALL dicfs(nfree, nnz, b, bdiag, bcol_ptr, brow_ind,  &
             l, ldiag, lcol_ptr, lrow_ind, nv, alpha)
  
!        Compute the gradient grad q(x[k]) = g + A*(x[k] - x[0]),
!        of q at x[k] for the free variables.
!        Recall that w contains  A*(x[k] - x[0]).
!        Compute the norm of the reduced gradient Z'*g.
  
  DO j = 1, nfree
    gfree(j) = w(indfree(j)) + g(indfree(j))
    wa(j) = g(indfree(j))
  END DO
  gfnorm = dnrm2(nfree, wa, 1)
  
!        Solve the trust region subproblem in the free variables
!        to generate a direction p[k]. Store p[k] in the array w.
  
  tol = rtol*gfnorm
  stol = zero
  
  CALL dtrpcg(nfree, b, bdiag, bcol_ptr, brow_ind, gfree, delta,  &
              l, ldiag, lcol_ptr, lrow_ind, tol, stol, itermax, w, itertr,  &
              infotr)
  
  iters = iters + itertr
  CALL dstrsol(nfree, l, ldiag, lcol_ptr, lrow_ind, w, 'T')
  
!        Use a projected search to obtain the next iterate.
!        The projected search algorithm stores s[k] in w.
  
  DO j = 1, nfree
    wa(j) = x(indfree(j))
    wa(n+j) = xl(indfree(j))
    wa(2*n+j) = xu(indfree(j))
  END DO
  
  CALL dprsrch(nfree, wa, wa(n+1:), wa(2*n+1:),  &
               b, bdiag, bcol_ptr, brow_ind, gfree, w)
  
!        Update the minimizer and the step.
!        Note that s now contains x[k+1] - x[0].
  
  DO j = 1, nfree
    x(indfree(j)) = wa(j)
    s(indfree(j)) = s(indfree(j)) + w(j)
  END DO
  
!        Compute A*(x[k+1] - x[0]) and store in w.
  
  CALL dssyax(n, a, adiag, acol_ptr, arow_ind, s, w)
  
!        Compute the gradient grad q(x[k+1]) = g + A*(x[k+1] - x[0])
!        of q at x[k+1] for the free variables.
  
  DO j = 1, nfree
    gfree(j) = w(indfree(j)) + g(indfree(j))
  END DO
  gfnormf = dnrm2(nfree, gfree, 1)
  
!        Convergence and termination test.
!        We terminate if the preconditioned conjugate gradient method
!        encounters a direction of negative curvature, or
!        if the step is at the trust region bound.
  
  IF (gfnormf <= rtol*gfnorm) THEN
    info = 1
    RETURN
  ELSE IF (infotr == 3 .OR. infotr == 4) THEN
    info = 2
    RETURN
  ELSE IF (iters > itermax) THEN
    info = 3
    RETURN
  END IF
  
END DO

RETURN

END SUBROUTINE dspcg



SUBROUTINE dbreakpt(n, x, xl, xu, w, nbrpt, brptmin, brptmax)
 
! 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)   :: xl(:)
REAL (dp), INTENT(IN)   :: xu(:)
REAL (dp), INTENT(IN)   :: w(:)
INTEGER, INTENT(OUT)    :: nbrpt
REAL (dp), INTENT(OUT)  :: brptmin
REAL (dp), INTENT(OUT)  :: brptmax


!  **********

!  Subroutine dbreakpt

!  This subroutine computes the number of break-points, and
!  the minimal and maximal break-points of the projection of
!  x + alpha*w on the n-dimensional interval [xl,xu].

!  The subroutine statement is

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

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

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

!    nbrpt is an integer variable.
!      On entry nbrpt need not be specified.
!      On exit nbrpt is the number of break points.

!    brptmin is a REAL (dp) variable
!      On entry brptmin need not be specified.
!      On exit brptmin is minimal break-point.

!    brptmax is a REAL (dp) variable
!      On entry brptmax need not be specified.
!      On exit brptmax is maximal break-point.

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

!  **********


INTEGER   :: i
REAL (dp) :: brpt

nbrpt = 0
DO i = 1, n
  IF (x(i) < xu(i) .AND. w(i) > zero) THEN
    nbrpt = nbrpt + 1
    brpt =  (xu(i) - x(i))/w(i)
    IF (nbrpt == 1) THEN
      brptmin = brpt
      brptmax = brpt
    ELSE
      brptmin = MIN(brpt,brptmin)
      brptmax = MAX(brpt,brptmax)
    END IF
  ELSE IF (x(i) > xl(i) .AND. w(i) < zero) THEN
    nbrpt = nbrpt + 1
    brpt = (xl(i) - x(i))/w(i)
    IF (nbrpt == 1) THEN
      brptmin = brpt
      brptmax = brpt
    ELSE
      brptmin = MIN(brpt,brptmin)
      brptmax = MAX(brpt,brptmax)
    END IF
  END IF
END DO

!     Handle the exceptional case.

IF (nbrpt == 0) THEN
  brptmin = zero
  brptmax = zero
END IF

RETURN

END SUBROUTINE dbreakpt



SUBROUTINE dgpstep(n, x, xl, xu, alpha, w, s)
 
! 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)   :: xl(:)
REAL (dp), INTENT(IN)   :: xu(:)
REAL (dp), INTENT(IN)   :: alpha
REAL (dp), INTENT(IN)   :: w(:)
REAL (dp), INTENT(OUT)  :: s(:)

!  **********

!  Subroutine dgpstep

!  This subroutine computes the gradient projection step

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

!  where P is the projection on the n-dimensional interval [xl,xu].

!  The subroutine statement is

!    subroutine  dgpstep(n,x,xl,xu,alpha,w,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.

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

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

!    s is a REAL (dp) array of dimension n.
!      On entry s need not be specified.
!      On exit s contains the gradient projection step.

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

!  **********

INTEGER   :: i

!     This computation of the gradient projection step avoids
!     rounding errors for the components that are feasible.

DO i = 1, n
  IF (x(i) + alpha*w(i) < xl(i)) THEN
    s(i) = xl(i) - x(i)
  ELSE IF (x(i) + alpha*w(i) > xu(i)) THEN
    s(i) = xu(i) - x(i)
  ELSE
    s(i) = alpha*w(i)
  END IF
END DO

RETURN

END SUBROUTINE dgpstep



SUBROUTINE dmid(n, x, xl, xu)
 
! 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(:)

!  **********

!  Subroutine dmid

!  This subroutine computes the projection of x
!  on the n-dimensional interval [xl,xu].

!  The subroutine statement is

!    subroutine dmid(n, x, xl, xu)

!  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 projection of x on [xl,xu].