lm.f90

module lm
use precision, only : fdp

! Levenberg_Marquardt MINPACK routines which are used by both LMDIF & LMDER
! under BSD-like license
!

implicit none

integer, parameter :: lmdp = fdp

private


public :: lmdp, lmdif1, lmdif, lmder1, lmder, enorm
public :: decompactify, compactify


contains


  function decompactify(p,pmin,pmax,ndim)
    implicit none
    integer, intent(in) :: ndim
    real(fdp), dimension(ndim) :: decompactify
    real(fdp), dimension(ndim), intent(in) :: p,pmin,pmax
    real(fdp), dimension(ndim) :: t

    integer :: i

    if (any(p.gt.pmax).or.any(p.lt.pmin)) then
       stop 'decompactify: parameter out of bounds!'
    endif

    t = (p - 0.5_fdp*(pmin+pmax)) / (0.5_fdp*(pmax-pmin))

    if (any(abs(t).ge.1._fdp)) stop 'decompactify: pb!'

    do i=1,ndim
       if (t(i).eq.1._fdp) t(i)=1._fdp-epsilon(1._fdp)
       if (t(i).eq.-1._fdp) t(i) = -1._fdp + epsilon(1._fdp)
    enddo

    decompactify = t/(1._fdp-abs(t))

  end function decompactify


  function compactify(pinf,pmin,pmax,ndim)
    implicit none
    integer, intent(in) :: ndim
    real(fdp), dimension(ndim) :: compactify
    real(fdp), dimension(ndim),intent(in) :: pinf,pmin,pmax
    real(fdp), dimension(ndim) :: t

    t = pinf/(1._fdp+abs(pinf))

    if (any(abs(t).ge.1._fdp)) stop 'compactify pb!'

    compactify = 0.5_fdp*(pmin+pmax) + 0.5_fdp*(pmax-pmin)*t

  end function compactify




  
  SUBROUTINE lmdif1(fcn1, m, n, x, fvec, tol, info, iwa)

    INTEGER, INTENT(IN)        :: m
    INTEGER, INTENT(IN)        :: n
    REAL (lmdp), INTENT(IN OUT)  :: x(:)
    REAL (lmdp), INTENT(OUT)     :: fvec(:)
    REAL (lmdp), INTENT(IN)      :: tol
    INTEGER, INTENT(OUT)       :: info
    INTEGER, INTENT(OUT)       :: iwa(:)




    include 'lm.h'

    !  **********

    !  subroutine lmdif1

    !  The purpose of lmdif1 is to minimize the sum of the squares of m nonlinear
    !  functions in n variables by a modification of the Levenberg-Marquardt
    !  algorithm.  This is done by using the more general least-squares
    !  solver lmdif.  The user must provide a subroutine which calculates the
    !  functions.  The jacobian is then calculated by a forward-difference
    !  approximation.

    !  the subroutine statement is

    !    subroutine lmdif1(fcn, m, n, x, fvec, tol, info, iwa)

    !  where

    !    fcn is the name of the user-supplied subroutine which calculates
    !      the functions.  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, iflag)
    !      integer m, n, iflag
    !      REAL (lmdp) x(n), fvec(m)
    !      ----------
    !      calculate the functions at x and return this vector in fvec.
    !      ----------
    !      return
    !      end

    !      the value of iflag should not be changed by fcn unless
    !      the user wants to terminate execution of lmdif1.
    !      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.

    !    tol is a nonnegative input variable.  Termination occurs when the
    !      algorithm estimates either that the relative error in the sum of
    !      squares is at most tol or that the relative error between x and the
    !      solution is at most tol.

    !    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  algorithm estimates that the relative error
    !                in the sum of squares is at most tol.

    !      info = 2  algorithm estimates that the relative error
    !                between x and the solution is at most tol.

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

    !      info = 4  fvec is orthogonal to the columns of the
    !                jacobian to machine precision.

    !      info = 5  number of calls to fcn has reached or exceeded 200*(n+1).

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

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

    !    iwa is an integer work array of length n.

    !    wa is a work array of length lwa.

    !    lwa is a positive integer input variable not less than m*n+5*n+m.

    !  subprograms called

    !    user-supplied ...... fcn

    !    minpack-supplied ... lmdif

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

    !  **********
    INTEGER   :: maxfev, mode, nfev, nprint
    REAL (lmdp) :: epsfcn, ftol, gtol, xtol, fjac(m,n)
    REAL (lmdp), PARAMETER :: factor = 100._lmdp, zero = 0.0_lmdp

    info = 0

    !     check the input parameters for errors.

    IF (n <= 0 .OR. m < n .OR. tol < zero) GO TO 10

    !     call lmdif.

    maxfev = 300*(n + 1)
    ftol = tol
    xtol = tol
    gtol = zero
    epsfcn = zero
    mode = 1
    nprint = 0
    CALL lmdif(fcn1, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn,   &
         mode, factor, nprint, info, nfev, fjac, iwa)
    IF (info == 8) info = 4

