dmain.f90

牛顿优化算法源fortran代码

PROGRAM dmain
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 15:45:25
! Latest version - 5 July 1999

USE tron
USE coloring
IMPLICIT NONE

INTEGER, PARAMETER  :: nmax=90000, nnzmax=10*nmax
INTEGER, PARAMETER  :: nread=1, nwrite=2
!     **********

!     Driver for bound-constrained problems.

!     Subprograms called

!       USER ........... dminfg, dminhs, dminsp, dminxb

!       MINPACK-2 ...... dtron, dgpnrm2, dsphesd, dmid, dtimer

!       Level 1 BLAS ... dnrm2

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

!     **********
REAL (dp), PARAMETER :: zero=0.0_dp, one=1.0_dp

CHARACTER (LEN=60) :: task

LOGICAL   :: search
INTEGER   :: arow_ind(nnzmax), acol_ind(nnzmax), brow_ind(nnzmax)
INTEGER   :: arow_ptr(nmax+1), acol_ptr(nmax+1), bcol_ptr(nmax+1),  &
             lrow_ind(nnzmax), lcol_ptr(nmax+1)
INTEGER   :: n, itermax
REAL (dp) :: f, delta
REAL (dp) :: x(nmax), xl(nmax), xu(nmax), g(nmax), hs(nmax)
REAL (dp) :: adiag(nmax), bdiag(nmax), ldiag(nmax)
REAL (dp) :: a(nnzmax), b(nnzmax), l(nnzmax)

!     Tolerances.

REAL (dp) :: cgtol, frtol, fatol, fmin

!     Summary information.

INTEGER   :: iterscg, nbind, nfev, nfree, ngev, nhev, nhsev

!     Evaluation of the Hessian matrix.

INTEGER   :: maxgrp, numgrp
INTEGER   :: listp(nmax), ngrp(nmax)
REAL (dp) :: y(nmax), eta(nmax)

!     Test problems.

CHARACTER (LEN=6) :: prob
CHARACTER (LEN=2) :: ch
INTEGER           :: i, info, j, maxfev, nnz, nx, ny
REAL (dp)         :: gnorm, gnorm0, par, gtol

!     Timing.

REAL :: ttimes, ttimef, times, timef
REAL :: fgtime, htime, ttime

OPEN (nread, FILE='tron.dat', STATUS='OLD')
OPEN (nwrite, FILE='tron.inf')

DO
  
  READ (nread,*) prob, n, nx, ny, par
  WRITE (*,*)    prob, n, nx, ny, par
  
  IF (prob(1:4) == 'STOP') THEN
    CLOSE(nread)
    STOP
  END IF
  
!        Generate the initial point and project into [xl,xu].
  
  ch = 'XS'
  CALL dminfg(n, nx, ny, x, f, g, ch, prob, par)
  CALL dminxb(n, nx, ny, xl, xu, prob)
  CALL dmid(n, x, xl, xu)
  
!        Initialize variables.
  
  nfev = 0
  ngev = 0
  nhev = 0
  nhsev = 0
  fgtime = zero
  htime = zero
  
!        Set parameters.
  
  itermax = n
  maxfev = 1000
  fatol = zero
  frtol = 1.d-12
  fmin = -1.0D+32
  cgtol = 0.1_dp
  gtol = 1.0D-5
  
  CALL dtimer(ttimes)
  
!        Calculate the sparsity pattern.
  
  CALL dminsp(n, nx, ny, nnz, arow_ind, acol_ind, prob)
  CALL dsetsp(n, nnz, arow_ind, acol_ind, acol_ptr, arow_ptr, 1,  &
              info, listp, ngrp, maxgrp)
  IF (info <= 0) THEN
    WRITE (nwrite,*) 'ERROR: INFO IN SUBROUTINE SETSP IS ', info
    STOP
  END IF
  
!        Start the iteration.
  
  task = 'START'
  search = .TRUE.
  DO WHILE (search)
    
!           Function evaluation.
    
    IF (task == 'F' .OR. task == 'START') THEN
      CALL dtimer(times)
      ch = 'F'
      CALL dminfg(n, nx, ny, x, f, g, ch, prob, par)
      nfev = nfev + 1
      CALL dtimer(timef)
      fgtime = fgtime + (timef - times)
    END IF
    
!           Evaluate the gradient and the Hessian matrix.
    
    IF (task == 'GH' .OR. task == 'START') THEN
      CALL dtimer(times)
      ch = 'G'
      CALL dminfg(n, nx, ny, x, f, g, ch, prob, par)
      ngev = ngev + 1
      CALL dtimer(timef)
      fgtime = fgtime + (timef - times)
      
