lm.f90

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

!  subroutine lmder

!  the purpose of lmder is to minimize the sum of the squares of
!  m nonlinear functions in n variables by a modification of
!  the levenberg-marquardt algorithm. the user must provide a
!  subroutine which calculates the functions and the jacobian.

!  the subroutine statement is

!    subroutine lmder(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol,
!                     maxfev,diag,mode,factor,nprint,info,nfev,
!                     njev,ipvt,qtf,wa1,wa2,wa3,wa4)

!  where

!    fcn is the name of the user-supplied subroutine which
!      calculates the functions and the jacobian. fcn must
!      be declared in an external statement in the user
!      calling program, and should be written as follows.

!      subroutine fcn(m,n,x,fvec,fjac,ldfjac,iflag)
!      integer m,n,ldfjac,iflag
!      REAL (dp) x(:),fvec(m),fjac(ldfjac,n)
!      ----------
!      if iflag = 1 calculate the functions at x and
!      return this vector in fvec. do not alter fjac.
!      if iflag = 2 calculate the jacobian at x and
!      return this matrix in fjac.  Do not alter fvec.
!      ----------
!      return
!      end

!      the value of iflag should not be changed by fcn unless
!      the user wants to terminate execution of lmder.
!      in this case set iflag to a negative integer.

!    m is a positive integer input variable set to the number
!      of functions.

!    n is a positive integer input variable set to the number
!      of variables. n must not exceed m.

!    x is an array of length n. on input x must contain
!      an initial estimate of the solution vector. on output x
!      contains the final estimate of the solution vector.

!    fvec is an output array of length m which contains
!      the functions evaluated at the output x.

!    fjac is an output m by n array. the upper n by n submatrix
!      of fjac contains an upper triangular matrix r with
!      diagonal elements of nonincreasing magnitude such that

!             t     t           t
!            p *(jac *jac)*p = r *r

!      where p is a permutation matrix and jac is the final calculated
!      jacobian.  Column j of p is column ipvt(j) (see below) of the
!      identity matrix.  The lower trapezoidal part of fjac contains
!      information generated during the computation of r.

!    ldfjac is a positive integer input variable not less than m
!      which specifies the leading dimension of the array fjac.

!    ftol is a nonnegative input variable.  Termination occurs when both
!      the actual and predicted relative reductions in the sum of squares
!      are at most ftol.   Therefore, ftol measures the relative error
!      desired in the sum of squares.

!    xtol is a nonnegative input variable. termination
!      occurs when the relative error between two consecutive
!      iterates is at most xtol. therefore, xtol measures the
!      relative error desired in the approximate solution.

!    gtol is a nonnegative input variable.  Termination occurs when the
!      cosine of the angle between fvec and any column of the jacobian is
!      at most gtol in absolute value.  Therefore, gtol measures the
!      orthogonality desired between the function vector and the columns
!      of the jacobian.

!    maxfev is a positive integer input variable.  Termination occurs when
!      the number of calls to fcn with iflag = 1 has reached maxfev.

!    diag is an array of length n.  If mode = 1 (see below), diag is
!      internally set.  If mode = 2, diag must contain positive entries
!      that serve as multiplicative scale factors for the variables.

!    mode is an integer input variable.  if mode = 1, the
!      variables will be scaled internally.  if mode = 2,
!      the scaling is specified by the input diag.  other
!      values of mode are equivalent to mode = 1.

!    factor is a positive input variable used in determining the
!      initial step bound. this bound is set to the product of
!      factor and the euclidean norm of diag*x if nonzero, or else
!      to factor itself. in most cases factor should lie in the
!      interval (.1,100.).100. is a generally recommended value.

!    nprint is an integer input variable that enables controlled printing
!      of iterates if it is positive.  In this case, fcn is called with
!      iflag = 0 at the beginning of the first iteration and every nprint
!      iterations thereafter and immediately prior to return, with x, fvec,
!      and fjac available for printing.  fvec and fjac should not be
!      altered.  If nprint is not positive, no special calls of fcn with
!      iflag = 0 are made.

!    info is an integer output variable.  If the user has terminated
!      execution, info is set to the (negative) value of iflag.
!      See description of fcn.  Otherwise, info is set as follows.

!      info = 0  improper input parameters.

!      info = 1  both actual and predicted relative reductions
!                in the sum of squares are at most ftol.

!      info = 2  relative error between two consecutive iterates
!                is at most xtol.

!      info = 3  conditions for info = 1 and info = 2 both hold.