10  RETURN

    !     last card of subroutine lmdif1.

  END SUBROUTINE lmdif1



  SUBROUTINE lmdif(fcn1, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn,  &
       mode, factor, nprint, info, nfev, fjac, ipvt)


    INTEGER, INTENT(IN)        :: m
    INTEGER, INTENT(IN)        :: n
    REAL (lmdp), INTENT(IN OUT)  :: x(:)
    REAL (lmdp), INTENT(OUT)     :: fvec(:)
    REAL (lmdp), INTENT(IN)      :: ftol
    REAL (lmdp), INTENT(IN)      :: xtol
    REAL (lmdp), INTENT(IN OUT)  :: gtol
    INTEGER, INTENT(IN OUT)    :: maxfev
    REAL (lmdp), INTENT(IN OUT)  :: epsfcn
    INTEGER, INTENT(IN)        :: mode
    REAL (lmdp), INTENT(IN)      :: factor
    INTEGER, INTENT(IN)        :: nprint
    INTEGER, INTENT(OUT)       :: info
    INTEGER, INTENT(OUT)       :: nfev
    REAL (lmdp), INTENT(OUT)     :: fjac(:,:)    ! fjac(ldfjac,n)
    INTEGER, INTENT(OUT)       :: ipvt(:)

    include 'lm.h'


    !  **********

    !  subroutine lmdif

    !  The purpose of lmdif 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.  The jacobian is then calculated by a forward-difference
    !  approximation.

    !  the subroutine statement is

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

    ! N.B. 7 of these arguments have been removed in this version.

    !  where

    !    fcn is the name of the user-supplied subroutine which calculates the
    !      functions.  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, iflag)
    !      integer m, n, iflag
    !      REAL (lmdp) x(:), fvec(m)
    !      ----------
    !      calculate the functions at x and return this vector in fvec.
    !      ----------
    !      return
    !      end

    !      the value of iflag should not be changed by fcn unless
    !      the user wants to terminate execution of lmdif.
    !      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.

    !    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 is at least maxfev by the end of an iteration.

    !    epsfcn is an input variable used in determining a suitable step length
    !      for the forward-difference approximation.  This approximation assumes
    !      that the relative errors in the functions are of the order of epsfcn.
    !      If epsfcn is less than the machine precision, it is assumed that the
    !      relative errors in the functions are of the order of the machine
    !      precision.

    !    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 and fvec available
    !      for printing.  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 <= 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 has reached or exceeded 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.

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

    !    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,fdjac2,lmpar,qrfac

    !    fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod

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

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

    !     epsmch is the machine precision.

    epsmch = EPSILON(zero)

    info = 0
    iflag = 0
    nfev = 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 fcn1(m, n, x, fvec, 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 fdjac2(fcn1, m, n, x, fvec, fjac, iflag, epsfcn)
    nfev = nfev + n
    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 fcn1(m, n, x, fvec, 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 fcn1(m, n, wa2, wa4, 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)
       DO  i = 1, j
          wa3(i) = wa3(i) + fjac(i,j)*temp
       END DO
    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 fcn1(m, n, x, fvec, iflag)
    RETURN

    !     last card of subroutine lmdif.

  END SUBROUTINE lmdif



  SUBROUTINE lmder1(fcn2, m, n, x, fvec, fjac, tol, info, ipvt)


    INTEGER, INTENT(IN)        :: m
    INTEGER, INTENT(IN)        :: n
    REAL (lmdp), INTENT(IN OUT)  :: x(:)
    REAL (lmdp), INTENT(OUT)     :: fvec(:)
    REAL (lmdp), INTENT(IN OUT)  :: fjac(:,:)    ! fjac(ldfjac,n)
    REAL (lmdp), INTENT(IN)      :: tol
    INTEGER, INTENT(OUT)       :: info
    INTEGER, INTENT(IN OUT)    :: ipvt(:)


    include 'lm.h'


    !  **********

    !  subroutine lmder1

    !  The purpose of lmder1 is to minimize the sum of the squares of
    !  m nonlinear functions in n variables by a modification of the
    !  levenberg-marquardt algorithm.  This is done by using the more
    !  general least-squares solver lmder.  The user must provide a
    !  subroutine which calculates the functions and the jacobian.

    !  the subroutine statement is

    !    subroutine lmder1(fcn, m, n, x, fvec, fjac, tol, info, ipvt)

    !  where

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

    !      subroutine fcn(m, n, x, fvec, fjac, iflag)
    !      integer   :: m, n, ldfjac, iflag
    !      REAL (lmdp) :: x(:), fvec(:), fjac(:,:)
    !      ----------
    !      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 lmder1.
    !      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.

    !    tol is a nonnegative input variable. termination occurs
    !      when the algorithm estimates either that the relative
    !      error in the sum of squares is at most tol or that
    !      the relative error between x and the solution is at most tol.

    !    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  algorithm estimates that the relative error
    !                in the sum of squares is at most tol.

    !      info = 2  algorithm estimates that the relative error
    !                between x and the solution is at most tol.

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

    !      info = 4  fvec is orthogonal to the columns of the
    !                jacobian to machine precision.

    !      info = 5  number of calls to fcn with iflag = 1 has reached 100*(n+1).

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

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

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

    !    wa is a work array of length lwa.

    !    lwa is a positive integer input variable not less than 5*n+m.

    !  subprograms called

    !    user-supplied ...... fcn

    !    minpack-supplied ... lmder

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

    !  **********
    INTEGER   :: maxfev, mode, nfev, njev, nprint
    REAL (lmdp) :: ftol, gtol, xtol
    REAL (lmdp), PARAMETER :: factor = 100._lmdp, zero = 0.0_lmdp

    info = 0

    !     check the input parameters for errors.

    IF ( n <= 0 .OR. m < n .OR. tol < zero ) GO TO 10

    !     call lmder.

    maxfev = 100*(n + 1)
    ftol = tol
    xtol = tol
    gtol = zero
    mode = 1
    nprint = 0
    CALL lmder(fcn2, m, n, x, fvec, fjac, ftol, xtol, gtol, maxfev,  &
         mode, factor, nprint, info, nfev, njev, ipvt)
    IF (info == 8) info = 4

10  RETURN

    !     last card of subroutine lmder1.

  END SUBROUTINE lmder1



  SUBROUTINE lmder(fcn2, m, n, x, fvec, fjac, ftol, xtol, gtol, maxfev, &
       mode, factor, nprint, info, nfev, njev, ipvt)


    INTEGER, INTENT(IN)        :: m
    INTEGER, INTENT(IN)        :: n
    REAL (lmdp), INTENT(IN OUT)  :: x(:)
    REAL (lmdp), INTENT(OUT)     :: fvec(m)
    REAL (lmdp), INTENT(OUT)     :: fjac(:,:)    ! fjac(ldfjac,n)
    REAL (lmdp), INTENT(IN)      :: ftol
    REAL (lmdp), INTENT(IN)      :: xtol
    REAL (lmdp), INTENT(IN OUT)  :: gtol
    INTEGER, INTENT(IN OUT)    :: maxfev
    INTEGER, INTENT(IN)        :: mode
    REAL (lmdp), INTENT(IN)      :: factor
    INTEGER, INTENT(IN)        :: nprint
    INTEGER, INTENT(OUT)       :: info
    INTEGER, INTENT(OUT)       :: nfev
    INTEGER, INTENT(OUT)       :: njev
    INTEGER, INTENT(OUT)       :: ipvt(:)

    include 'lm.h'


    !  **********

    !  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 (lmdp) 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 (lmdp) :: actred, delta, dirder, epsmch, fnorm, fnorm1, gnorm,  &
         par, pnorm, prered, ratio, sum, temp, temp1, temp2, xnorm
    REAL (lmdp) :: diag(n), qtf(n), wa1(n), wa2(n), wa3(n), wa4(m)
    REAL (lmdp), PARAMETER :: one = 1.0_lmdp, p1 = 0.1_lmdp, p5 = 0.5_lmdp,  &
         p25 = 0.25_lmdp, p75 = 0.75_lmdp, p0001 = 0.0001_lmdp, &
         zero = 0.0_lmdp

    !     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 fcn2(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 fcn2(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 fcn2(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 fcn2(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 fcn2(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)


    INTEGER, INTENT(IN)        :: n
    REAL (lmdp), INTENT(IN OUT)  :: r(:,:)
    INTEGER, INTENT(IN)        :: ipvt(:)
    REAL (lmdp), INTENT(IN)      :: diag(:)
    REAL (lmdp), INTENT(IN)      :: qtb(:)
    REAL (lmdp), INTENT(IN)      :: delta
    REAL (lmdp), INTENT(OUT)     :: par
    REAL (lmdp), INTENT(OUT)     :: x(:)
    REAL (lmdp), INTENT(OUT)     :: sdiag(:)

    !  **********

    !  subroutine lmpar

    !  Given an m by n matrix a, an n by n nonsingular diagonal matrix d,
    !  an m-vector b, and a positive number delta, the problem is to determine a
    !  value for the parameter par such that if x solves the system

    !        a*x = b ,     sqrt(par)*d*x = 0 ,

    !  in the least squares sense, and dxnorm is the euclidean
    !  norm of d*x, then either par is zero and

    !        (dxnorm-delta) <= 0.1*delta ,

    !  or par is positive and

    !        abs(dxnorm-delta) <= 0.1*delta .

    !  This subroutine completes the solution of the problem if it is provided
    !  with the necessary information from the r factorization, with column
    !  qpivoting, of a.  That is, if a*p = q*r, where p is a permutation matrix,
    !  q has orthogonal columns, and r is an upper triangular matrix with diagonal
    !  elements of nonincreasing magnitude, then lmpar expects the full upper
    !  triangle of r, the permutation matrix p, and the first n components of
    !  (q transpose)*b.
    !  On output lmpar also provides an upper triangular matrix s such that

    !         t   t                   t
    !        p *(a *a + par*d*d)*p = s *s .

    !  s is employed within lmpar and may be of separate interest.

    !  Only a few iterations are generally needed for convergence of the algorithm.
    !  If, however, the limit of 10 iterations is reached, then the output par
    !  will contain the best value obtained so far.

    !  the subroutine statement is

    !    subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, wa1,wa2)

    !  where

    !    n is a positive integer input variable set to the order of r.

    !    r is an n by n array. on input the full upper triangle
    !      must contain the full upper triangle of the matrix r.
    !      On output the full upper triangle is unaltered, and the
    !      strict lower triangle contains the strict upper triangle
    !      (transposed) of the upper triangular matrix s.

    !    ldr is a positive integer input variable not less than n
    !      which specifies the leading dimension of the array r.

    !    ipvt is an integer input array of length n which defines the
    !      permutation matrix p such that a*p = q*r. column j of p
    !      is column ipvt(j) of the identity matrix.

    !    diag is an input array of length n which must contain the
    !      diagonal elements of the matrix d.

    !    qtb is an input array of length n which must contain the first
    !      n elements of the vector (q transpose)*b.

    !    delta is a positive input variable which specifies an upper
    !      bound on the euclidean norm of d*x.

    !    par is a nonnegative variable. on input par contains an
    !      initial estimate of the levenberg-marquardt parameter.
    !      on output par contains the final estimate.

    !    x is an output array of length n which contains the least
    !      squares solution of the system a*x = b, sqrt(par)*d*x = 0,
    !      for the output par.

    !    sdiag is an output array of length n which contains the
    !      diagonal elements of the upper triangular matrix s.

    !    wa1 and wa2 are work arrays of length n.

    !  subprograms called

    !    minpack-supplied ... dpmpar,enorm,qrsolv

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

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

    !  **********
    INTEGER   :: iter, j, jm1, jp1, k, l, nsing
    REAL (lmdp) :: dxnorm, dwarf, fp, gnorm, parc, parl, paru, sum, temp
    REAL (lmdp) :: wa1(n), wa2(n)
    REAL (lmdp), PARAMETER :: p1 = 0.1_lmdp, p001 = 0.001_lmdp, zero = 0.0_lmdp

    !     dwarf is the smallest positive magnitude.

    dwarf = TINY(zero)

    !     compute and store in x the gauss-newton direction. if the
    !     jacobian is rank-deficient, obtain a least squares solution.

    nsing = n
    DO  j = 1, n
       wa1(j) = qtb(j)
       IF (r(j,j) == zero .AND. nsing == n) nsing = j - 1
       IF (nsing < n) wa1(j) = zero
    END DO

    DO  k = 1, nsing
       j = nsing - k + 1
       wa1(j) = wa1(j)/r(j,j)
       temp = wa1(j)
       jm1 = j - 1
       wa1(1:jm1) = wa1(1:jm1) - r(1:jm1,j)*temp
    END DO

    DO  j = 1, n
       l = ipvt(j)
       x(l) = wa1(j)
    END DO

    !     initialize the iteration counter.
    !     evaluate the function at the origin, and test
    !     for acceptance of the gauss-newton direction.

    iter = 0
    wa2(1:n) = diag(1:n)*x(1:n)
    dxnorm = enorm(n, wa2)
    fp = dxnorm - delta
    IF (fp <= p1*delta) GO TO 220

    !     if the jacobian is not rank deficient, the newton
    !     step provides a lower bound, parl, for the zero of
    !     the function.  Otherwise set this bound to zero.

    parl = zero
    IF (nsing < n) GO TO 120
    DO  j = 1, n
       l = ipvt(j)
       wa1(j) = diag(l)*(wa2(l)/dxnorm)
    END DO
    DO  j = 1, n
       sum = DOT_PRODUCT( r(1:j-1,j), wa1(1:j-1) )
       wa1(j) = (wa1(j) - sum)/r(j,j)
    END DO
    temp = enorm(n,wa1)
    parl = ((fp/delta)/temp)/temp

    !     calculate an upper bound, paru, for the zero of the function.

120 DO  j = 1, n
       sum = DOT_PRODUCT( r(1:j,j), qtb(1:j) )
       l = ipvt(j)
       wa1(j) = sum/diag(l)
    END DO
    gnorm = enorm(n,wa1)
    paru = gnorm/delta
    IF (paru == zero) paru = dwarf/MIN(delta,p1)

    !     if the input par lies outside of the interval (parl,paru),
    !     set par to the closer endpoint.

    par = MAX(par,parl)
    par = MIN(par,paru)
    IF (par == zero) par = gnorm/dxnorm

    !     beginning of an iteration.

150 iter = iter + 1

    !        evaluate the function at the current value of par.

    IF (par == zero) par = MAX(dwarf, p001*paru)
    temp = SQRT(par)
    wa1(1:n) = temp*diag(1:n)
    CALL qrsolv(n, r, ipvt, wa1, qtb, x, sdiag)
    wa2(1:n) = diag(1:n)*x(1:n)
    dxnorm = enorm(n, wa2)
    temp = fp
    fp = dxnorm - delta

    !        if the function is small enough, accept the current value
    !        of par. also test for the exceptional cases where parl
    !        is zero or the number of iterations has reached 10.

    IF (ABS(fp) <= p1*delta .OR. parl == zero .AND. fp <= temp  &
         .AND. temp < zero .OR. iter == 10) GO TO 220

    !        compute the newton correction.

    DO  j = 1, n
       l = ipvt(j)
       wa1(j) = diag(l)*(wa2(l)/dxnorm)
    END DO
    DO  j = 1, n
       wa1(j) = wa1(j)/sdiag(j)
       temp = wa1(j)
       jp1 = j + 1
       wa1(jp1:n) = wa1(jp1:n) - r(jp1:n,j)*temp
    END DO
    temp = enorm(n,wa1)
    parc = ((fp/delta)/temp)/temp

    !        depending on the sign of the function, update parl or paru.

    IF (fp > zero) parl = MAX(parl,par)
    IF (fp < zero) paru = MIN(paru,par)

    !        compute an improved estimate for par.

    par = MAX(parl, par+parc)

    !        end of an iteration.

    GO TO 150

    !     termination.

220 IF (iter == 0) par = zero
    RETURN

    !     last card of subroutine lmpar.

  END SUBROUTINE lmpar



  SUBROUTINE qrfac(m, n, a, pivot, ipvt, rdiag, acnorm)


    INTEGER, INTENT(IN)        :: m
    INTEGER, INTENT(IN)        :: n
    REAL (lmdp), INTENT(IN OUT)  :: a(:,:)
    LOGICAL, INTENT(IN)        :: pivot
    INTEGER, INTENT(OUT)       :: ipvt(:)
    REAL (lmdp), INTENT(OUT)     :: rdiag(:)
    REAL (lmdp), INTENT(OUT)     :: acnorm(:)

    !  **********

    !  subroutine qrfac

    !  This subroutine uses Householder transformations with column pivoting
    !  (optional) to compute a qr factorization of the m by n matrix a.
    !  That is, qrfac determines an orthogonal matrix q, a permutation matrix p,
    !  and an upper trapezoidal matrix r with diagonal elements of nonincreasing
    !  magnitude, such that a*p = q*r.  The householder transformation for
    !  column k, k = 1,2,...,min(m,n), is of the form

    !                        t
    !        i - (1/u(k))*u*u

    !  where u has zeros in the first k-1 positions.  The form of this
    !  transformation and the method of pivoting first appeared in the
    !  corresponding linpack subroutine.

    !  the subroutine statement is

    !    subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa)

    ! N.B. 3 of these arguments have been omitted in this version.

    !  where

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

    !    n is a positive integer input variable set to the number of columns of a.

    !    a is an m by n array.  On input a contains the matrix for
    !      which the qr factorization is to be computed.  On output
    !      the strict upper trapezoidal part of a contains the strict
    !      upper trapezoidal part of r, and the lower trapezoidal
    !      part of a contains a factored form of q (the non-trivial
    !      elements of the u vectors described above).

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

    !    pivot is a logical input variable.  If pivot is set true,
    !      then column pivoting is enforced.  If pivot is set false,
    !      then no column pivoting is done.

    !    ipvt is an integer output array of length lipvt.  ipvt
    !      defines the permutation matrix p such that a*p = q*r.
    !      Column j of p is column ipvt(j) of the identity matrix.
    !      If pivot is false, ipvt is not referenced.

    !    lipvt is a positive integer input variable.  If pivot is false,
    !      then lipvt may be as small as 1.  If pivot is true, then
    !      lipvt must be at least n.

    !    rdiag is an output array of length n which contains the
    !      diagonal elements of r.

    !    acnorm is an output array of length n which contains the norms of the
    !      corresponding columns of the input matrix a.
    !      If this information is not needed, then acnorm can coincide with rdiag.

    !    wa is a work array of length n.  If pivot is false, then wa
    !      can coincide with rdiag.

    !  subprograms called

    !    minpack-supplied ... dpmpar,enorm

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

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

    !  **********
    INTEGER   :: i, j, jp1, k, kmax, minmn
    REAL (lmdp) :: ajnorm, epsmch, sum, temp, wa(n)
    REAL (lmdp), PARAMETER :: one = 1.0_lmdp, p05 = 0.05_lmdp, zero = 0.0_lmdp

    !     epsmch is the machine precision.

    epsmch = EPSILON(zero)

    !     compute the initial column norms and initialize several arrays.

    DO  j = 1, n
       acnorm(j) = enorm(m,a(1:,j))
       rdiag(j) = acnorm(j)
       wa(j) = rdiag(j)
       IF (pivot) ipvt(j) = j
    END DO

    !     Reduce a to r with Householder transformations.

    minmn = MIN(m,n)
    DO  j = 1, minmn
       IF (.NOT.pivot) GO TO 40

       !        Bring the column of largest norm into the pivot position.

       kmax = j
       DO  k = j, n
          IF (rdiag(k) > rdiag(kmax)) kmax = k
       END DO
       IF (kmax == j) GO TO 40
       DO  i = 1, m
          temp = a(i,j)
          a(i,j) = a(i,kmax)
          a(i,kmax) = temp
       END DO
       rdiag(kmax) = rdiag(j)
       wa(kmax) = wa(j)
       k = ipvt(j)
       ipvt(j) = ipvt(kmax)
       ipvt(kmax) = k

       !     Compute the Householder transformation to reduce the
       !     j-th column of a to a multiple of the j-th unit vector.

40     ajnorm = enorm(m-j+1, a(j:,j))
       IF (ajnorm == zero) CYCLE
       IF (a(j,j) < zero) ajnorm = -ajnorm
       a(j:m,j) = a(j:m,j)/ajnorm
       a(j,j) = a(j,j) + one

       !     Apply the transformation to the remaining columns and update the norms.

       jp1 = j + 1
       DO  k = jp1, n
          sum = DOT_PRODUCT( a(j:m,j), a(j:m,k) )
          temp = sum/a(j,j)
          a(j:m,k) = a(j:m,k) - temp*a(j:m,j)
          IF (.NOT.pivot .OR. rdiag(k) == zero) CYCLE
          temp = a(j,k)/rdiag(k)
          rdiag(k) = rdiag(k)*SQRT(MAX(zero, one-temp**2))
          IF (p05*(rdiag(k)/wa(k))**2 > epsmch) CYCLE
          rdiag(k) = enorm(m-j, a(jp1:,k))
          wa(k) = rdiag(k)
       END DO
       rdiag(j) = -ajnorm
    END DO
    RETURN

    !     last card of subroutine qrfac.

  END SUBROUTINE qrfac



  SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag)


    INTEGER, INTENT(IN)        :: n
    REAL (lmdp), INTENT(IN OUT)  :: r(:,:)
    INTEGER, INTENT(IN)        :: ipvt(:)
    REAL (lmdp), INTENT(IN)      :: diag(:)
    REAL (lmdp), INTENT(IN)      :: qtb(:)
    REAL (lmdp), INTENT(OUT)     :: x(:)
    REAL (lmdp), INTENT(OUT)     :: sdiag(:)

    !  **********

    !  subroutine qrsolv

    !  Given an m by n matrix a, an n by n diagonal matrix d, and an m-vector b,
    !  the problem is to determine an x which solves the system

    !        a*x = b ,     d*x = 0 ,

    !  in the least squares sense.

    !  This subroutine completes the solution of the problem if it is provided
    !  with the necessary information from the qr factorization, with column
    !  pivoting, of a.  That is, if a*p = q*r, where p is a permutation matrix,
    !  q has orthogonal columns, and r is an upper triangular matrix with diagonal
    !  elements of nonincreasing magnitude, then qrsolv expects the full upper
    !  triangle of r, the permutation matrix p, and the first n components of
    !  (q transpose)*b.  The system a*x = b, d*x = 0, is then equivalent to

    !               t       t
    !        r*z = q *b ,  p *d*p*z = 0 ,

    !  where x = p*z. if this system does not have full rank,
    !  then a least squares solution is obtained.  On output qrsolv
    !  also provides an upper triangular matrix s such that

    !         t   t               t
    !        p *(a *a + d*d)*p = s *s .

    !  s is computed within qrsolv and may be of separate interest.

    !  the subroutine statement is

    !    subroutine qrsolv(n, r, ldr, ipvt, diag, qtb, x, sdiag, wa)

    ! N.B. Arguments LDR and WA have been removed in this version.

    !  where

    !    n is a positive integer input variable set to the order of r.

    !    r is an n by n array.  On input the full upper triangle must contain
    !      the full upper triangle of the matrix r.
    !      On output the full upper triangle is unaltered, and the strict lower
    !      triangle contains the strict upper triangle (transposed) of the
    !      upper triangular matrix s.

    !    ldr is a positive integer input variable not less than n
    !      which specifies the leading dimension of the array r.

    !    ipvt is an integer input array of length n which defines the
    !      permutation matrix p such that a*p = q*r.  Column j of p
    !      is column ipvt(j) of the identity matrix.

    !    diag is an input array of length n which must contain the
    !      diagonal elements of the matrix d.

    !    qtb is an input array of length n which must contain the first
    !      n elements of the vector (q transpose)*b.

    !    x is an output array of length n which contains the least
    !      squares solution of the system a*x = b, d*x = 0.

    !    sdiag is an output array of length n which contains the
    !      diagonal elements of the upper triangular matrix s.

    !    wa is a work array of length n.

    !  subprograms called

    !    fortran-supplied ... ABS,SQRT

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

    !  **********
    INTEGER   :: i, j, k, kp1, l, nsing
    REAL (lmdp) :: COS, cotan, qtbpj, SIN, sum, TAN, temp, wa(n)
    REAL (lmdp), PARAMETER :: p5 = 0.5_lmdp, p25 = 0.25_lmdp, zero = 0.0_lmdp

    !     Copy r and (q transpose)*b to preserve input and initialize s.
    !     In particular, save the diagonal elements of r in x.

    DO  j = 1, n
       r(j:n,j) = r(j,j:n)
       x(j) = r(j,j)
       wa(j) = qtb(j)
    END DO

    !     Eliminate the diagonal matrix d using a givens rotation.

    DO  j = 1, n

       !        Prepare the row of d to be eliminated, locating the
       !        diagonal element using p from the qr factorization.

       l = ipvt(j)
       IF (diag(l) == zero) CYCLE
       sdiag(j:n) = zero
       sdiag(j) = diag(l)

       !     The transformations to eliminate the row of d modify only a single
       !     element of (q transpose)*b beyond the first n, which is initially zero.

       qtbpj = zero
       DO  k = j, n

          !        Determine a givens rotation which eliminates the
          !        appropriate element in the current row of d.

          IF (sdiag(k) == zero) CYCLE
          IF (ABS(r(k,k)) < ABS(sdiag(k))) THEN
             cotan = r(k,k)/sdiag(k)
             SIN = p5/SQRT(p25 + p25*cotan**2)
             COS = SIN*cotan
          ELSE
             TAN = sdiag(k)/r(k,k)
             COS = p5/SQRT(p25 + p25*TAN**2)
             SIN = COS*TAN
          END IF

          !        Compute the modified diagonal element of r and
          !        the modified element of ((q transpose)*b,0).

          r(k,k) = COS*r(k,k) + SIN*sdiag(k)
          temp = COS*wa(k) + SIN*qtbpj
          qtbpj = -SIN*wa(k) + COS*qtbpj
          wa(k) = temp

          !        Accumulate the tranformation in the row of s.

          kp1 = k + 1
          DO  i = kp1, n
             temp = COS*r(i,k) + SIN*sdiag(i)
             sdiag(i) = -SIN*r(i,k) + COS*sdiag(i)
             r(i,k) = temp
          END DO
       END DO

       !     Store the diagonal element of s and restore
       !     the corresponding diagonal element of r.

       sdiag(j) = r(j,j)
       r(j,j) = x(j)
    END DO

    !     Solve the triangular system for z.  If the system is singular,
    !     then obtain a least squares solution.

    nsing = n
    DO  j = 1, n
       IF (sdiag(j) == zero .AND. nsing == n) nsing = j - 1
       IF (nsing < n) wa(j) = zero
    END DO

    DO  k = 1, nsing
       j = nsing - k + 1
       sum = DOT_PRODUCT( r(j+1:nsing,j), wa(j+1:nsing) )
       wa(j) = (wa(j) - sum)/sdiag(j)
    END DO

    !     Permute the components of z back to components of x.

    DO  j = 1, n
       l = ipvt(j)
       x(l) = wa(j)
    END DO
    RETURN

    !     last card of subroutine qrsolv.

  END SUBROUTINE qrsolv



  FUNCTION enorm(n,x) RESULT(fn_val)

    INTEGER, INTENT(IN)    :: n
    REAL (lmdp), INTENT(IN)  :: x(:)
    REAL (lmdp)              :: fn_val

    !  **********

    !  function enorm

    !  given an n-vector x, this function calculates the euclidean norm of x.

    !  the euclidean norm is computed by accumulating the sum of squares in
    !  three different sums.  The sums of squares for the small and large
    !  components are scaled so that no overflows occur.  Non-destructive
    !  underflows are permitted.  Underflows and overflows do not occur in the
    !  computation of the unscaled sum of squares for the intermediate
    !  components.  The definitions of small, intermediate and large components
    !  depend on two constants, rdwarf and rgiant.  The main restrictions on
    !  these constants are that rdwarf**2 not underflow and rgiant**2 not
    !  overflow.  The constants given here are suitable for every known computer.

    !  the function statement is

    !    REAL (lmdp) function enorm(n,x)

    !  where

    !    n is a positive integer input variable.

    !    x is an input array of length n.

    !  subprograms called

    !    fortran-supplied ... ABS,SQRT

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

    !  **********
    INTEGER   :: i
    REAL (lmdp) :: agiant, floatn, s1, s2, s3, xabs, x1max, x3max
    REAL (lmdp), PARAMETER :: one = 1.0_lmdp, zero = 0.0_lmdp, rdwarf = 3.834E-20_lmdp,  &
         rgiant = 1.304E+19_lmdp

    s1 = zero
    s2 = zero
    s3 = zero
    x1max = zero
    x3max = zero
    floatn = n
    agiant = rgiant/floatn
    DO  i = 1, n
       xabs = ABS(x(i))
       IF (xabs > rdwarf .AND. xabs < agiant) GO TO 70
       IF (xabs <= rdwarf) GO TO 30

       !              sum for large components.

       IF (xabs <= x1max) GO TO 10
       s1 = one + s1*(x1max/xabs)**2
       x1max = xabs
       GO TO 20

10     s1 = s1 + (xabs/x1max)**2

20     GO TO 60

       !              sum for small components.

30     IF (xabs <= x3max) GO TO 40
       s3 = one + s3*(x3max/xabs)**2
       x3max = xabs
       GO TO 60

40     IF (xabs /= zero) s3 = s3 + (xabs/x3max)**2

60     CYCLE

       !           sum for intermediate components.

70     s2 = s2 + xabs**2
    END DO

    !     calculation of norm.

    IF (s1 == zero) GO TO 100
    fn_val = x1max*SQRT(s1 + (s2/x1max)/x1max)
    GO TO 120

100 IF (s2 == zero) GO TO 110
    IF (s2 >= x3max) fn_val = SQRT(s2*(one + (x3max/s2)*(x3max*s3)))
    IF (s2 < x3max) fn_val = SQRT(x3max*((s2/x3max) + (x3max*s3)))
    GO TO 120

110 fn_val = x3max*SQRT(s3)

120 RETURN

    !     last card of function enorm.

  END FUNCTION enorm



  SUBROUTINE fdjac2(fcn1, m, n, x, fvec, fjac, iflag, epsfcn)

    INTEGER, INTENT(IN)        :: m
    INTEGER, INTENT(IN)        :: n
    REAL (lmdp), INTENT(IN OUT)  :: x(n)
    REAL (lmdp), INTENT(IN)      :: fvec(m)
    REAL (lmdp), INTENT(OUT)     :: fjac(:,:)    ! fjac(ldfjac,n)
    INTEGER, INTENT(IN OUT)    :: iflag
    REAL (lmdp), INTENT(IN)      :: epsfcn

    include 'lm.h'

    !  **********

    !  subroutine fdjac2

    !  this subroutine computes a forward-difference approximation
    !  to the m by n jacobian matrix associated with a specified
    !  problem of m functions in n variables.

    !  the subroutine statement is

    !    subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)

    !  where

    !    fcn is the name of the user-supplied subroutine which calculates the
    !      functions.  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,iflag)
    !      integer m,n,iflag
    !      REAL (lmdp) x(n),fvec(m)
    !      ----------
    !      calculate the functions at x and
    !      return this vector in fvec.
    !      ----------
    !      return
    !      end

    !      the value of iflag should not be changed by fcn unless
    !      the user wants to terminate execution of fdjac2.
    !      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 input array of length n.

    !    fvec is an input array of length m which must contain the
    !      functions evaluated at x.

    !    fjac is an output m by n array which contains the
    !      approximation to the jacobian matrix evaluated at x.

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

    !    iflag is an integer variable which can be used to terminate
    !      the execution of fdjac2.  see description of fcn.

    !    epsfcn is an input variable used in determining a suitable step length
    !      for the forward-difference approximation.  This approximation assumes
    !      that the relative errors in the functions are of the order of epsfcn.
    !      If epsfcn is less than the machine precision, it is assumed that the
    !      relative errors in the functions are of the order of the machine
    !      precision.

    !    wa is a work array of length m.

    !  subprograms called

    !    user-supplied ...... fcn

    !    minpack-supplied ... dpmpar

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

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

    !  **********
    INTEGER   :: j
    REAL (lmdp) :: eps, epsmch, h, temp, wa(m)
    REAL (lmdp), PARAMETER :: zero = 0.0_lmdp

    !     epsmch is the machine precision.

    epsmch = EPSILON(zero)

    eps = SQRT(MAX(epsfcn, epsmch))
    DO  j = 1, n
       temp = x(j)
       h = eps*ABS(temp)
       IF (h == zero) h = eps
       x(j) = temp + h
       CALL fcn1(m, n, x, wa, iflag)
       IF (iflag < 0) EXIT
       x(j) = temp
       fjac(1:m,j) = (wa(1:m) - fvec(1:m))/h
    END DO

    RETURN

    !     last card of subroutine fdjac2.

  END SUBROUTINE fdjac2


end module lm