!              Evaluate the Hessian matrix.
      
      CALL dtimer(times)
      DO i = 1, n
        hs(i) = zero
        eta(i) = one
      END DO
      DO numgrp = 1, maxgrp
        DO j = 1, n
          IF (ngrp(j) == numgrp) hs(j) = one
        END DO
        CALL dminhs(n, nx, ny, x, hs, y, prob, par)
        CALL dsphesd(n, arow_ind, acol_ind, arow_ptr, acol_ptr,  &
                     listp, ngrp, maxgrp, numgrp, eta, y, a, adiag)
        DO j = 1, n
          IF (ngrp(j) == numgrp) hs(j) = zero
        END DO
      END DO
      nhev = nhev + 1
      nhsev = nhsev + maxgrp
      CALL dtimer(timef)
      htime = htime + (timef - times)
      
    END IF
    
!           Initialize the trust region bound.
    
    IF (task == 'START') THEN
      gnorm0 = dnrm2(n, g, 1)
      delta = dnrm2(n, g, 1)
    END IF
    
!           Stopping criteria.
    
    IF (task == 'GH' .OR. task == 'START') THEN
      gnorm = dgpnrm2(n, x, xl, xu, g)
      IF (gnorm <= gtol*gnorm0) THEN
        search = .false.
        task = 'CONVERGENCE: GTOL TEST SATISFIED'
      END IF
    END IF
    
    IF (nfev > maxfev) THEN
      search = .false.
      task = 'ERROR: NFEV > MAXFEV'
    END IF
    
!           Call the optimizer.
    
    IF (search) THEN
      CALL dtron(n, x, xl, xu, f, g, a, adiag, acol_ptr, arow_ind,  &
                 frtol, fatol, fmin, cgtol, itermax, delta, task,   &
                 b, bdiag, bcol_ptr, brow_ind,  &
                 l, ldiag, lcol_ptr, lrow_ind, iterscg)
    END IF
    
!           Exit search if the algorithm has converged.
    
    IF (task(1:4) == 'CONV') search = .false.
    
  END DO
  
  CALL dtimer(ttimef)
  
  WRITE (*, *) task
  
!        Summary information.
  
  nfree = 0
  nbind = 0
  DO i = 1, n
    IF (xl(i) < x(i) .AND. x(i) < xu(i)) THEN
      nfree = nfree + 1
    ELSE IF ((x(i) == xl(i) .AND. g(i) >= zero) .OR.  &
             (x(i) == xu(i) .AND. g(i) <= zero) .OR. (xl(i) == xu(i))) THEN
      nbind = nbind + 1
    END IF
  END DO
  
  ch = 'FG'
  CALL dminfg(n, nx, ny, x, f, g, ch, prob, par)
  gnorm = dgpnrm2(n, x, xl, xu, g)
  WRITE (nwrite,1000) prob, n, maxgrp, DBLE(nnz)/n, nfree, n-nfree, nbind,  &
                      nfev, ngev, nhev, nhsev, iterscg, f, gnorm
  
!        Timing information.
  
  ttime = ttimef - ttimes
  fgtime = 100. * fgtime/ttime
  htime  = 100. * htime/ttime
  
  WRITE (nwrite,2000) ttime, fgtime, htime, task
  
END DO
STOP

1000 FORMAT (' Problem ',  a6,                                    //  &
             ' Number of variables                         ', i12/   &
             ' Number of coloring groups                   ', i12/   &
             ' Average number of nonzeros in the strictly  '    /   &
             ' lower triangular part of the Hessian matrix ', f12.2/   &
             ' Number of free variables                    ', i12/   &
             ' Number of active variables                  ', i12/   &
             ' Number of binding variables                 ', i12/   &
             ' Number of function evaluations              ', i12/   &
             ' Number of gradient evaluations              ', i12/   &
             ' Number of Hessian evaluations               ', i12/   &
             ' Number of Hessian-vector evaluations        ', i12/   &
             ' Number of conjugate gradient iterations     ', i12 //   &
             ' Function value at final iterate          '   , g15.8/   &
             ' Projected gradient at final iterate      '   , g15.3 /)

2000 FORMAT (' Total execution time                        ', f12.2/   &
             ' Percentage in function evaluations          ', f12.1/   &
             ' Percentage in Hessian evaluations           ', f12.1 //   &
             ' Exit message     '                           , a60 /)

