t_lmdif.f90

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

!  PROBLEM NUMBER (NPROB).

!  SUBPROGRAMS CALLED

!    MINPACK-SUPPLIED ... SSQFCN

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

!  **********

! COMMON /refnum/ nprob,nfev,njev

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
END INTERFACE

CALL ssqfcn (m, n, x, fvec, nprob)
IF (iflag == 1) nfev = nfev + 1
IF (iflag == 2) njev = njev + 1
RETURN

!     LAST CARD OF INTERFACE SUBROUTINE FCN.

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

!     SUBROUTINE SSQFCN

!     THIS SUBROUTINE DEFINES THE FUNCTIONS OF EIGHTEEN NONLINEAR
!     LEAST SQUARES PROBLEMS. THE ALLOWABLE VALUES OF (M,N) FOR
!     FUNCTIONS 1,2 AND 3 ARE VARIABLE BUT WITH M >= N.
!     FOR FUNCTIONS 4,5,6,7,8,9 AND 10 THE VALUES OF (M,N) ARE
!     (2,2),(3,3),(4,4),(2,2),(15,3),(11,4) AND (16,3), RESPECTIVELY.
!     FUNCTION 11 (WATSON) HAS M = 31 WITH N USUALLY 6 OR 9.
!     HOWEVER, ANY N, N = 2,...,31, IS PERMITTED.
!     FUNCTIONS 12,13 AND 14 HAVE N = 3,2 AND 4, RESPECTIVELY, BUT
!     ALLOW ANY M >= N, WITH THE USUAL CHOICES BEING 10,10 AND 20.
!     FUNCTION 15 (CHEBYQUAD) ALLOWS M AND N VARIABLE WITH M >= N.
!     FUNCTION 16 (BROWN) ALLOWS N VARIABLE WITH M = N.
!     FOR FUNCTIONS 17 AND 18, THE VALUES OF (M,N) ARE
!     (33,5) AND (65,11), RESPECTIVELY.

!     THE SUBROUTINE STATEMENT IS

!       SUBROUTINE SSQFCN(M, N, X, FVEC, NPROB)

!     WHERE

!       M AND N ARE POSITIVE INTEGER INPUT VARIABLES. N MUST NOT
!         EXCEED M.

!       X IS AN INPUT ARRAY OF LENGTH N.

!       FVEC IS AN OUTPUT ARRAY OF LENGTH M WHICH CONTAINS THE NPROB
!         FUNCTION EVALUATED AT X.

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

!     SUBPROGRAMS CALLED

!       FORTRAN-SUPPLIED ... DATAN,DCOS,EXP,DSIN,SQRT,DSIGN

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

!     **********
INTEGER   :: i, iev, j, nm1
REAL (dp) :: div, dx, prod, total, s1, s2, temp, ti, tmp1, tmp2, tmp3, tmp4, tpi

REAL (dp), PARAMETER :: zero = 0.0_dp, zp25 = 0.25_dp, zp5 = 0.5_dp,  &
    one = 1.0_dp, two = 2.0_dp, five = 5.0_dp, eight = 8.0_dp, ten = 10._dp, &
    c13 = 13._dp, c14 = 14._dp, c29 = 29._dp, c45 = 45._dp
REAL (dp), PARAMETER :: v(11) = (/  &
    4.0D0, 2.0D0, 1.0D0, 5.0E-1_dp, 2.5E-1_dp, 1.67E-1_dp, 1.25E-1_dp,  &
    1.0E-1_dp, 8.33E-2_dp, 7.14E-2_dp, 6.25E-2_dp /)
REAL (dp), PARAMETER :: y1(15) = (/  &
    1.4E-1_dp, 1.8E-1_dp, 2.2E-1_dp, 2.5E-1_dp, 2.9E-1_dp, 3.2E-1_dp,  &
    3.5E-1_dp, 3.9E-1_dp, 3.7E-1_dp, 5.8E-1_dp, 7.3E-1_dp, 9.6E-1_dp,  &
    1.34_dp, 2.1_dp, 4.39_dp /)
REAL (dp), PARAMETER :: y2(11) = (/   &
    1.957E-1_dp, 1.947E-1_dp, 1.735E-1_dp, 1.6E-1_dp, 8.44E-2_dp, 6.27E-2_dp, &
    4.56E-2_dp, 3.42E-2_dp, 3.23E-2_dp, 2.35E-2_dp, 2.46E-2_dp /)
REAL (dp), PARAMETER :: y3(16) = (/   &
    3.478D4, 2.861D4, 2.365D4, 1.963D4, 1.637D4, 1.372D4, 1.154D4, 9.744D3,  &
    8.261D3, 7.03D3, 6.005D3, 5.147D3, 4.427D3, 3.82D3, 3.307D3, 2.872D3 /)