!      info = 4  the cosine of the angle between fvec and any column of
!                the jacobian is at most gtol in absolute value.

!      info = 5  number of calls to fcn with iflag = 1 has reached maxfev.

!      info = 6  ftol is too small.  No further reduction in
!                the sum of squares is possible.

!      info = 7  xtol is too small.  No further improvement in
!                the approximate solution x is possible.

!      info = 8  gtol is too small.  fvec is orthogonal to the
!                columns of the jacobian to machine precision.

!    nfev is an integer output variable set to the number of
!      calls to fcn with iflag = 1.

!    njev is an integer output variable set to the number of
!      calls to fcn with iflag = 2.

!    ipvt is an integer output array of length n.  ipvt
!      defines a permutation matrix p such that jac*p = q*r,
!      where jac is the final calculated jacobian, q is
!      orthogonal (not stored), and r is upper triangular
!      with diagonal elements of nonincreasing magnitude.
!      column j of p is column ipvt(j) of the identity matrix.

!    qtf is an output array of length n which contains
!      the first n elements of the vector (q transpose)*fvec.

!    wa1, wa2, and wa3 are work arrays of length n.

!    wa4 is a work array of length m.

!  subprograms called

!    user-supplied ...... fcn

!    minpack-supplied ... dpmpar,enorm,lmpar,qrfac

!    fortran-supplied ... ABS,MAX,MIN,SQRT,mod

!  argonne national laboratory. minpack project. march 1980.
!  burton s. garbow, kenneth e. hillstrom, jorge j. more

!  **********
INTEGER   :: i, iflag, iter, j, l
REAL (dp) :: actred, delta, dirder, epsmch, fnorm, fnorm1, gnorm,  &
             par, pnorm, prered, ratio, sum, temp, temp1, temp2, xnorm
REAL (dp) :: diag(n), qtf(n), wa1(n), wa2(n), wa3(n), wa4(m)
REAL (dp), PARAMETER :: one = 1.0_dp, p1 = 0.1_dp, p5 = 0.5_dp,  &
                        p25 = 0.25_dp, p75 = 0.75_dp, p0001 = 0.0001_dp, &
                        zero = 0.0_dp

!     epsmch is the machine precision.

epsmch = EPSILON(zero)

info = 0
iflag = 0
nfev = 0
njev = 0

!     check the input parameters for errors.

IF (n <= 0 .OR. m < n .OR. ftol < zero .OR. xtol < zero .OR. gtol < zero  &
    .OR. maxfev <= 0 .OR. factor <= zero) GO TO 300
IF (mode /= 2) GO TO 20
DO  j = 1, n
  IF (diag(j) <= zero) GO TO 300
END DO

!     evaluate the function at the starting point and calculate its norm.

20 iflag = 1
CALL fcn(m, n, x, fvec, fjac, iflag)
nfev = 1
IF (iflag < 0) GO TO 300
fnorm = enorm(m, fvec)

!     initialize levenberg-marquardt parameter and iteration counter.

par = zero
iter = 1

!     beginning of the outer loop.

!        calculate the jacobian matrix.

30 iflag = 2
CALL fcn(m, n, x, fvec, fjac, iflag)
njev = njev + 1
IF (iflag < 0) GO TO 300

!        if requested, call fcn to enable printing of iterates.

IF (nprint <= 0) GO TO 40
iflag = 0
IF (MOD(iter-1,nprint) == 0) CALL fcn(m, n, x, fvec, fjac, iflag)
IF (iflag < 0) GO TO 300

!        compute the qr factorization of the jacobian.

40 CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2)

!        on the first iteration and if mode is 1, scale according
!        to the norms of the columns of the initial jacobian.

IF (iter /= 1) GO TO 80
IF (mode == 2) GO TO 60
DO  j = 1, n
  diag(j) = wa2(j)
  IF (wa2(j) == zero) diag(j) = one
END DO

!        on the first iteration, calculate the norm of the scaled x
!        and initialize the step bound delta.

60 wa3(1:n) = diag(1:n)*x(1:n)
xnorm = enorm(n,wa3)
delta = factor*xnorm
IF (delta == zero) delta = factor

!        form (q transpose)*fvec and store the first n components in qtf.

80 wa4(1:m) = fvec(1:m)
DO  j = 1, n
  IF (fjac(j,j) == zero) GO TO 120
  sum = DOT_PRODUCT( fjac(j:m,j), wa4(j:m) )
  temp = -sum/fjac(j,j)
  DO  i = j, m
    wa4(i) = wa4(i) + fjac(i,j)*temp
  END DO
  120 fjac(j,j) = wa1(j)
  qtf(j) = wa4(j)