CONTAINS


SUBROUTINE dtimer(time)

REAL, INTENT(OUT) :: time

CALL CPU_TIME(time)
RETURN

END SUBROUTINE dtimer



SUBROUTINE dminfg(n, nx, ny, x, f, g, task, prob, par)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 15:49:09

INTEGER, INTENT(IN)                :: n
INTEGER, INTENT(IN)                :: nx
INTEGER, INTENT(IN)                :: ny
REAL (dp), INTENT(IN OUT)          :: x(:)
REAL (dp), INTENT(OUT)             :: f
REAL (dp), INTENT(OUT)             :: g(:)
CHARACTER (LEN=*), INTENT(IN OUT)  :: task
CHARACTER (LEN=*), INTENT(IN OUT)  :: prob
REAL (dp), INTENT(IN)              :: par

!  *********

!  Subroutine dminfg

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

!  The subroutine statement is

!    subroutine dminfg(n, nx, ny, x, f, g, task, prob, par)

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

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

!    f is a REAL (dp) variable.
!      On entry f need not be specified.
!      On exit f is set to the function evaluated at x if task = 'F' or 'FG'.

!    g is a REAL (dp) array of dimension n.
!      On entry g need not be specified.
!      On exit g contains the gradient evaluated at x if task = 'G' or 'FG'.

!    task is a character*60 variable.
!      On entry task specifies the action of the subroutine:

!         task               action
!         ----               ------
!          'F'     Evaluate the function at x.
!          'G'     Evaluate the gradient vector at x.
!          'FG'    Evaluate the function and the gradient at x.
!          'XS'    Set x to the standard starting point xs.

!      On exit task may be changed.
!         task can be changed in routine dljcfg.

!    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 probem-dependent parameter.
!      On exit par is unchanged.

!  Subprograms called

!    MINPACK-2 ... deptfg, dgl1fg, dgl2fg, dmsafg, dmsabc,
!                  dljcfg, dodcfg, dpjbfg, dsscfg

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

!  **********

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

! EXTERNAL deptfg, dgl1fg, dgl2fg, dmsafg, dmsabc, dljcfg, dodcfg, dpjbfg,  &
!          dsscfg

!     Select a problem.

SELECT CASE ( prob(1:4) )
  CASE ( 'DEPT' )
    CALL deptfg(nx, ny, x, f, g, task, par)
  CASE( 'DPJB' )
    CALL dpjbfg(nx, ny, x, f, g, task, par, ten)
  CASE( 'DMSA' )
    ALLOCATE( bottom(nx+2), top(nx+2), left(ny+2), right(ny+2) )
    CALL dmsabc(nx, ny, bottom, top, left, right)
    CALL dmsafg(nx, ny, x, f, g, task, bottom, top, left, right)
    DEALLOCATE( bottom, top, left, right )
  CASE( 'DODC' )
    CALL dodcfg(nx, ny, x, f, g, task, par)
  CASE( 'DSSC' )
    CALL dsscfg(nx, ny, x, f, g, task, par)
  CASE( 'DGL1' )
    CALL dgl1fg(n, x, f, g, task, par)
  CASE( 'DGL2' )
    CALL dgl2fg(nx, ny, x, f, g, task, INT(par))
  CASE( 'DLJ2' )
    CALL dljcfg(n, x, f, g, task, 2, n/2)
  CASE( 'DLJ3' )
    CALL dljcfg(n, x, f, g, task, 3, n/3)
  CASE DEFAULT
    prob = 'ERROR'
END SELECT

RETURN
END SUBROUTINE dminfg



SUBROUTINE dminhs(n, nx, ny, x, s, hs, prob, par)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 15:49:09

INTEGER, INTENT(IN)                :: n
INTEGER, INTENT(IN)                :: nx
INTEGER, INTENT(IN)                :: ny
REAL (dp), INTENT(IN OUT)          :: x(:)
REAL (dp), INTENT(IN OUT)          :: s(:)
REAL (dp), INTENT(OUT)             :: hs(:)
CHARACTER (LEN=*), INTENT(IN OUT)  :: prob
REAL (dp), INTENT(IN OUT)          :: par

!  **********

!  Subroutine dminhs

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

!  The subroutine statement is

!    subroutine dminhs(n, nx, ny, x, s, hs, prob, par)

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

!    x is a REAL (dp) array of dimension n.
!      On entry x specifies the vector x.
!      On exit x may be changed.
!         Routine dgl2hs can change x & s.