REAL (dp), PARAMETER :: y4(33) = (/   &
    8.44E-1_dp, 9.08E-1_dp, 9.32E-1_dp, 9.36E-1_dp, 9.25E-1_dp, 9.08E-1_dp,  &
    8.81E-1_dp, 8.5E-1_dp, 8.18E-1_dp, 7.84E-1_dp, 7.51E-1_dp, 7.18E-1_dp,   &
    6.85E-1_dp, 6.58E-1_dp, 6.28E-1_dp, 6.03E-1_dp, 5.8E-1_dp, 5.58E-1_dp,   &
    5.38E-1_dp, 5.22E-1_dp, 5.06E-1_dp, 4.9E-1_dp, 4.78E-1_dp, 4.67E-1_dp,   &
    4.57E-1_dp, 4.48E-1_dp, 4.38E-1_dp, 4.31E-1_dp, 4.24E-1_dp, 4.2E-1_dp,   &
    4.14E-1_dp, 4.11E-1_dp, 4.06E-1_dp /)
REAL (dp), PARAMETER :: y5(65) = (/   &
    1.366_dp, 1.191_dp, 1.112_dp, 1.013_dp, 9.91E-1_dp, 8.85E-1_dp, 8.31E-1_dp,   &
    8.47E-1_dp, 7.86E-1_dp, 7.25E-1_dp, 7.46E-1_dp, 6.79E-1_dp, 6.08E-1_dp,  &
    6.55E-1_dp, 6.16E-1_dp, 6.06E-1_dp, 6.02E-1_dp, 6.26E-1_dp, 6.51E-1_dp,  &
    7.24E-1_dp, 6.49E-1_dp, 6.49E-1_dp, 6.94E-1_dp, 6.44E-1_dp, 6.24E-1_dp,  &
    6.61E-1_dp, 6.12E-1_dp, 5.58E-1_dp, 5.33E-1_dp, 4.95E-1_dp, 5.0E-1_dp,   &
    4.23E-1_dp, 3.95E-1_dp, 3.75E-1_dp, 3.72E-1_dp, 3.91E-1_dp, 3.96E-1_dp,  &
    4.05E-1_dp, 4.28E-1_dp, 4.29E-1_dp, 5.23E-1_dp, 5.62E-1_dp, 6.07E-1_dp,  &
    6.53E-1_dp, 6.72E-1_dp, 7.08E-1_dp, 6.33E-1_dp, 6.68E-1_dp, 6.45E-1_dp,  &
    6.32E-1_dp, 5.91E-1_dp, 5.59E-1_dp, 5.97E-1_dp, 6.25E-1_dp, 7.39E-1_dp,  &
    7.1E-1_dp, 7.29E-1_dp, 7.2E-1_dp, 6.36E-1_dp, 5.81E-1_dp, 4.28E-1_dp,    &
    2.92E-1_dp, 1.62E-1_dp,  9.8E-2_dp, 5.4E-2_dp /)

