t_lmdif.f90

开发的lm算法,很有用的一种优化算法. 对非线性优化有很大用处

MODULE common_refnum
IMPLICIT NONE

! COMMON /refnum/ nprob,nfev,njev
INTEGER, SAVE :: nprob, nfev, njev

END MODULE common_refnum



PROGRAM test_lmdif
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-12-11  Time: 00:05:31
 
!  **********

!  THIS PROGRAM TESTS CODES FOR THE LEAST-SQUARES SOLUTION OF
!  M NONLINEAR EQUATIONS IN N VARIABLES. IT CONSISTS OF A DRIVER
!  AND AN INTERFACE SUBROUTINE FCN. THE DRIVER READS IN DATA,
!  CALLS THE NONLINEAR LEAST-SQUARES SOLVER, AND FINALLY PRINTS
!  OUT INFORMATION ON THE PERFORMANCE OF THE SOLVER. THIS IS
!  ONLY A SAMPLE DRIVER, MANY OTHER DRIVERS ARE POSSIBLE. THE
!  INTERFACE SUBROUTINE FCN IS NECESSARY TO TAKE INTO ACCOUNT THE
!  FORMS OF CALLING SEQUENCES USED BY THE FUNCTION AND JACOBIAN
!  SUBROUTINES IN THE VARIOUS NONLINEAR LEAST-SQUARES SOLVERS.

!  SUBPROGRAMS CALLED

!    USER-SUPPLIED ...... FCN

!    MINPACK-SUPPLIED ... DPMPAR,ENORM,INITPT,LMDIF1,SSQFCN

!    FORTRAN-SUPPLIED ... DSQRT

!  ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.
!  BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE

!  **********

USE Levenberg_Marquardt
USE common_refnum
IMPLICIT NONE
INTEGER   :: i, ic, info, k, m, n, ntries
INTEGER   :: iwa(40), ma(60), na(60), nf(60), nj(60), np(60), nx(60)
REAL (dp) :: factor, fnorm1, fnorm2, tol
REAL (dp) :: fnm(60), fvec(65), x(40)

! EXTERNAL fcn

INTERFACE
  SUBROUTINE fcn(m, n, x, fvec, iflag)
    IMPLICIT NONE
    INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
    INTEGER, INTENT(IN)        :: m, n
    REAL (dp), INTENT(IN)      :: x(:)
    REAL (dp), INTENT(IN OUT)  :: fvec(:)
    INTEGER, INTENT(IN OUT)    :: iflag
  END SUBROUTINE fcn

  SUBROUTINE ssqfcn (m, n, x, fvec, nprob)
    USE Levenberg_Marquardt
    IMPLICIT NONE
    INTEGER, INTENT(IN)     :: m, n
    REAL (dp), INTENT(IN)   :: x(:)
    REAL (dp), INTENT(OUT)  :: fvec(:)
    INTEGER, INTENT(IN)     :: nprob
  END SUBROUTINE ssqfcn
END INTERFACE

! COMMON /refnum/ nprob,nfev,njev

!     LOGICAL INPUT UNIT IS ASSUMED TO BE NUMBER 5.
!     LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6.

INTEGER, PARAMETER   :: nread = 5, nwrite = 6

REAL (dp), PARAMETER :: one = 1.0_dp, ten = 10.0_dp

tol = SQRT( EPSILON(one) )
ic = 0
10 READ (nread,50) nprob, n, m, ntries
IF (nprob <= 0) GO TO 30
factor = one
DO  k=1,ntries
  ic = ic+1
  CALL initpt (n, x, nprob, factor)
  CALL ssqfcn (m, n, x, fvec, nprob)
  fnorm1 = enorm(m, fvec)
  WRITE (nwrite,60) nprob, n, m
  nfev = 0
  njev = 0
  CALL lmdif1 (fcn, m, n, x, fvec, tol, info, iwa)
  CALL ssqfcn (m, n, x, fvec, nprob)
  fnorm2 = enorm(m,fvec)
  np(ic) = nprob
  na(ic) = n
  ma(ic) = m
  nf(ic) = nfev
  njev = njev/n
  nj(ic) = njev
  nx(ic) = info
  fnm(ic) = fnorm2
  WRITE (nwrite,70) fnorm1, fnorm2, nfev, njev, info, x(1:n)
  factor = ten*factor
