dmain.f90

牛顿优化算法源fortran代码

!    row_ptr is an integer array of dimension n + 1.
!      On entry row_ptr must contain pointers to the rows of A.
!         The nonzeros in row j of A must be in positions
!         row_ptr(j), ... , row_ptr(j+1) - 1 of col_ind.
!      On exit row_ptr 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 of row_ind.
!      On exit col_ptr is unchanged.

!    listp is an integer array of length n.
!      On input listp need not be specified.
!      On output listp specifies a permutation of the matrix.
!         Element (i,j) of the matrix is the (listp(i), listp(j))
!         element of the permuted matrix.

!    ngrp is an integer array of length n.
!      On entry ngrp need not be specified.
!      On exit ngrp specifies the partition of the columns of A.
!         Column j belongs to group ngrp(j).

!    maxgrp is an integer variable.
!      On entry maxgrp need not be specified.  WRONG.
!      On exit maxgrp specifies the number of groups in
!         the partition of the columns of A.

!    numgrp is an integer variable.
!      On input numgrp must be set to a group number.
!      On output numgrp is unchanged.

!    eta is a REAL (dp) variable.
!      On input eta is the difference parameter vector.
!      On output eta is unchanged.

!    fhesd is a REAL (dp) array of length n.
!      On input fhesd contains an approximation to A*d, where A
!        is the Hessian matrix and d is the difference vector for group numgrp.
!      On output fhesd is unchanged.

!    fhes is a REAL (dp) array of length nnz.
!      On input fhes need not be specified.
!      On output fhes is overwritten. When numgrp = maxgrp the
!         array fhes contains an approximation to the Hessian matrix
!         in compressed column storage. The elements in column j of
!         the strict lower triangular part of the Hessian matrix are

!            fhes(k), k = col_ptr(j),...,col_ptr(j+1)-1,

!         and the row indices for these elements are

!            row_ind(k), k = col_ptr(j),...,col_ptr(j+1)-1.

!    diag is a REAL (dp) array of length n.
!      On input diag need not be specified.
!      On output diag is overwritten.  When numgrp = maxgrp the array diag
!         contains the diagonal elements of an approximation to the Hessian
!         matrix.

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

!  **********

INTEGER   :: i, ip, irow, j, jp, k, l, numg, numl
INTEGER   :: iwk(n)
REAL (dp) :: sum

!     Store the i-th element of gradient difference fhesd corresponding to
!     group numgrp if there is a position (i,j) such that ngrp(j) = numgrp and
!     (i,j) is mapped onto the lower triangular part of the permuted matrix.

DO j = 1, n
  IF (ngrp(j) == numgrp) THEN
    diag(j) = fhesd(j)/eta(j)
    numl = listp(j)
    DO ip = row_ptr(j), row_ptr(j+1)-1
      i = col_ind(ip)
      IF (listp(i) > numl) THEN
        DO jp = col_ptr(i), col_ptr(i+1)-1
          IF (row_ind(jp) == j) THEN
            fhes(jp) = fhesd(i)
            EXIT
          END IF
        END DO
      END IF
    END DO
    DO jp = col_ptr(j), col_ptr(j+1)-1
      i = row_ind(jp)
      IF (listp(i) >= numl) fhes(jp) = fhesd(i)
    END DO
  END IF
END DO

!     Exit if this is not the last group.

IF (numgrp < maxgrp) RETURN

!     Mark all column indices j such that (i,j) is mapped onto
!     the lower triangular part of the permuted matrix.

DO i = 1, n
  numl = listp(i)
  DO ip = row_ptr(i), row_ptr(i+1)-1
    j = col_ind(ip)
    IF (numl >= listp(j)) col_ind(ip) = -col_ind(ip)
  END DO
  DO jp = col_ptr(i), col_ptr(i+1)-1
    j = row_ind(jp)
    IF (numl > listp(j)) row_ind(jp) = -row_ind(jp)
  END DO
END DO

!     Invert the array listp.

DO j = 1, n
  iwk(listp(j)) = j
END DO
listp(1:n) = iwk(1:n)

!     Determine the lower triangular part of the original matrix.

DO irow = n, 1, -1
  i = listp(irow)
  
!        Find the positions of the elements in the i-th row of the lower
!        triangular part of the original matrix that have already been
!        determined.
  
  DO ip = row_ptr(i), row_ptr(i+1)-1
    j = col_ind(ip)
    IF (j > 0) THEN
      DO jp = col_ptr(j), col_ptr(j+1)-1
        IF (row_ind(jp) == i) THEN
          iwk(j) = jp
          EXIT
        END IF
      END DO
    END IF
  END DO
  
!        Determine the elements in the i-th row of the lower
!        triangular part of the original matrix which get mapped
!        onto the lower triangular part of the permuted matrix.
  
  DO k = row_ptr(i), row_ptr(i+1)-1
    j = -col_ind(k)
    IF (j > 0) THEN
      col_ind(k) = j
      