!     FUNCTION ROUTINE SELECTOR.

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

    total = SUM( x(1:n) )
    temp = two*total/DBLE(m) + one
    DO  i=1,m
      fvec(i) = -temp
      IF (i <= n) fvec(i) = fvec(i) + x(i)
    END DO

  CASE (    2)     !     LINEAR FUNCTION - RANK 1.

    total = zero
    DO  j=1,n
      total = total + DBLE(j)*x(j)
    END DO
    DO  i=1,m
      fvec(i) = DBLE(i)*total - one
    END DO

  CASE (    3)     !     LINEAR FUNCTION - RANK 1 WITH ZERO COLUMNS AND ROWS.

    total = zero
    nm1 = n-1
    DO  j=2,nm1
      total = total + DBLE(j)*x(j)
    END DO
    DO  i=1,m
      fvec(i) = DBLE(i-1)*total - one
    END DO
    fvec(m) = -one

  CASE (    4)     !     ROSENBROCK FUNCTION.

    fvec(1) = ten*(x(2) - x(1)**2)
    fvec(2) = one - x(1)

  CASE (    5)     !     HELICAL VALLEY FUNCTION.

    tpi = eight*ATAN(one)
    tmp1 = SIGN(zp25,x(2))
    IF (x(1) > zero) tmp1 = ATAN(x(2)/x(1))/tpi
    IF (x(1) < zero) tmp1 = ATAN(x(2)/x(1))/tpi + zp5
    tmp2 = SQRT(x(1)**2 + x(2)**2)
    fvec(1) = ten*(x(3) - ten*tmp1)
    fvec(2) = ten*(tmp2 - one)
    fvec(3) = x(3)

  CASE (    6)     !     POWELL SINGULAR FUNCTION.

    fvec(1) = x(1) + ten*x(2)
    fvec(2) = SQRT(five)*(x(3) - x(4))
    fvec(3) = (x(2) - two*x(3))**2
    fvec(4) = SQRT(ten)*(x(1) - x(4))**2

  CASE (    7)     !     FREUDENSTEIN AND ROTH FUNCTION.

    fvec(1) = -c13 + x(1) + ((five - x(2))*x(2) - two)*x(2)
    fvec(2) = -c29 + x(1) + ((one + x(2))*x(2) - c14)*x(2)

  CASE (    8)     !     BARD FUNCTION.

    DO  i=1,15
      tmp1 = DBLE(i)
      tmp2 = DBLE(16-i)
      tmp3 = tmp1
      IF (i > 8) tmp3 = tmp2
      fvec(i) = y1(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3))
    END DO

  CASE (    9)     !     KOWALIK AND OSBORNE FUNCTION.

    DO  i=1,11
      tmp1 = v(i)*(v(i) + x(2))
      tmp2 = v(i)*(v(i) + x(3)) + x(4)
      fvec(i) = y2(i) - x(1)*tmp1/tmp2
    END DO

  CASE (   10)     !     MEYER FUNCTION.

    DO  i=1,16
      temp = five*DBLE(i) + c45 + x(3)
      tmp1 = x(2)/temp
      tmp2 = EXP(tmp1)
      fvec(i) = x(1)*tmp2 - y3(i)
    END DO

  CASE (   11)     !     WATSON FUNCTION.

    DO  i=1,29
      div = DBLE(i)/c29
      s1 = zero
      dx = one
      DO  j=2,n
        s1 = s1 + DBLE(j-1)*dx*x(j)
        dx = div*dx
      END DO
      s2 = zero
      dx = one
      DO  j=1,n
        s2 = s2 + dx*x(j)
        dx = div*dx
      END DO
      fvec(i) = s1 - s2**2 - one
    END DO
    fvec(30) = x(1)
    fvec(31) = x(2) - x(1)**2 - one

  CASE (   12)     !     BOX 3-DIMENSIONAL FUNCTION.

    DO  i=1,m
      temp = DBLE(i)
      tmp1 = temp/ten
      fvec(i) = EXP(-tmp1*x(1)) - EXP(-tmp1*x(2)) + (EXP(-temp) -  &
                EXP(- tmp1))*x(3)
    END DO

  CASE (   13)     !     JENNRICH AND SAMPSON FUNCTION.

    DO  i=1,m
      temp = DBLE(i)
      fvec(i) = two + two*temp - EXP(temp*x(1)) - EXP(temp*x(2))
    END DO

  CASE (   14)     !     BROWN AND DENNIS FUNCTION.

    DO  i=1,m
      temp = DBLE(i)/five
      tmp1 = x(1) + temp*x(2) - EXP(temp)
      tmp2 = x(3) + SIN(temp)*x(4) - COS(temp)
      fvec(i) = tmp1**2 + tmp2**2
    END DO

  CASE (   15)     !     CHEBYQUAD FUNCTION.

    fvec(1:m) = zero
    DO  j=1,n
      tmp1 = one
      tmp2 = two*x(j) - one
      temp = two*tmp2
      DO  i=1,m
        fvec(i) = fvec(i) + tmp2
        ti = temp*tmp2 - tmp1
        tmp1 = tmp2
        tmp2 = ti
      END DO
    END DO
    dx = one/DBLE(n)
    iev = -1
    DO  i=1,m
      fvec(i) = dx*fvec(i)
      IF (iev > 0) fvec(i) = fvec(i) + one/(DBLE(i)**2 - one)
      iev = -iev
    END DO

  CASE (   16)     !     BROWN ALMOST-LINEAR FUNCTION.

    total = -DBLE(n+1)
    prod = one
    DO  j=1,n
      total = total + x(j)
      prod = x(j)*prod
    END DO
    DO  i=1,n
      fvec(i) = x(i) + total
    END DO
    fvec(n) = prod - one

  CASE (   17)     !     OSBORNE 1 FUNCTION.

DO  i=1,33
  temp = ten*DBLE(i-1)
  tmp1 = EXP(-x(4)*temp)
  tmp2 = EXP(-x(5)*temp)
  fvec(i) = y4(i) - (x(1) + x(2)*tmp1 + x(3)*tmp2)
END DO

  CASE (   18)     !     OSBORNE 2 FUNCTION.

    DO  i=1,65
      temp = DBLE(i-1)/ten
      tmp1 = EXP(-x(5)*temp)
      tmp2 = EXP(-x(6)*(temp-x(9))**2)
      tmp3 = EXP(-x(7)*(temp-x(10))**2)
      tmp4 = EXP(-x(8)*(temp-x(11))**2)
      fvec(i) = y5(i) - (x(1)*tmp1 + x(2)*tmp2 + x(3)*tmp3 + x(4)*tmp4)
    END DO

END SELECT

RETURN

!     LAST CARD OF SUBROUTINE SSQFCN.

END SUBROUTINE ssqfcn