module fvn_misc
use fvn_common
implicit none

! Muller
interface fvn_muller
    module procedure fvn_z_muller
end interface fvn_muller

public
private :: lmdif,lmpar,qrfac,qrsolv,enorm,fdjac2
! These are made private to not interfere with a
! possibly linked minpack

contains
!
! Muller
!
!
!
! William Daniau 2007
!
! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f
! http://plato.asu.edu/ftp/other_software/muller.f
!
! it can be used as a replacement for imsl routine dzanly with minor changes
!
!-----------------------------------------------------------------------
!
!   purpose             - zeros of an analytic complex function
!                           using the muller method with deflation
!
!   usage               - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
!                           infer,ier)
!
!   arguments    f      - a complex function subprogram, f(z), written
!                           by the user specifying the equation whose
!                           roots are to be found.  f must appear in
!                           an external statement in the calling pro-
!                           gram.
!                eps    - 1st stopping criterion.  let fp(z)=f(z)/p
!                           where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))
!                           and z(1),...,z(k-1) are previously found
!                           roots.  if ((cdabs(f(z)).le.eps) .and.
!                           (cdabs(fp(z)).le.eps)), then z is accepted
!                           as a root. (input)
!                eps1   - 2nd stopping criterion.  a root is accepted
!                           if two successive approximations to a given
!                           root agree within eps1. (input)
!                             note. if either or both of the stopping
!                             criteria are fulfilled, the root is
!                             accepted.
!                kn     - the number of known roots which must be stored
!                           in x(1),...,x(kn), prior to entry to muller
!                nguess - the number of initial guesses provided. these
!                           guesses must be stored in x(kn+1),...,
!                           x(kn+nguess).  nguess must be set equal
!                           to zero if no guesses are provided. (input)
!                n      - the number of new roots to be found by
!                           muller (input)
!                x      - a complex vector of length kn+n.  x(1),...,
!                           x(kn) on input must contain any known
!                           roots.  x(kn+1),..., x(kn+n) on input may,
!                           on user option, contain initial guesses for
!                           the n new roots which are to be computed.
!                           if the user does not provide an initial
!                           guess, zero is used.
!                           on output, x(kn+1),...,x(kn+n) contain the
!                           approximate roots found by muller.
!                itmax  - the maximum allowable number of iterations
!                           per root (input)
!                infer  - an integer vector of length kn+n.  on
!                           output infer(j) contains the number of
!                           iterations used in finding the j-th root
!                           when convergence was achieved.  if
!                           convergence was not obtained in itmax
!                           iterations, infer(j) will be greater than
!                           itmax (output).
!                ier    - error parameter (output)
!                         warning error
!                           ier = 33 indicates failure to converge with-
!                             in itmax iterations for at least one of
!                             the (n) new roots.
!
!
!   remarks      muller always returns the last approximation for root j
!                in x(j). if the convergence criterion is satisfied,
!                then infer(j) is less than or equal to itmax. if the
!                convergence criterion is not satisified, then infer(j)
!                is set to either itmax+1 or itmax+k, with k greater
!                than 1. infer(j) = itmax+1 indicates that muller did
!                not obtain convergence in the allowed number of iter-
!                ations. in this case, the user may wish to set itmax
!                to a larger value. infer(j) = itmax+k means that con-
!                vergence was obtained (on iteration k) for the defla-
!                ted function
!                              fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
!
!                but failed for f(z). in this case, better initial
!                guesses might help or, it might be necessary to relax
!                the convergence criterion.
!
!-----------------------------------------------------------------------
!
subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
     implicit none
      double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq,eps1w
      double complex ::   d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
                          tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
                          zero,p1,one,four,p5

      double complex, external :: f
      integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
                    knpng,jk,ick,nn,lm1,errcode
      double complex :: x(kn+n)
      integer :: infer(kn+n)


      data                zero/(0.0d0,0.0d0)/,p1/(0.1d0,0.0d0)/, &
                          one/(1.0d0,0.0d0)/,four/(4.0d0,0.0d0)/, &
                          p5/(0.5d0,0.0d0)/, &
                          rzero/0.0d0/,rten/10.0d0/,rhun/100.0d0/, &
                          ax/0.1d0/,ickmax/3/,rp01/0.01d0/

            ier = 0
            if (n .lt. 1) then ! What the hell are doing here then ...
                return
            end if
            !eps1 = rten **(-nsig)
            eps1w = min(eps1,rp01)

            knp1 = kn+1
            knpn = kn+n
            knpng = kn+nguess
            do i=1,knpn
                infer(i) = 0
                if (i .gt. knpng) x(i) = zero
            end do
            l= knp1

            ic=0
rloop:      do while (l<=knpn)   ! Main loop over new roots
                jk = 0
                ick = 0
                xl = x(l)
icloop:         do
                    ic = 0
                    h = ax
                    h = p1*h
                    if (cdabs(xl) .gt. ax) h = p1*xl
!                                  first three points are
!                                    xl+h,  xl-h,  xl
                    rt = xl+h
                    call deflated_work(errcode)
                    if (errcode == 1) then
                        exit icloop
                    end if

                    z0 = fprt
                    y0 = frt
                    x0 = rt
                    rt = xl-h
                    call deflated_work(errcode)
                    if (errcode == 1) then
                        exit icloop
                    end if

                    z1 = fprt
                    y1 = frt
                    h = xl-rt
                    d = h/(rt-x0)
                    rt = xl

                    call deflated_work(errcode)
                    if (errcode == 1) then
                        exit icloop
                    end if

 
                    z2 = fprt
                    y2 = frt
!                                  begin main algorithm
 iloop:             do
                        dd = one + d
                        t1 = z0*d*d
                        t2 = z1*dd*dd
                        xx = z2*dd
                        t3 = z2*d
                        bi = t1-t2+xx+t3
                        den = bi*bi-four*(xx*t1-t3*(t2-xx))
!                                  use denominator of maximum amplitude 
                        t1 = sqrt(den)
                        qz = rhun*max(cdabs(bi),cdabs(t1))
                        t2 = bi + t1
                        tpq = cdabs(t2)+qz
                        if (tpq .eq. qz) t2 = zero
                        t3 = bi - t1
                        tpq = cdabs(t3) + qz
                        if (tpq .eq. qz) t3 = zero
                        den = t2
                        qz = cdabs(t3)-cdabs(t2)
                        if (qz .gt. rzero) den = t3
!                                  test for zero denominator            
                        if (cdabs(den) .eq. rzero) then
                            call trans_rt()
                            call deflated_work(errcode)
                            if (errcode == 1) then
                                exit icloop
                            end if
                            z2 = fprt
                            y2 = frt
                            cycle iloop
                        end if


                        d = -xx/den
                        d = d+d
                        h = d*h
                        rt = rt + h
!                                  check convergence of the first kind  
                        if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then
                            if (ic .ne. 0) then
                                exit icloop
                            end if
                            ic = 1
                            z0 = y1
                            z1 = y2
                            z2 = f(rt)
                            xl = rt
                            ick = ick+1
                            if (ick .le. ickmax) then
                                cycle iloop 
                            end if
!                                  warning error, itmax = maximum
                            jk = itmax + jk
                            ier = 33
                        end if
                        if (ic .ne. 0) then
                            cycle icloop
                        end if
                        call deflated_work(errcode)
                        if (errcode == 1) then
                            exit icloop
                        end if

                        do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
                            !   take remedial action to induce
                            !   convergence
                            d = d*p5
                            h = h*p5
                            rt = rt-h
                            call deflated_work(errcode)
                            if (errcode == 1) then
                                exit icloop
                            end if
                        end do
                        z0 = z1
                        z1 = z2
                        z2 = fprt
                        y0 = y1
                        y1 = y2
                        y2 = frt
                    end do iloop
                end do icloop
        x(l) = rt
        infer(l) = jk
        l = l+1
      end do rloop

      contains
        subroutine trans_rt()
           tem = rten*eps1w
           if (cdabs(rt) .gt. ax) tem = tem*rt
           rt = rt+tem
           d = (h+tem)*d/h
           h = h+tem
        end subroutine trans_rt
 
        subroutine deflated_work(errcode)
            ! errcode=0 => no errors
            ! errcode=1 => jk>itmax or convergence of second kind achieved
            integer :: errcode,flag
 
            flag=1
    loop1:  do while(flag==1)
                errcode=0
                jk = jk+1
                if (jk .gt. itmax) then
                    ier=33
                    errcode=1
                    return
                end if
                frt = f(rt)
                fprt = frt
                if (l /= 1) then
                    lm1 = l-1
                    do i=1,lm1
                        tem = rt - x(i)
                        if (cdabs(tem) .eq. rzero) then
                        !if (ic .ne. 0) go to 15 !! ?? possible?
                            call trans_rt()
                            cycle loop1
                        end if
                        fprt = fprt/tem
                    end do
                end if
                flag=0
            end do loop1
 
            if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
                errcode=1
                return
            end if

        end subroutine deflated_work

      end subroutine