!              Determine the (i,j) element.
      
      numg = ngrp(j)
      sum = zero
      DO ip = row_ptr(i), row_ptr(i+1)-1
        l = ABS(col_ind(ip))
        IF (ngrp(l) == numg .AND. l /= j) sum = sum + fhes(iwk(l))*eta(l)
      END DO
      DO jp = col_ptr(i), col_ptr(i+1)-1
        l = ABS(row_ind(jp))
        IF (ngrp(l) == numg .AND. l /= j) sum = sum + fhes(jp)*eta(l)
      END DO
      
!              Store the (i,j) element.
      
      DO jp = col_ptr(j), col_ptr(j+1)-1
        IF (row_ind(jp) == i) THEN
          fhes(jp) = (fhes(jp) - sum)/eta(j)
          EXIT
        END IF
      END DO
    END IF
  END DO
  
  
!        Determine the elements in the i-th row of the strict upper
!        triangular part of the original matrix which get mapped
!        onto the lower triangular part of the permuted matrix.
  
  DO k = col_ptr(i), col_ptr(i+1)-1
    j = -row_ind(k)
    IF (j > 0) THEN
      row_ind(k) = j
      
!              Determine the (i,j) element.
      
      numg = ngrp(j)
      sum = zero
      DO ip = row_ptr(i), row_ptr(i+1)-1
        l = ABS(col_ind(ip))
        IF (ngrp(l) == numg) sum = sum + fhes(iwk(l))*eta(l)
      END DO
      DO jp = col_ptr(i), col_ptr(i+1)-1
        l = ABS(row_ind(jp))
        IF (ngrp(l) == numg .AND. l /= j) sum = sum + fhes(jp)*eta(l)
      END DO
      
!              Store the (i,j) element.
      
      fhes(k) = (fhes(k) - sum)/eta(j)
    END IF
  END DO
END DO

!     Re-invert the array listp.

DO j = 1, n
  iwk(listp(j)) = j
END DO
listp(1:n) = iwk(1:n)

RETURN
END SUBROUTINE dsphesd



SUBROUTINE dsetsp(n, nnz, row_ind, col_ind, col_ptr, row_ptr,  &
                  method, info, listp, ngrp, maxgrp)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 11:18:32

INTEGER, INTENT(IN)      :: n
INTEGER, INTENT(OUT)     :: nnz
INTEGER, INTENT(OUT)     :: row_ind(:)
INTEGER, INTENT(IN OUT)  :: col_ind(:)
INTEGER, INTENT(IN OUT)  :: col_ptr(:)    ! col_ptr(n+1)
INTEGER, INTENT(IN OUT)  :: row_ptr(:)    ! row_ptr(n+1)
INTEGER, INTENT(IN)      :: method
INTEGER, INTENT(OUT)     :: info
INTEGER, INTENT(OUT)     :: listp(:)
INTEGER, INTENT(OUT)     :: ngrp(:)
INTEGER, INTENT(OUT)     :: maxgrp

!  **********

!  Subroutine dsetsp

!  Given the non-zero elements of a symmetric matrix A in coordinate
!  format, this subroutine computes the coloring information for
!  determining A from a lower triangular substitution method.

!  The sparsity pattern of the matrix A is specified by the
!  arrays row_ind and col_ind.  On input the indices for the
!  non-zero elements in the lower triangular part of A are

!        (row_ind(k),col_ind(k)), k = 1,2,...,nnz.

!  The (row_ind(k),col_ind(k)) pairs may be specified in any order.
!  Duplicate input pairs are permitted, but they are eliminated.
!  Diagonal elements must be part of the sparsity pattern.
!  Any pair (row_ind(k),col_ind(k)), where row_ind(k) is less than
!  col_ind(k), is replaced by the pair (col_ind(k),row_ind(k)).

!  The input coordinate format is changed to a storage format.
!  On output the strict lower triangular part of A is stored in both
!  compressed column storage and compressed row storage.

!  The information required to determine the matrix from a lower
!  triangular substitution method is obtained from subroutine dssm.

!  The subroutine statement is

!    subroutine dsetsp(n, nnz, row_ind, col_ind, col_ptr, row_ptr,
!                      method, info, listp, ngrp, maxgrp, iwa)

!  where

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

!    nnz is an integer variable.
!      On entry nnz is the number of non-zeros entries in the
!        coordinate format.
!      On exit nnz is the number of non-zeroes in the strict
!        lower triangular part of the matrix A.

!    row_ind is an integer array of length nnz.
!      On entry row_ind must contain the row indices of the non-zero
!         elements of A in coordinate format.
!      On exit row_ind contains row indices for the strict
!         lower triangular part of A in compressed column storage.

