📄 t_lmder.f90
字号:
MODULE common_refnum
IMPLICIT NONE
! COMMON /refnum/ nprob,nfev,njev
INTEGER, SAVE :: nprob, nfev, njev
END MODULE common_refnum
PROGRAM test_lmder
! Code converted using TO_F90 by Alan Miller
! Date: 1999-12-09 Time: 23:25:02
! **********
! 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,LMDER1,SSQFCN
! FORTRAN-SUPPLIED ... SQRT
! 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) :: fjac(65,40), fnm(60), fvec(65), x(40)
! LOGICAL INPUT UNIT IS ASSUMED TO BE NUMBER 5.
! LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6.
INTEGER, PARAMETER :: nread = 5, nwrite = 6
! EXTERNAL fcn
INTERFACE
SUBROUTINE fcn(m, n, x, fvec, fjac, 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(:)
REAL (dp), INTENT(OUT) :: fjac(:,:)
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
SUBROUTINE ssqjac (m, n, x, fjac, nprob)
USE Levenberg_Marquardt
IMPLICIT NONE
INTEGER, INTENT(IN) :: m, n, nprob
REAL (dp), INTENT(IN) :: x(:)
REAL (dp), INTENT(OUT) :: fjac(:,:)
END SUBROUTINE ssqjac
END INTERFACE
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 lmder1 (fcn, m, n, x, fvec, fjac, 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
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 LMDER1'/)
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_lmder
SUBROUTINE fcn (m, n, x, fvec, fjac, 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(:)
REAL (dp), INTENT(OUT) :: fjac(:,:)
INTEGER, INTENT(IN OUT) :: iflag
INTERFACE
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
SUBROUTINE ssqjac (m, n, x, fjac, nprob)
USE Levenberg_Marquardt
IMPLICIT NONE
INTEGER, INTENT(IN) :: m, n, nprob
REAL (dp), INTENT(IN) :: x(:)
REAL (dp), INTENT(OUT) :: fjac(:,:)
END SUBROUTINE ssqjac
END INTERFACE
! **********
! 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 AND JACOBIAN SUBROUTINES SSQFCN AND SSQJAC WITH
! THE APPROPRIATE VALUE OF PROBLEM NUMBER (NPROB).
! SUBPROGRAMS CALLED
! MINPACK-SUPPLIED ... SSQFCN,SSQJAC
! ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.
! BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
! **********
IF (iflag == 1) CALL ssqfcn (m, n, x, fvec, nprob)
IF (iflag == 2) CALL ssqjac (m, n, x, fjac, nprob)
IF (iflag == 1) nfev = nfev + 1
IF (iflag == 2) njev = njev + 1
RETURN
! LAST CARD OF INTERFACE SUBROUTINE FCN.
END SUBROUTINE fcn
SUBROUTINE ssqjac (m, n, x, fjac, nprob)
USE Levenberg_Marquardt
IMPLICIT NONE
INTEGER, INTENT(IN) :: m, n, nprob
REAL (dp), INTENT(IN) :: x(:)
REAL (dp), INTENT(OUT) :: fjac(:,:)
! **********
! SUBROUTINE SSQJAC
! THIS SUBROUTINE DEFINES THE JACOBIAN MATRICES OF EIGHTEEN
! NONLINEAR LEAST SQUARES PROBLEMS. THE PROBLEM DIMENSIONS ARE
! AS DESCRIBED IN THE PROLOGUE COMMENTS OF SSQFCN.
! THE SUBROUTINE STATEMENT IS
! SUBROUTINE SSQJAC(M,N,X,FJAC,LDFJAC,NPROB)
! WHERE
! M AND N ARE POSITIVE INTEGER INPUT VARIABLES. N MUST NOT EXCEED M.
! X IS AN INPUT ARRAY OF LENGTH N.
! FJAC IS AN M BY N OUTPUT ARRAY WHICH CONTAINS THE JACOBIAN
! MATRIX OF THE NPROB FUNCTION EVALUATED AT X.
! LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN M
! WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC.
! NPROB IS A POSITIVE INTEGER VARIABLE WHICH DEFINES THE
! NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 18.
! SUBPROGRAMS CALLED
! FORTRAN-SUPPLIED ... DATAN,DCOS,DEXP,DSIN,SQRT
! ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.
! BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
! **********
INTEGER :: i, j, k, mm1, nm1
REAL (dp) :: div, dx, prod, s2, temp, ti, tmp1, tmp2, tmp3, tmp4, tpi
REAL (dp), PARAMETER :: v(11) = (/ 4.0D0, 2.0D0, 1.0D0, 5.0D-1, 2.5D-1, &
1.67D-1, 1.25D-1, 1.0D-1, 8.33D-2, &
7.14D-2, 6.25D-2 /)
REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp, two = 2.0_dp, &
three = 3.0_dp, four = 4.0_dp, five = 5.0_dp, &
eight = 8.0_dp, ten = 10.0_dp, c14 = 14.0_dp, &
c20 = 20.0_dp, c29 = 29.0_dp, c45 = 45.0_dp, &
c100 = 100.0_dp
! JACOBIAN ROUTINE SELECTOR.
SELECT CASE ( nprob )
CASE ( 1) ! LINEAR FUNCTION - FULL RANK.
temp = two / DBLE(m)
fjac(1:m,1:n) = -temp
DO j=1,n
fjac(j,j) = fjac(j,j) + one
END DO
CASE ( 2) ! LINEAR FUNCTION - RANK 1.
DO j=1,n
DO i=1,m
fjac(i,j) = i * j
END DO
END DO
CASE ( 3) ! LINEAR FUNCTION - RANK 1 WITH ZERO COLUMNS AND ROWS.
fjac(1:m,1:n) = zero
nm1 = n-1
mm1 = m-1
DO j=2,nm1
DO i=2,mm1
fjac(i,j) = (i-1) * j
END DO
END DO
CASE ( 4) ! ROSENBROCK FUNCTION.
fjac(1,1) = -c20*x(1)
fjac(1,2) = ten
fjac(2,1) = -one
fjac(2,2) = zero
CASE ( 5) ! HELICAL VALLEY FUNCTION.
tpi = eight*ATAN(one)
temp = x(1)**2 + x(2)**2
tmp1 = tpi*temp
tmp2 = SQRT(temp)
fjac(1,1) = c100*x(2)/tmp1
fjac(1,2) = -c100*x(1)/tmp1
fjac(1,3) = ten
fjac(2,1) = ten*x(1)/tmp2
fjac(2,2) = ten*x(2)/tmp2
fjac(2,3) = zero
fjac(3,1) = zero
fjac(3,2) = zero
fjac(3,3) = one
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -