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