END DO
GO TO 10

30 WRITE (nwrite,80) ic
WRITE (nwrite,90)
DO  i=1,ic
  WRITE (nwrite,100) np(i), na(i), ma(i), nf(i), nj(i), nx(i), fnm(i)
END DO
STOP
50 FORMAT (4I5)
60 FORMAT (////'      PROBLEM', i5, '      DIMENSIONS', 2I5//)
70 FORMAT ('      INITIAL L2 NORM OF THE RESIDUALS', g15.7//  &
           '      FINAL L2 NORM OF THE RESIDUALS  ', g15.7//  &
           '      NUMBER OF FUNCTION EVALUATIONS  ', i10//  &
           '      NUMBER OF JACOBIAN EVALUATIONS  ', i10//  &
           '      EXIT PARAMETER', t39, i10//  &
           '      FINAL APPROXIMATE SOLUTION'// (t6, 5g15.7))
80 FORMAT (' SUMMARY OF ', i3, ' CALLS TO LMDIF1'/)
90 FORMAT (' NPROB   N    M   NFEV  NJEV  INFO  FINAL L2 NORM'/)
100 FORMAT (3I5, 3I6, ' ', g15.7)


CONTAINS


SUBROUTINE initpt (n, x, nprob, factor)
INTEGER, INTENT(IN)     :: n, nprob
REAL (dp), INTENT(IN)   :: factor
REAL (dp), INTENT(OUT)  :: x(:)

!  **********

!  SUBROUTINE INITPT

!  THIS SUBROUTINE SPECIFIES THE STANDARD STARTING POINTS FOR THE
!  FUNCTIONS DEFINED BY SUBROUTINE SSQFCN. THE SUBROUTINE RETURNS
!  IN X A MULTIPLE (FACTOR) OF THE STANDARD STARTING POINT. FOR
!  THE 11TH FUNCTION THE STANDARD STARTING POINT IS ZERO, SO IN
!  THIS CASE, IF FACTOR IS NOT UNITY, THEN THE SUBROUTINE RETURNS
!  THE VECTOR  X(J) = FACTOR, J=1,...,N.

!  THE SUBROUTINE STATEMENT IS

!    SUBROUTINE INITPT(N,X,NPROB,FACTOR)

!  WHERE

!    N IS A POSITIVE INTEGER INPUT VARIABLE.

!    X IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS THE STANDARD
!      STARTING POINT FOR PROBLEM NPROB MULTIPLIED BY FACTOR.

!    NPROB IS A POSITIVE INTEGER INPUT VARIABLE WHICH DEFINES THE
!      NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18.

!    FACTOR IS AN INPUT VARIABLE WHICH SPECIFIES THE MULTIPLE OF
!      THE STANDARD STARTING POINT.  IF FACTOR IS UNITY, NO
!      MULTIPLICATION IS PERFORMED.

!  ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.
!  BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE

!  **********

INTEGER   :: j
REAL (dp) :: h
REAL (dp), PARAMETER :: zero = 0.0_dp, half = 0.5_dp, one = 1.0_dp,  &
                        two = 2.0_dp, three = 3.0_dp, five = 5.0_dp, &
                        seven = 7.0_dp, ten = 10.0_dp, twenty = 20.0_dp, &
                        twntf = 25.0_dp
REAL (dp), PARAMETER :: c1 = 1.2_dp, c2 = 0.25_dp, c3 = 0.39_dp,    &
                        c4 = 0.415_dp, c5 = 0.02_dp, c6 = 4000._dp, &
                        c7 = 250._dp, c8 = 0.3_dp, c9 = 0.4_dp,     &
                        c10 = 1.5_dp, c11 = 0.01_dp, c12 = 1.3_dp,  &
                        c13 = 0.65_dp, c14 = 0.7_dp, c15 = 0.6_dp,  &
                        c16 = 4.5_dp, c17 = 5.5_dp

!     SELECTION OF INITIAL POINT.

SELECT CASE ( nprob )
  CASE (    1:3)   !     LINEAR FUNCTION - FULL RANK OR RANK 1.

    x(1:n) = one

  CASE (    4)     !     ROSENBROCK FUNCTION.

    x(1) = -c1
    x(2) = one

  CASE (    5)     !     HELICAL VALLEY FUNCTION.

    x(1) = -one
    x(2) = zero
    x(3) = zero

  CASE (    6)     !     POWELL SINGULAR FUNCTION.

    x(1) = three
    x(2) = -one
    x(3) = zero
    x(4) = one

  CASE (    7)     !     FREUDENSTEIN AND ROTH FUNCTION.

x(1) = half
x(2) = -two

  CASE (    8)     !     BARD FUNCTION.

x(1) = one
x(2) = one
x(3) = one

  CASE (    9)     !     KOWALIK AND OSBORNE FUNCTION.

x(1) = c2
x(2) = c3
x(3) = c4
x(4) = c3

  CASE (   10)     !     MEYER FUNCTION.

    x(1) = c5
    x(2) = c6
    x(3) = c7

  CASE (   11)     !     WATSON FUNCTION.

    x(1:n) = zero

  CASE (   12)     !     BOX 3-DIMENSIONAL FUNCTION.

    x(1) = zero
    x(2) = ten
    x(3) = twenty

  CASE (   13)     !     JENNRICH AND SAMPSON FUNCTION.

    x(1) = c8
    x(2) = c9

  CASE (   14)     !     BROWN AND DENNIS FUNCTION.

    x(1) = twntf
    x(2) = five
    x(3) = -five
    x(4) = -one

  CASE (   15)     !     CHEBYQUAD FUNCTION.

    h = one / DBLE(n+1)
    DO  j=1,n
      x(j) = DBLE(j)*h
    END DO

  CASE (   16)     !     BROWN ALMOST-LINEAR FUNCTION.

    x(1:n) = half

  CASE (   17)     !     OSBORNE 1 FUNCTION.

    x(1) = half
    x(2) = c10
    x(3) = -one
    x(4) = c11
    x(5) = c5

  CASE (   18)     !     OSBORNE 2 FUNCTION.

    x(1) = c12
    x(2) = c13
    x(3) = c13
    x(4) = c14
    x(5) = c15
    x(6) = three
    x(7) = five
    x(8) = seven
    x(9) = two
    x(10) = c16
    x(11) = c17

END SELECT

!     COMPUTE MULTIPLE OF INITIAL POINT.

IF (factor == one) GO TO 250
IF (nprob == 11) GO TO 230
x(1:n) = factor*x(1:n)
GO TO 250

230 x(1:n) = factor

250 RETURN

!     LAST CARD OF SUBROUTINE INITPT.

END SUBROUTINE initpt

!     LAST CARD OF DRIVER.

END PROGRAM test_lmdif



SUBROUTINE fcn (m, n, x, fvec, iflag)
USE common_refnum
IMPLICIT NONE
INTEGER, PARAMETER         :: dp = SELECTED_REAL_KIND(12, 60)
INTEGER, INTENT(IN)        :: m, n
REAL (dp), INTENT(IN)      :: x(:)
REAL (dp), INTENT(IN OUT)  :: fvec(:)
INTEGER, INTENT(IN OUT)    :: iflag

!  **********

!  THE CALLING SEQUENCE OF FCN SHOULD BE IDENTICAL TO THE
!  CALLING SEQUENCE OF THE FUNCTION SUBROUTINE IN THE NONLINEAR
!  LEAST-SQUARES SOLVER. FCN SHOULD ONLY CALL THE TESTING
!  FUNCTION SUBROUTINE SSQFCN WITH THE APPROPRIATE VALUE OF