dmain.f90

牛顿优化算法源fortran代码

!    s is a REAL (dp) array of dimension n.
!      On entry s specifies a vector s.
!      On exit s may be changed.

!    hs is a REAL (dp) array of dimension n.
!      On entry hs need not be specified.
!      On exit hs contains the product H*s where H is the Hessian matrix at x.

!    prob is a character*6 variable.
!      On entry prob specifies the problem.
!      On exit prob is set to 'ERROR' if prob is not an acceptable problem name.
!         Otherwise prob is unchanged.

!    par is a REAL (dp) variable.
!      On entry par specifies a problem-dependent parameter.
!      On exit par is unchanged.

!  Subprograms called

!    MINPACK-2 ... depths, dgl1hs, dgl2hs, dmsahs, dmsabc,
!                  dodchs, dpjbhs, dsschs

!  MINPACK-2 Project. March 1999.
!  Argonne National Laboratory.
!  Brett M. Averick and Jorge J. More'.

!  **********

REAL (dp), PARAMETER   :: b=10.0_dp
REAL (dp), ALLOCATABLE :: bottom(:), top(:), left(:), right(:)

! EXTERNAL depths, dgl1hs, dgl2hs, dmsahs, dmsabc, dodchs, dpjbhs, dsschs

!     Select a problem.

SELECT CASE ( prob(1:4) )
  CASE ( 'DEPT' )
    CALL depths(nx, ny, s, hs)
  CASE( 'DGL1' )
    CALL dgl1hs(n, x, s, hs, par)
  CASE( 'DGL2' )
    CALL dgl2hs(nx, ny, x, s, hs, INT(par))
  CASE( 'DMSA' )
    ALLOCATE( bottom(nx+2), top(nx+2), left(ny+2), right(ny+2) )
    CALL dmsabc(nx, ny, bottom, top, left, right)
    CALL dmsahs(nx, ny, x, s, hs, bottom, top, left, right)
    DEALLOCATE( bottom, top, left, right )
  CASE( 'DODC' )
    CALL dodchs(nx, ny, x, s, hs, par)
  CASE( 'DPJB' )
    CALL dpjbhs(nx, ny, s, hs, par, b)
  CASE( 'DSSC' )
    CALL dsschs(nx, ny, x, s, hs, par)
  CASE DEFAULT
    prob = 'ERROR'
END SELECT

RETURN
END SUBROUTINE dminhs



SUBROUTINE dminsp(n, nx, ny, nnz, indrow, indcol, prob)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 15:49:10

INTEGER, INTENT(IN)                :: n
INTEGER, INTENT(IN)                :: nx
INTEGER, INTENT(IN)                :: ny
INTEGER, INTENT(OUT)               :: nnz
INTEGER, INTENT(OUT)               :: indrow(:)
INTEGER, INTENT(OUT)               :: indcol(:)
CHARACTER (LEN=*), INTENT(IN OUT)  :: prob

!  *********

!  Subroutine dminsp

!  This subroutine computes the sparsity pattern for
!  the minimization problem from the MINPACK-2 test problem
!  collection specified by the character variable prob.

!  The subroutine statement is

!    dminsp(n, nx, ny, nnz, indrow, indcol, prob)

!  where

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

!    nx is an integer variable.
!      On entry nx is the number of grid points in the first
!         coordinate direction.
!      On exit nx is unchanged.

!    ny is an integer variable.
!      On entry ny is the number of grid points in the second
!         coordinate direction. If the problem is formulated in
!         one spatial dimension, ny = 1.
!      On exit ny is unchanged.

!    nnz is an integer variable.
!      On entry nnz need not be specified.
!      On exit nnz is set to the number of nonzeros in the
!         lower triangle of the Hessian matrix.

!    indrow is an integer array of dimension at least nnz.
!      On entry indrow need not be specified.
!      On exit indrow contains the row indices of the nonzeros
!        in the lower triangle of the Hessian matrix.

!    indcol is an integer array of dimension at least nnz.
!      On entry indcol need not be specified.
!      On exit indcol contains the column indices of the nonzeros
!         in the lower triangle of the Hessian matrix.

!    prob is a character*6 variable.
!      On entry prob specifies the problem.
!      On exit prob is set to 'ERROR' if prob is not an
!         acceptable problem name. Otherwise prob is unchanged.

!  Subprograms called

!    MINPACK-2 ... deptsp, dgl1sp, dgl2sp, dmsasp, dmsabc,
!                  dodcsp, dpjbsp, dsscsp