!    col_ind is an integer array of length nnz.
!      On entry col_ind must contain the column indices of the
!         non-zero elements of A in coordinate format.
!      On exit col_ind contains column indices for the strict
!         lower triangular part of A in compressed row storage.

!    row_ptr is an integer array of length n + 1.
!      On entry row_ptr need not be specified.
!      On exit row_ptr must contain pointers to the rows of A.
!         The non-zeros in row j of A must be in positions
!         row_ptr(j), ... , row_ptr(j+1) - 1 of col_ind.

!    col_ptr is an integer array of length n + 1.
!      On entry col_ptr need not be specified.
!      On exit col_ptr must contain pointers to the columns of A.
!         The non-zeros in column j of A must be in positions
!         col_ptr(j), ... , col_ptr(j+1) - 1 of row_ind.

!    method is an integer variable.
!      On input with method = 1, the direct method is used to
!        determine the partition and symmetric permutation.
!        Otherwise, the indirect method is used.

!    info is an integer variable.
!      On input info need not be specified.
!      On output info is set as follows.

!         info = 1  Normal termination.

!         info = 0  Input n or nnz are not positive.

!         info < 0  There is an error in the sparsity pattern.
!                   If k = -info then row_ind(k) or col_ind(k) is
!                   not an integer between 1 and n, or the k-th
!                   diagonal element is not in the sparsity pattern.

!    listp is an integer array of length n.
!      On input listp need not be specified.
!      On output listp specifies a permutation of the matrix.
!         Element (i,j) of the matrix is the (listp(i),listp(j))
!         element of the permuted matrix.

!    ngrp is an integer array of length n.
!      On entry ngrp need not be specified.
!      On exit ngrp specifies the partition of the columns
!         of A. Column j belongs to group ngrp(j).

!    maxgrp is an integer variable.
!      On entry maxgrp need not be specified.
!      On exit maxgrp specifies the number of groups in
!         the partition of the columns of A.

!  Subprograms called

!    MINPACK-2  ......  dssm

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

!  **********

INTEGER :: i, j, mingrp, ndiag

!     Subroutine dssm first checks the sparsity structure.
!     If there are no errors, the input format is changed to
!     compressed column storage. Finally, the information required
!     to determine the matrix from matrix-vector products is obtained.

CALL dssm(n, nnz, row_ind, col_ind, method, listp,  &
          ngrp, maxgrp, mingrp, info, row_ptr, col_ptr)

!     Exit if there are errors on input.

IF (info <= 0) RETURN

!     Change the sparsity structure to exclude the diagonal entries.

ndiag = 0
DO j = 1, n
  DO i = col_ptr(j), col_ptr(j+1) - 1
    IF (row_ind(i) == j) THEN
      ndiag = ndiag + 1
    ELSE
      row_ind(i-ndiag) =  row_ind(i)
    END IF
  END DO
  col_ptr(j) = col_ptr(j) - (j - 1)
END DO
col_ptr(n+1) = col_ptr(n+1) - n
nnz = nnz - n

RETURN
END SUBROUTINE dsetsp



SUBROUTINE deptfg(nx, ny, x, f, fgrad, task, c)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 15:49:07

INTEGER, INTENT(IN)            :: nx
INTEGER, INTENT(IN)            :: ny
REAL (dp), INTENT(IN OUT)      :: x(:)    ! x(nx*ny)
REAL (dp), INTENT(OUT)         :: f
REAL (dp), INTENT(OUT)         :: fgrad(:)    ! fgrad(nx*ny)
CHARACTER (LEN=*), INTENT(IN)  :: task
REAL (dp), INTENT(IN)          :: c

!  **********

!  Subroutine deptfg

!  This subroutine computes the function and gradient of the
!  elastic-plastic torsion problem.

!  The subroutine statement is

!    subroutine deptfg(nx, ny, x, f, fgrad, task, c)

!  where

!    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.
!      On exit ny is unchanged.

!    x is a REAL (dp) array of dimension nx*ny.
!      On entry x specifies the vector x if task = 'F', 'G', or 'FG'.
!         Otherwise x need not be specified.
!      On exit x is unchanged if task = 'F', 'G', or 'FG'.  Otherwise
!         x is set according to task.

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

!    fgrad is a REAL (dp) array of dimension nx*ny.
!      On entry fgrad need not be specified.
!      On exit fgrad 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 at x.
!          'FG'    Evaluate the function and the gradient at x.
!          'XS'    Set x to the standard starting point xs.
!          'XL'    Set x to the lower bound xl.
!          'XU'    Set x to the upper bound xu.

!      On exit task is unchanged.

!    c is a REAL (dp) variable.
!      On entry c is the angle of twist per unit length.
!      On exit c is unchanged.

!  MINPACK-2 Project. March 1999.
!  Argonne National Laboratory and University of Minnesota.