END DO

!        compute the norm of the scaled gradient.

gnorm = zero
IF (fnorm == zero) GO TO 170
DO  j = 1, n
  l = ipvt(j)
  IF (wa2(l) == zero) CYCLE
  sum = zero
  DO  i = 1, j
    sum = sum + fjac(i,j)*(qtf(i)/fnorm)
  END DO
  gnorm = MAX(gnorm,ABS(sum/wa2(l)))
END DO

!        test for convergence of the gradient norm.

170 IF (gnorm <= gtol) info = 4
IF (info /= 0) GO TO 300

!        rescale if necessary.

IF (mode == 2) GO TO 200
DO  j = 1, n
  diag(j) = MAX(diag(j), wa2(j))
END DO

!        beginning of the inner loop.

!           determine the levenberg-marquardt parameter.

200 CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2)

!           store the direction p and x + p. calculate the norm of p.

DO  j = 1, n
  wa1(j) = -wa1(j)
  wa2(j) = x(j) + wa1(j)
  wa3(j) = diag(j)*wa1(j)
END DO
pnorm = enorm(n, wa3)

!           on the first iteration, adjust the initial step bound.

IF (iter == 1) delta = MIN(delta,pnorm)

!           evaluate the function at x + p and calculate its norm.

iflag = 1
CALL fcn(m, n, wa2, wa4, fjac, iflag)
nfev = nfev + 1
IF (iflag < 0) GO TO 300
fnorm1 = enorm(m, wa4)

!           compute the scaled actual reduction.

actred = -one
IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2

!           compute the scaled predicted reduction and
!           the scaled directional derivative.

DO  j = 1, n
  wa3(j) = zero
  l = ipvt(j)
  temp = wa1(l)
  wa3(1:j) = wa3(1:j) + fjac(1:j,j)*temp
END DO
temp1 = enorm(n,wa3)/fnorm
temp2 = (SQRT(par)*pnorm)/fnorm
prered = temp1**2 + temp2**2/p5
dirder = -(temp1**2 + temp2**2)

!           compute the ratio of the actual to the predicted reduction.

ratio = zero
IF (prered /= zero) ratio = actred/prered

!           update the step bound.

IF (ratio <= p25) THEN
  IF (actred >= zero) temp = p5
  IF (actred < zero) temp = p5*dirder/(dirder + p5*actred)
  IF (p1*fnorm1 >= fnorm .OR. temp < p1) temp = p1
  delta = temp*MIN(delta, pnorm/p1)
  par = par/temp
ELSE
  IF (par /= zero .AND. ratio < p75) GO TO 260
  delta = pnorm/p5
  par = p5*par
END IF

!           test for successful iteration.

260 IF (ratio < p0001) GO TO 290

!           successful iteration. update x, fvec, and their norms.

DO  j = 1, n
  x(j) = wa2(j)
  wa2(j) = diag(j)*x(j)
END DO
fvec(1:m) = wa4(1:m)
xnorm = enorm(n,wa2)
fnorm = fnorm1
iter = iter + 1

!           tests for convergence.

290 IF (ABS(actred) <= ftol .AND. prered <= ftol .AND. p5*ratio <= one) info = 1
IF (delta <= xtol*xnorm) info = 2
IF (ABS(actred) <= ftol .AND. prered <= ftol  &
    .AND. p5*ratio <= one .AND. info == 2) info = 3
IF (info /= 0) GO TO 300

!           tests for termination and stringent tolerances.

IF (nfev >= maxfev) info = 5
IF (ABS(actred) <= epsmch .AND. prered <= epsmch  &
    .AND. p5*ratio <= one) info = 6
IF (delta <= epsmch*xnorm) info = 7
IF (gnorm <= epsmch) info = 8
IF (info /= 0) GO TO 300

!           end of the inner loop. repeat if iteration unsuccessful.

IF (ratio < p0001) GO TO 200

!        end of the outer loop.

GO TO 30

!     termination, either normal or user imposed.

300 IF (iflag < 0) info = iflag
iflag = 0
IF (nprint > 0) CALL fcn(m, n, x, fvec, fjac, iflag)
RETURN

!     last card of subroutine lmder.

END SUBROUTINE lmder



SUBROUTINE lmpar(n, r, ipvt, diag, qtb, delta, par, x, sdiag)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-12-09  Time: 12:46:12

! N.B. Arguments LDR, WA1 & WA2 have been removed.

INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: r(:,:)