!  MINPACK-2 Project. March 1999.
!  Argonne National Laboratory.
!  Brett M. Averick and Jorge J. More'.

!  **********

! EXTERNAL deptsp, dgl1sp, dgl2sp, dmsasp, dodcsp, dpjbsp, dsscsp

!     Select a problem.

SELECT CASE ( prob(1:4) )
  CASE ( 'DEPT' )
    CALL deptsp(nx, ny, nnz, indrow, indcol)
  CASE( 'DGL1' )
    CALL dgl1sp(n, nnz, indrow, indcol)
  CASE( 'DGL2' )
    CALL dgl2sp(nx, ny, nnz, indrow, indcol)
  CASE( 'DMSA' )
    CALL dmsasp(nx, ny, nnz, indrow, indcol)
  CASE( 'DODC' )
    CALL dodcsp(nx, ny, nnz, indrow, indcol)
  CASE( 'DPJB' )
    CALL dpjbsp(nx, ny, nnz, indrow, indcol)
  CASE( 'DSSC' )
    CALL dsscsp(nx, ny, nnz, indrow, indcol)
  CASE DEFAULT
    prob = 'ERROR'
END SELECT

RETURN
END SUBROUTINE dminsp



SUBROUTINE dminxb(n, nx, ny, xl, xu, prob)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 15:49:11

INTEGER, INTENT(IN)            :: n
INTEGER, INTENT(IN)            :: nx
INTEGER, INTENT(IN)            :: ny
REAL (dp), INTENT(OUT)         :: xl(:)
REAL (dp), INTENT(OUT)         :: xu(:)
CHARACTER (LEN=*), INTENT(IN)  :: prob

!  *********

!  Subroutine dminxb

!  This subroutine computes the lower and upper bounds for
!  the minimization problem from the MINPACK-2 test problem
!  collection specified by the character variable prob.

!  The subroutine statement is

!    dminxb(n, nx, ny, xl, xu, prob)

!  where

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

!    nx is an integer variable.
!      On entry nx is the number of grid points in the first
!         coordinate direction.
!      On exit nx is unchanged.

!    ny is an integer variable.
!      On entry ny is the number of grid points in the second
!         coordinate direction. If the problem is formulated in
!         one spatial dimension, ny = 1.
!      On exit ny is unchanged.

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

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

!    prob is a character*6 variable.
!      On entry prob specifies the problem.
!      On exit prob is unchanged.

!  MINPACK-2 Project. March 1999.
!  Argonne National Laboratory.
!  Brett M. Averick and Jorge J. More'.

!  **********

REAL (dp), PARAMETER :: p0001=1.0D-4, p001=1.0D-3, p01=1.0D-2, p1=1.0D-1

REAL (dp), PARAMETER :: oned2=1.0D2, xbmax=1.0D20

INTEGER   :: i, j, k
REAL (dp) :: hx, hy
REAL (dp) :: tmp

!     Select a problem.

IF (prob(1:4) == 'DEPT') THEN
  hx = one/DBLE(nx+1)
  hy = one/DBLE(ny+1)
  DO j = 1, ny
    tmp = DBLE(MIN(j,ny-j+1))*hy
    DO i = 1, nx
      k = nx*(j-1) + i
      xu(k) = MIN(DBLE(MIN(i,nx-i+1))*hx,tmp)
      xl(k) = -xu(k)
    END DO
  END DO
ELSE IF (prob(1:5) == 'DMSA1') THEN
  DO i = 1, n
    xl(i) = -4*p1
    xu(i) =  4*p1
  END DO
ELSE IF (prob(1:5) == 'DMSA2') THEN
  DO i = 1, n
    xl(i) = -2*p1
    xu(i) =  2*p1
  END DO
ELSE IF (prob(1:5) == 'DMSA3') THEN
  DO i = 1, n
    xl(i) = -p1
    xu(i) =  p1
  END DO
ELSE IF (prob(1:4) == 'DPJB') THEN
  DO i = 1, n
    xl(i) = zero
    xu(i) = oned2
  END DO
ELSE IF (prob(1:5) == 'DSSC1') THEN
  DO i = 1, n
    xl(i) = p1
    xu(i) = one
  END DO
ELSE IF (prob(1:5) == 'DSSC2') THEN
  DO i = 1, n
    xl(i) = p01
    xu(i) = one
  END DO