!
!
!
! Non linear least square using Levenberg-Marquardt algorithm and
! a finite difference jacobian
!
! This uses MINPACK Routines (http://www.netlib.org/minpack)
! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
!
!  The purpose of fvn_lm 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.
!
!  call fvn_lm(fcn,m,n,a,info,tol)
!       
!       fcn : fcn is the user-supplied subroutine which calculates
!        the functions.  fcn must follow the following interface that must
!        be declared in the calling subroutine :
!
!         interface
!             subroutine fcn(m,n,a,fvec,iflag)
!                 use fvn_common
!                 integer(ip_kind), intent(in) :: m
!                 integer(ip_kind), intent(in) :: n
!                 real(dp_kind), dimension(:), intent(in) :: a
!                 real(dp_kind), dimension(:), intent(inout) :: fvec
!                 integer(ip_kind), intent(inout) :: iflag
!             end subroutine
!         end interface
!
!          This is the function which calculate the differences for which which square sum 
!          will be minimized outputing this difference in vector fvec.
!           Parameters of fcn are as follows :
!               m : positive integer input variable set to the number of functions
!               (number of measurement points)
!               n : positive integer input variable set to the number of variables
!               (number of parameters in the function to fit)
!               a : vector of length n containing parameters for which fcn should
!               perform the calculation
!               fvec : vector of length m containing the resulting evaluation
!               iflag : integer normally not used, can be used to exit the
!               the algorithm by setting it to a negative value
!  
!       m : positive integer input variable set to the number of functions
!               (number of measurement points)
!       n : positive integer input variable set to the number of variables
!               (number of parameters in the function to fit)
!       a : vector of length n, on input must contains an initial guess (or unity vector)
!           and on output the solution vector
!       info : is an output positive integer
!           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.
!       tol : is an optional positive value. 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. If not provided default value is :
!               sqrt(epsilon(0.0d0))
!
subroutine fvn_lm(fcn,m,n,a,info,tol)
    integer(ip_kind), intent(in) :: m
    integer(ip_kind), intent(in) :: n
    real(dp_kind), dimension(:), intent(inout) :: a
    integer(ip_kind), intent(out) :: info
    real(dp_kind), intent(in), optional :: tol

    real(dp_kind) :: rtol
    real(dp_kind), dimension(:), allocatable :: fvec
    integer(ip_kind), dimension(:), allocatable :: iwa

    interface
        subroutine fcn(m,n,a,fvec,iflag)
            use fvn_common
            integer(ip_kind), intent(in) :: m
            integer(ip_kind), intent(in) :: n
            real(dp_kind), dimension(:), intent(in) :: a
            real(dp_kind), dimension(:), intent(inout) :: fvec
            integer(ip_kind), intent(inout) :: iflag
        end subroutine
    end interface

    integer(ip_kind) :: maxfev, mode, nfev, nprint
    real(dp_kind) :: epsfcn, ftol, gtol, xtol, fjac(m,n)
    real(dp_kind), parameter :: factor = 100._8, zero = 0.0_8

    allocate(fvec(m),iwa(n))

    rtol=sqrt(epsilon(0.d0))
    if (present(tol)) rtol=tol

    info = 0

    !     check the input parameters for errors.

    if (n <= 0 .or. m < n .or. rtol < zero) return

    !     call lmdif.

    maxfev = 200*(n + 1)
    ftol = rtol
    xtol = rtol
    gtol = zero
    epsfcn = zero
    mode = 1
    nprint = 0

    call lmdif(fcn, m, n, a, fvec, ftol, xtol, gtol, maxfev, epsfcn,  &
           mode, factor, nprint, info, nfev, fjac, iwa)

    if (info == 8) info = 4

    deallocate(fvec,iwa)
end subroutine


! lmdif
! MINPACK Subroutine (http://www.netlib.org/minpack)
! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
SUBROUTINE lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn,  &
                 mode, factor, nprint, info, nfev, fjac, ipvt)
 
! N.B. Arguments LDFJAC, DIAG, QTF, WA1, WA2, WA3 & WA4 have been removed.
INTEGER, PARAMETER         :: dp = dp_kind
INTEGER, INTENT(IN)        :: m
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: x(:)
REAL (dp), INTENT(OUT)     :: fvec(:)
REAL (dp), INTENT(IN)      :: ftol
REAL (dp), INTENT(IN)      :: xtol
REAL (dp), INTENT(IN OUT)  :: gtol
INTEGER, INTENT(IN OUT)    :: maxfev
REAL (dp), INTENT(IN OUT)  :: epsfcn
INTEGER, INTENT(IN)        :: mode
REAL (dp), INTENT(IN)      :: factor
INTEGER, INTENT(IN)        :: nprint
INTEGER, INTENT(OUT)       :: info
INTEGER, INTENT(OUT)       :: nfev
REAL (dp), INTENT(OUT)     :: fjac(:,:)    ! fjac(ldfjac,n)
INTEGER, INTENT(OUT)       :: ipvt(:)

! EXTERNAL fcn

INTERFACE
  SUBROUTINE fcn(m, n, x, fvec, iflag)
    use fvn_common
    INTEGER(ip_kind), INTENT(IN)        :: m, n
    REAL (dp_kind), INTENT(IN)      :: x(:)
    REAL (dp_kind), INTENT(IN OUT)  :: fvec(:)
    INTEGER(ip_kind), INTENT(IN OUT)    :: iflag
  END SUBROUTINE fcn
END INTERFACE

!  **********

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

!     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, 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(fcn, 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 fcn(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 fcn(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 fcn(m, n, x, fvec, iflag)
RETURN

!     last card of subroutine lmdif.

END SUBROUTINE lmdif


!  **********

! lmpar
! MINPACK Subroutine (http://www.netlib.org/minpack)
! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
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, PARAMETER :: dp = dp_kind
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: r(:,:)
INTEGER, INTENT(IN)        :: ipvt(:)
REAL (dp), INTENT(IN)      :: diag(:)
REAL (dp), INTENT(IN)      :: qtb(:)
REAL (dp), INTENT(IN)      :: delta
REAL (dp), INTENT(OUT)     :: par
REAL (dp), INTENT(OUT)     :: x(:)
REAL (dp), 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 (dp) :: dxnorm, dwarf, fp, gnorm, parc, parl, paru, sum, temp
REAL (dp) :: wa1(n), wa2(n)
REAL (dp), PARAMETER :: p1 = 0.1_dp, p001 = 0.001_dp, zero = 0.0_dp

!     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


! qrfac
! MINPACK Subroutine (http://www.netlib.org/minpack)
! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
SUBROUTINE qrfac(m, n, a, pivot, ipvt, rdiag, acnorm)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-12-09  Time: 12:46:17

! N.B. Arguments LDA, LIPVT & WA have been removed.
INTEGER, PARAMETER :: dp = dp_kind
INTEGER, INTENT(IN)        :: m
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: a(:,:)
LOGICAL, INTENT(IN)        :: pivot
INTEGER, INTENT(OUT)       :: ipvt(:)
REAL (dp), INTENT(OUT)     :: rdiag(:)
REAL (dp), 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 (dp) :: ajnorm, epsmch, sum, temp, wa(n)
REAL (dp), PARAMETER :: one = 1.0_dp, p05 = 0.05_dp, zero = 0.0_dp

!     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


! qrsolv
! MINPACK Subroutine (http://www.netlib.org/minpack)
! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag)
 
! N.B. Arguments LDR & WA have been removed.
INTEGER, PARAMETER :: dp = dp_kind
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: r(:,:)
INTEGER, INTENT(IN)        :: ipvt(:)
REAL (dp), INTENT(IN)      :: diag(:)
REAL (dp), INTENT(IN)      :: qtb(:)
REAL (dp), INTENT(OUT)     :: x(:)
REAL (dp), 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 (dp) :: COS, cotan, qtbpj, SIN, sum, TAN, temp, wa(n)
REAL (dp), PARAMETER :: p5 = 0.5_dp, p25 = 0.25_dp, zero = 0.0_dp

!     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


! enorm
! MINPACK Subroutine (http://www.netlib.org/minpack)
! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
FUNCTION enorm(n,x) RESULT(fn_val)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-12-09  Time: 12:45:34
INTEGER, PARAMETER :: dp = dp_kind
INTEGER, INTENT(IN)    :: n
REAL (dp), INTENT(IN)  :: x(:)
REAL (dp)              :: 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 (dp) 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 (dp) :: agiant, floatn, s1, s2, s3, xabs, x1max, x3max
REAL (dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp, rdwarf = 3.834E-20_dp,  &
                        rgiant = 1.304E+19_dp

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


! fdjac2
! MINPACK Subroutine (http://www.netlib.org/minpack)
! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au
SUBROUTINE fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-12-09  Time: 12:45:44

! N.B. Arguments LDFJAC & WA have been removed.
INTEGER, PARAMETER :: dp = dp_kind
INTEGER, INTENT(IN)        :: m
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: x(n)
REAL (dp), INTENT(IN)      :: fvec(m)
REAL (dp), INTENT(OUT)     :: fjac(:,:)    ! fjac(ldfjac,n)
INTEGER, INTENT(IN OUT)    :: iflag
REAL (dp), INTENT(IN)      :: epsfcn

INTERFACE
  SUBROUTINE fcn(m, n, x, fvec, iflag)
    use fvn_common
    INTEGER(ip_kind), INTENT(IN)        :: m, n
    REAL (dp_kind), INTENT(IN)      :: x(:)
    REAL (dp_kind), INTENT(IN OUT)  :: fvec(:)
    INTEGER(ip_kind), INTENT(IN OUT)    :: iflag
  END SUBROUTINE fcn
END INTERFACE

!  **********

!  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 (dp) 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 (dp) :: eps, epsmch, h, temp, wa(m)
REAL (dp), PARAMETER :: zero = 0.0_dp

!     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 fcn(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 fvn_misc