ELSE IF (prob(1:5) == 'DSSC3') THEN
  DO i = 1, n
    xl(i) = p001
    xu(i) = one
  END DO
ELSE IF (prob(1:5) == 'DSSC4') THEN
  DO i = 1, n
    xl(i) = p0001
    xu(i) = one
  END DO
ELSE
  DO i = 1, n
    xl(i) = -xbmax
    xu(i) =  xbmax
  END DO
END IF

RETURN
END SUBROUTINE dminxb



FUNCTION dgpnrm2(n, x, xl, xu, g) RESULT(norm)
 
! 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)  :: g(:)
REAL (dp)              :: norm

!  **********

!  Function dgpnrm2

!  This function computes the Euclidean norm of the projected gradient at x.

!  The function statement is

!    function dgpnrm2(n, x, xl, xu, g)

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

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

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

!  **********

INTEGER :: i

norm = zero
DO i = 1, n
  IF (xl(i) /= xu(i)) THEN
    IF (x(i) == xl(i)) THEN
      norm = norm + MIN(g(i), zero)**2
    ELSE IF (x(i) == xu(i)) THEN
      norm = norm + MAX(g(i), zero)**2
    ELSE
      norm = norm + g(i)**2
    END IF
  END IF
END DO
norm = SQRT(norm)

RETURN

END FUNCTION dgpnrm2



SUBROUTINE dsphesd(n, row_ind, col_ind, row_ptr, col_ptr, listp, ngrp,  &
                   maxgrp, numgrp, eta, fhesd, fhes, diag)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 11:18:32

INTEGER, INTENT(IN)      :: n
INTEGER, INTENT(IN OUT)  :: row_ind(:)
INTEGER, INTENT(IN OUT)  :: col_ind(:)
INTEGER, INTENT(IN)      :: row_ptr(:)    ! row_ptr(n+1)
INTEGER, INTENT(IN)      :: col_ptr(:)    ! col_ptr(n+1)
INTEGER, INTENT(OUT)     :: listp(:)
INTEGER, INTENT(IN)      :: ngrp(:)
INTEGER, INTENT(IN)      :: maxgrp
INTEGER, INTENT(IN)      :: numgrp
REAL (dp), INTENT(IN)    :: eta(:)
REAL (dp), INTENT(IN)    :: fhesd(:)
REAL (dp), INTENT(OUT)   :: fhes(:)
REAL (dp), INTENT(OUT)   :: diag(:)

!  **********

!  Subroutine dsphesd

!  This subroutine computes an approximation to the (symmetric) Hessian matrix
!  of a function by a substitution method.
!  The lower triangular part of the approximation is stored in compressed
!  column storage.

!  This subroutine requires a symmetric permutation of the Hessian matrix and
!  a partition of the columns of the Hessian matrix consistent with the
!  determination of the Hessian matrix by a lower triangular substitution
!  method.   This information can be provided by subroutine dssm.

!  The symmetric permutation of the Hessian matrix is defined by the array
!  listp.  This array is only used internally.

!  The partition of the Hessian matrix is defined by the array ngrp by setting
!  ngrp(j) to the group number of column j.
!  The user must provide an approximation to the columns of the Hessian matrix
!  in each group by specifying a difference parameter vector eta and an
!  approximation to A*d where A is the Hessian matrix and the vector d is
!  defined by the following section of code.

!        do j = 1, n
!           d(j) = 0.0
!           if (ngrp(j) .eq. numgrp) d(j) = eta(j)
!        end do

!  In the above code numgrp is a group number and eta(j) is the
!  difference parameter used to approximate column j of the
!  Hessian matrix. Suitable values for eta(j) must be provided.

!  As mentioned above, an approximation to A*d must be provided.
!  For example, if grad f(x) is the gradient of the function at x, then

!        grad f(x+d) - grad f(x)

!  corresponds to the forward difference approximation.

!  The lower triangular substitution method requires that the approximations
!  to A*d for all the groups be stored in special locations of the array fhes.
!  This is done by calls with numgrp = 1, 2, ... ,maxgrp. On the call with
!  numgrp = maxgrp, the array fhes is overwritten with the approximation to the
!  lower triangular part of the Hessian matrix.

!  The subroutine statement is

!    subroutine dsphesd(n, row_ind, col_ind, row_ptr, col_ptr, listp, ngrp,
!                       maxgrp, numgrp, eta, fhesd, fhes, diag)

!  where

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

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