module minpack

  use const_and_precisions, only : wp_
  implicit none

contains

  subroutine hybrj1(fcn,n,x,fvec,fjac,ldfjac,tol,info,wa,lwa)
    use const_and_precisions, only : zero, one
    implicit none
!   arguments
    integer, intent(in) :: n, ldfjac, lwa
    integer, intent(out) :: info
    real(wp_), intent(in) :: tol
    real(wp_), intent(out) :: wa(lwa)
    real(wp_), intent(inout) :: fvec(n), fjac(ldfjac,n), x(n)
!   **********
!  
!   subroutine hybrj1
!  
!   the purpose of hybrj1 is to find a zero of a system of
!   n nonlinear functions in n variables by a modification
!   of the powell hybrid method. this is done by using the
!   more general nonlinear equation solver hybrj. the user
!   must provide a subroutine which calculates the functions
!   and the jacobian.
!  
!   the subroutine statement is
!  
!     subroutine hybrj1(fcn,n,x,fvec,fjac,ldfjac,tol,info,wa,lwa)
!  
!   where
!  
!     fcn is the name of the user-supplied subroutine which
!       calculates the functions and the jacobian. fcn must
!       be declared in an external statement in the user
!       calling program, and should be written as follows.
!  
!       subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
!       integer n,ldfjac,iflag
!       real(8) x(n),fvec(n),fjac(ldfjac,n)
!       ----------
!       if iflag = 1 calculate the functions at x and
!       return this vector in fvec. do not alter fjac.
!       if iflag = 2 calculate the jacobian at x and
!       return this matrix in fjac. do not alter fvec.
!       ---------
!       return
!       end
!  
!       the value of iflag should not be changed by fcn unless
!       the user wants to terminate execution of hybrj1.
!       in this case set iflag to a negative integer.
!  
!     n is a positive integer input variable set to the number
!       of functions and variables.
!  
!     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 n which contains
!       the functions evaluated at the output x.
!  
!     fjac is an output n by n array which contains the
!       orthogonal matrix q produced by the qr factorization
!       of the final approximate jacobian.
!  
!     ldfjac is a positive integer input variable not less than n
!       which specifies the leading dimension of the array fjac.
!  
!     tol is a nonnegative input variable. termination occurs
!       when the algorithm estimates that the relative error
!       between x and the solution is at most tol.
!  
!     info is an integer output variable. if the user has
!       terminated execution, info is set to the (negative)
!       value of iflag. see description of fcn. otherwise,
!       info is set as follows.
!  
!       info = 0   improper input parameters.
!  
!       info = 1   algorithm estimates that the relative error
!                  between x and the solution is at most tol.
!  
!       info = 2   number of calls to fcn with iflag = 1 has
!                  reached 100*(n+1).
!  
!       info = 3   tol is too small. no further improvement in
!                  the approximate solution x is possible.
!  
!       info = 4   iteration is not making good progress.
!  
!     wa is a work array of length lwa.
!  
!     lwa is a positive integer input variable not less than
!       (n*(n+13))/2.
!  
!   subprograms called
!  
!     user-supplied ...... fcn
!  
!     minpack-supplied ... hybrj
!  
!   argonne national laboratory. minpack project. march 1980.
!   burton s. garbow, kenneth e. hillstrom, jorge j. more
!  
!   **********
!   local variables
    integer :: j, lr, maxfev, mode, nfev, njev, nprint
    real(wp_) :: xtol
!   parameters
    real(wp_), parameter :: factor=1.0e2_wp_
    
    interface
      subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
        use const_and_precisions, only : wp_
        implicit none
        integer, intent(in) :: n,ldfjac,iflag
        real(wp_), intent(in) :: x(n)
        real(wp_), intent(inout) :: fvec(n),fjac(ldfjac,n)
      end subroutine fcn
    end interface

    info = 0
!  
!   check the input parameters for errors.
!  
    if (n <= 0 .or. ldfjac < n .or. tol < zero &
        .or. lwa < (n*(n + 13))/2) return
!  
!   call hybrj.
!  
    maxfev = 100*(n + 1)
    xtol = tol
    mode = 2
    do j = 1, n
      wa(j) = one
    end do
    nprint = 0
    lr = (n*(n + 1))/2
    call hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,wa(1),mode, &
               factor,nprint,info,nfev,njev,wa(6*n+1),lr,wa(n+1), &
               wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
    if (info == 5) info = 4
  end subroutine hybrj1

  subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, &
                     factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2, &
                     wa3,wa4)
    use const_and_precisions, only : zero, one, epsmch=>comp_eps
    implicit none
!   arguments
    integer, intent(in) :: n, ldfjac, maxfev, mode, nprint, lr
    integer, intent(out) :: info, nfev, njev
    real(wp_), intent(in) :: xtol, factor
    real(wp_), intent(out) :: fvec(n), fjac(ldfjac,n), r(lr), qtf(n), &
                                    wa1(n), wa2(n), wa3(n), wa4(n)
    real(wp_), intent(inout) :: x(n), diag(n)
!   **********
!  
!   subroutine hybrj
!  
!   the purpose of hybrj is to find a zero of a system of
!   n nonlinear functions in n variables by a modification
!   of the powell hybrid method. the user must provide a
!   subroutine which calculates the functions and the jacobian.
!  
!   the subroutine statement is
!  
!     subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,
!                      mode,factor,nprint,info,nfev,njev,r,lr,qtf,
!                      wa1,wa2,wa3,wa4)
!  
!   where
!  
!     fcn is the name of the user-supplied subroutine which
!       calculates the functions and the jacobian. fcn must
!       be declared in an external statement in the user
!       calling program, and should be written as follows.
!  
!       subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
!       integer n,ldfjac,iflag
!       real(8) x(n),fvec(n),fjac(ldfjac,n)
!       ----------
!       if iflag = 1 calculate the functions at x and
!       return this vector in fvec. do not alter fjac.
!       if iflag = 2 calculate the jacobian at x and
!       return this matrix in fjac. do not alter fvec.
!       ---------
!       return
!       end
!  
!       the value of iflag should not be changed by fcn unless
!       the user wants to terminate execution of hybrj.
!       in this case set iflag to a negative integer.
!  
!     n is a positive integer input variable set to the number
!       of functions and variables.
!  
!     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 n which contains
!       the functions evaluated at the output x.
!  
!     fjac is an output n by n array which contains the
!       orthogonal matrix q produced by the qr factorization
!       of the final approximate jacobian.
!  
!     ldfjac is a positive integer input variable not less than n
!       which specifies the leading dimension of the array fjac.
!  
!     xtol is a nonnegative input variable. termination
!       occurs when the relative error between two consecutive
!       iterates is at most xtol.
!  
!     maxfev is a positive integer input variable. termination
!       occurs when the number of calls to fcn with iflag = 1
!       has reached maxfev.
!  
!     diag is an array of length n. if mode = 1 (see
!       below), diag is internally set. if mode = 2, diag
!       must contain positive entries that serve as
!       multiplicative scale factors for the variables.
!  
!     mode is an integer input variable. if mode = 1, the
!       variables will be scaled internally. if mode = 2,
!       the scaling is specified by the input diag. other
!       values of mode are equivalent to mode = 1.
!  
!     factor is a positive input variable used in determining the
!       initial step bound. this bound is set to the product of
!       factor and the euclidean norm of diag*x if nonzero, or else
!       to factor itself. in most cases factor should lie in the
!       interval (.1,100.). 100. is a generally recommended value.
!  
!     nprint is an integer input variable that enables controlled
!       printing of iterates if it is positive. in this case,
!       fcn is called with iflag = 0 at the beginning of the first
!       iteration and every nprint iterations thereafter and
!       immediately prior to return, with x and fvec available
!       for printing. fvec and fjac should not be altered.
!       if nprint is not positive, no special calls of fcn
!       with iflag = 0 are made.
!  
!     info is an integer output variable. if the user has
!       terminated execution, info is set to the (negative)
!       value of iflag. see description of fcn. otherwise,
!       info is set as follows.
!  
!       info = 0   improper input parameters.
!  
!       info = 1   relative error between two consecutive iterates
!                  is at most xtol.
!  
!       info = 2   number of calls to fcn with iflag = 1 has
!                  reached maxfev.
!  
!       info = 3   xtol is too small. no further improvement in
!                  the approximate solution x is possible.
!  
!       info = 4   iteration is not making good progress, as
!                  measured by the improvement from the last
!                  five jacobian evaluations.
!  
!       info = 5   iteration is not making good progress, as
!                  measured by the improvement from the last
!                  ten iterations.
!  
!     nfev is an integer output variable set to the number of
!       calls to fcn with iflag = 1.
!  
!     njev is an integer output variable set to the number of
!       calls to fcn with iflag = 2.
!  
!     r is an output array of length lr which contains the
!       upper triangular matrix produced by the qr factorization
!       of the final approximate jacobian, stored rowwise.
!  
!     lr is a positive integer input variable not less than
!       (n*(n+1))/2.
!  
!     qtf is an output array of length n which contains
!       the vector (q transpose)*fvec.
!  
!     wa1, wa2, wa3, and wa4 are work arrays of length n.
!  
!   subprograms called
!  
!     user-supplied ...... fcn
!  
!     minpack-supplied ... dogleg,enorm,
!                          qform,qrfac,r1mpyq,r1updt
!  
!     fortran-supplied ... abs,dmax1,dmin1,mod
!  
!   argonne national laboratory. minpack project. march 1980.
!   burton s. garbow, kenneth e. hillstrom, jorge j. more
!  
!   **********
!   local variables
    integer :: i, iflag, iter, j, jm1, l, ncfail, ncsuc, nslow1, nslow2
    integer, dimension(1) :: iwa
    logical :: jeval, sing
    real(wp_) :: actred, delta, fnorm, fnorm1, pnorm, prered, &
            ratio, summ, temp, xnorm
!   parameters
    real(wp_), parameter :: p1 = 1.0e-1_wp_, p5 = 5.0e-1_wp_,   &
               p001 = 1.0e-3_wp_, p0001 = 1.0e-4_wp_
    
    interface
      subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
        use const_and_precisions, only : wp_
        implicit none
        integer, intent(in) :: n,ldfjac,iflag
        real(wp_), intent(in) :: x(n)
        real(wp_), intent(inout) :: fvec(n),fjac(ldfjac,n)
      end subroutine fcn
    end interface
!  
    info = 0
    iflag = 0
    nfev = 0
    njev = 0
!  
!   check the input parameters for errors.
!  
    if (n <= 0 .or. ldfjac < n .or. xtol < zero &
        .or. maxfev <= 0 .or. factor <= zero &
        .or. lr < (n*(n + 1))/2) go to 300
    if (mode == 2) then
      do j = 1, n
        if (diag(j) <= zero) go to 300
      end do
    end if
!  
!   evaluate the function at the starting point
!   and calculate its norm.
!  
    iflag = 1
    call fcn(n,x,fvec,fjac,ldfjac,iflag)
    nfev = 1
    if (iflag < 0) go to 300
    fnorm = enorm(n,fvec)
!  
!   initialize iteration counter and monitors.
!  
    iter = 1
    ncsuc = 0
    ncfail = 0
    nslow1 = 0
    nslow2 = 0
!  
!   beginning of the outer loop.
!  
    do
      jeval = .true.
!  
!     calculate the jacobian matrix.
!  
      iflag = 2
      call fcn(n,x,fvec,fjac,ldfjac,iflag)
      njev = njev + 1
      if (iflag < 0) go to 300
!  
!     compute the qr factorization of the jacobian.
!  
      call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
!  
!     on the first iteration and if mode is 1, scale according
!     to the norms of the columns of the initial jacobian.
!  
      if (iter == 1) then
        if (mode /= 2) then
          do j = 1, n
            diag(j) = wa2(j)
            if (wa2(j) == zero) diag(j) = one
          end do
        end if
!  
!       on the first iteration, calculate the norm of the scaled x
!       and initialize the step bound delta.
!  
        do j = 1, n
          wa3(j) = diag(j)*x(j)
        end do
        xnorm = enorm(n,wa3)
        delta = factor*xnorm
        if (delta == zero) delta = factor
      end if
!  
!     form (q transpose)*fvec and store in qtf.
!  
      do i = 1, n
        qtf(i) = fvec(i)
      end do
      do j = 1, n
        if (fjac(j,j) /= zero) then
          summ = zero
          do i = j, n
            summ = summ + fjac(i,j)*qtf(i)
          end do
          temp = -summ/fjac(j,j)
          do i = j, n
            qtf(i) = qtf(i) + fjac(i,j)*temp
          end do
        end if
      end do
!  
!     copy the triangular factor of the qr factorization into r.
!  
      sing = .false.
      do j = 1, n
        l = j
        jm1 = j - 1
        do i = 1, jm1
          r(l) = fjac(i,j)
          l = l + n - i
        end do
        r(l) = wa1(j)
        if (wa1(j) == zero) sing = .true.
      end do
!  
!     accumulate the orthogonal factor in fjac.
!  
      call qform(n,n,fjac,ldfjac,wa1)
!  
!     rescale if necessary.
!  
      if (mode /= 2) then
        do j = 1, n
          diag(j) = dmax1(diag(j),wa2(j))
        end do
      end if
!  
!     beginning of the inner loop.
!  
      do
!  
!       if requested, call fcn to enable printing of iterates.
!  
        if (nprint > 0) then
          iflag = 0
          if (mod(iter-1,nprint) == 0) call fcn(n,x,fvec,fjac,ldfjac,iflag)
          if (iflag < 0) go to 300
        end if
!  
!       determine the direction p.
!  
        call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
!  
!       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 = dmin1(delta,pnorm)
!  
!       evaluate the function at x + p and calculate its norm.
!  
        iflag = 1
        call fcn(n,wa2,wa4,fjac,ldfjac,iflag)
        nfev = nfev + 1
        if (iflag < 0) go to 300
        fnorm1 = enorm(n,wa4)
!  
!       compute the scaled actual reduction.
!  
        actred = -one
        if (fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
!  
!       compute the scaled predicted reduction.
!  
        l = 1
        do i = 1, n
          summ = zero
          do j = i, n
            summ = summ + r(l)*wa1(j)
            l = l + 1
          end do
          wa3(i) = qtf(i) + summ
        end do
        temp = enorm(n,wa3)
        prered = zero
        if (temp < fnorm) prered = one - (temp/fnorm)**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 < p1) then
          ncsuc = 0
          ncfail = ncfail + 1
          delta = p5*delta
        else
          ncfail = 0
          ncsuc = ncsuc + 1
          if (ratio >= p5 .or. ncsuc > 1) delta = dmax1(delta,pnorm/p5)
          if (abs(ratio-one) <= p1) delta = pnorm/p5
        end if
!  
!       test for successful iteration.
!  
        if (ratio >= p0001) then
!  
!         successful iteration. update x, fvec, and their norms.
!  
          do j = 1, n
            x(j) = wa2(j)
            wa2(j) = diag(j)*x(j)
            fvec(j) = wa4(j)
          end do
          xnorm = enorm(n,wa2)
          fnorm = fnorm1
          iter = iter + 1
        end if
!  
!       determine the progress of the iteration.
!  
        nslow1 = nslow1 + 1
        if (actred >= p001) nslow1 = 0
        if (jeval) nslow2 = nslow2 + 1
        if (actred >= p1) nslow2 = 0
!  
!       test for convergence.
!  
        if (delta <= xtol*xnorm .or. fnorm == zero) info = 1
        if (info /= 0) go to 300
!  
!       tests for termination and stringent tolerances.
!  
        if (nfev >= maxfev) info = 2
        if (p1*dmax1(p1*delta,pnorm) <= epsmch*xnorm) info = 3
        if (nslow2 == 5) info = 4
        if (nslow1 == 10) info = 5
        if (info /= 0) go to 300
!  
!       criterion for recalculating jacobian.
!  
        if (ncfail == 2) exit
!  
!       calculate the rank one modification to the jacobian
!       and update qtf if necessary.
!  
        do j = 1, n
          summ = zero
          do i = 1, n
            summ = summ + fjac(i,j)*wa4(i)
          end do
          wa2(j) = (summ - wa3(j))/pnorm
          wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
          if (ratio >= p0001) qtf(j) = summ
        end do
!  
!       compute the qr factorization of the updated jacobian.
!  
        call r1updt(n,n,r,lr,wa1,wa2,wa3,sing)
        call r1mpyq(n,n,fjac,ldfjac,wa2,wa3)
        call r1mpyq(1,n,qtf,1,wa2,wa3)
!  
!       end of the inner loop.
!  
        jeval = .false.
      end do
!  
!      end of the outer loop.
!  
    end do
    300 continue
!  
!   termination, either normal or user imposed.
!  
    if (iflag < 0) info = iflag
    iflag = 0
    if (nprint > 0) call fcn(n,x,fvec,fjac,ldfjac,iflag)
  end subroutine hybrj

  subroutine hybrj1mv(fcn,n,x,f0,fvec,fjac,ldfjac,tol,info,wa,lwa)
    use const_and_precisions, only : zero, one
    implicit none
!   arguments
    integer, intent(in) :: n, ldfjac, lwa
    integer, intent(out) :: info
    real(wp_), intent(in) :: tol,f0(n)
    real(wp_), intent(out) :: wa(lwa)
    real(wp_), intent(inout) :: fvec(n), fjac(ldfjac,n), x(n)
!   **********
!  
!   subroutine hybrj1mv
!  
!   the purpose of hybrj1mv is to find a zero of a system of
!   n nonlinear functions in n variables by a modification
!   of the powell hybrid method. this is done by using the
!   more general nonlinear equation solver hybrjmv. the user
!   must provide a subroutine which calculates the functions
!   and the jacobian.
!  
!   the subroutine statement is
!  
!     subroutine hybrj1mv(fcn,n,x,f0,fvec,fjac,ldfjac,tol,info,wa,lwa)
!  
!   where
!  
!     fcn is the name of the user-supplied subroutine which
!       calculates the functions and the jacobian. fcn must
!       be declared in an external statement in the user
!       calling program, and should be written as follows.
!  
!       subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
!       integer n,ldfjac,iflag
!       real(8) x(n),fvec(n),fjac(ldfjac,n)
!       ----------
!       if iflag = 1 calculate the functions at x and
!       return this vector in fvec. do not alter fjac.
!       if iflag = 2 calculate the jacobian at x and
!       return this matrix in fjac. do not alter fvec.
!       ---------
!       return
!       end
!  
!       the value of iflag should not be changed by fcn unless
!       the user wants to terminate execution of hybrj1mv.
!       in this case set iflag to a negative integer.
!  
!     n is a positive integer input variable set to the number
!       of functions and variables.
!  
!     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 n which contains
!       the functions evaluated at the output x.
!  
!     fjac is an output n by n array which contains the
!       orthogonal matrix q produced by the qr factorization
!       of the final approximate jacobian.
!  
!     ldfjac is a positive integer input variable not less than n
!       which specifies the leading dimension of the array fjac.
!  
!     tol is a nonnegative input variable. termination occurs
!       when the algorithm estimates that the relative error
!       between x and the solution is at most tol.
!  
!     info is an integer output variable. if the user has
!       terminated execution, info is set to the (negative)
!       value of iflag. see description of fcn. otherwise,
!       info is set as follows.
!  
!       info = 0   improper input parameters.
!  
!       info = 1   algorithm estimates that the relative error
!                  between x and the solution is at most tol.
!  
!       info = 2   number of calls to fcn with iflag = 1 has
!                  reached 100*(n+1).
!  
!       info = 3   tol is too small. no further improvement in
!                  the approximate solution x is possible.
!  
!       info = 4   iteration is not making good progress.
!  
!     wa is a work array of length lwa.
!  
!     lwa is a positive integer input variable not less than
!       (n*(n+13))/2.
!  
!   subprograms called
!  
!     user-supplied ...... fcn
!  
!     minpack-supplied ... hybrjmv
!  
!   argonne national laboratory. minpack project. march 1980.
!   burton s. garbow, kenneth e. hillstrom, jorge j. more
!  
!   **********
!   local variables
    integer :: j, lr, maxfev, mode, nfev, njev, nprint
    real(wp_) :: xtol
!   parameters
    real(wp_), parameter :: factor=1.0e2_wp_
    
    interface
      subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
        use const_and_precisions, only : wp_
        implicit none
        integer, intent(in) :: n,ldfjac,iflag
        real(wp_), intent(in) :: x(n),f0(n)
        real(wp_), intent(inout) :: fvec(n),fjac(ldfjac,n)
      end subroutine fcn
    end interface

    info = 0
!  
!   check the input parameters for errors.
!  
    if (n <= 0 .or. ldfjac < n .or. tol < zero &
        .or. lwa < (n*(n + 13))/2) return
!  
!   call hybrjmv.
!  
    maxfev = 100*(n + 1)
    xtol = tol
    mode = 2
    do j = 1, n
      wa(j) = one
    end do
    nprint = 0
    lr = (n*(n + 1))/2
    call hybrjmv(fcn,n,x,f0,fvec,fjac,ldfjac,xtol,maxfev,wa(1),mode, &
               factor,nprint,info,nfev,njev,wa(6*n+1),lr,wa(n+1), &
               wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
    if (info == 5) info = 4
  end subroutine hybrj1mv

  subroutine hybrjmv(fcn,n,x,f0,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, &
                     factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2, &
                     wa3,wa4)
    use const_and_precisions, only : zero, one, epsmch=>comp_eps
    implicit none
!   arguments
    integer, intent(in) :: n, ldfjac, maxfev, mode, nprint, lr
    integer, intent(out) :: info, nfev, njev
    real(wp_), intent(in) :: xtol, factor, f0(n)
    real(wp_), intent(out) :: fvec(n), fjac(ldfjac,n), r(lr), qtf(n), &
                                    wa1(n), wa2(n), wa3(n), wa4(n)
    real(wp_), intent(inout) :: x(n), diag(n)
!   **********
!  
!   subroutine hybrj
!  
!   the purpose of hybrj is to find a zero of a system of
!   n nonlinear functions in n variables by a modification
!   of the powell hybrid method. the user must provide a
!   subroutine which calculates the functions and the jacobian.
!  
!   the subroutine statement is
!  
!     subroutine hybrj(fcn,n,x,f0,fvec,fjac,ldfjac,xtol,maxfev,diag,
!                      mode,factor,nprint,info,nfev,njev,r,lr,qtf,
!                      wa1,wa2,wa3,wa4)
!  
!   where
!  
!     fcn is the name of the user-supplied subroutine which
!       calculates the functions and the jacobian. fcn must
!       be declared in an external statement in the user
!       calling program, and should be written as follows.
!  
!       subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
!       integer n,ldfjac,iflag
!       real(8) x(n),f0(n),fvec(n),fjac(ldfjac,n)
!       ----------
!       if iflag = 1 calculate the functions at x and
!       return this vector in fvec. do not alter fjac.
!       if iflag = 2 calculate the jacobian at x and
!       return this matrix in fjac. do not alter fvec.
!       ---------
!       return
!       end
!  
!       the value of iflag should not be changed by fcn unless
!       the user wants to terminate execution of hybrj.
!       in this case set iflag to a negative integer.
!  
!     n is a positive integer input variable set to the number
!       of functions and variables.
!  
!     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 n which contains
!       the functions evaluated at the output x.
!  
!     fjac is an output n by n array which contains the
!       orthogonal matrix q produced by the qr factorization
!       of the final approximate jacobian.
!  
!     ldfjac is a positive integer input variable not less than n
!       which specifies the leading dimension of the array fjac.
!  
!     xtol is a nonnegative input variable. termination
!       occurs when the relative error between two consecutive
!       iterates is at most xtol.
!  
!     maxfev is a positive integer input variable. termination
!       occurs when the number of calls to fcn with iflag = 1
!       has reached maxfev.
!  
!     diag is an array of length n. if mode = 1 (see
!       below), diag is internally set. if mode = 2, diag
!       must contain positive entries that serve as
!       multiplicative scale factors for the variables.
!  
!     mode is an integer input variable. if mode = 1, the
!       variables will be scaled internally. if mode = 2,
!       the scaling is specified by the input diag. other
!       values of mode are equivalent to mode = 1.
!  
!     factor is a positive input variable used in determining the
!       initial step bound. this bound is set to the product of
!       factor and the euclidean norm of diag*x if nonzero, or else
!       to factor itself. in most cases factor should lie in the
!       interval (.1,100.). 100. is a generally recommended value.
!  
!     nprint is an integer input variable that enables controlled
!       printing of iterates if it is positive. in this case,
!       fcn is called with iflag = 0 at the beginning of the first
!       iteration and every nprint iterations thereafter and
!       immediately prior to return, with x and fvec available
!       for printing. fvec and fjac should not be altered.
!       if nprint is not positive, no special calls of fcn
!       with iflag = 0 are made.
!  
!     info is an integer output variable. if the user has
!       terminated execution, info is set to the (negative)
!       value of iflag. see description of fcn. otherwise,
!       info is set as follows.
!  
!       info = 0   improper input parameters.
!  
!       info = 1   relative error between two consecutive iterates
!                  is at most xtol.
!  
!       info = 2   number of calls to fcn with iflag = 1 has
!                  reached maxfev.
!  
!       info = 3   xtol is too small. no further improvement in
!                  the approximate solution x is possible.
!  
!       info = 4   iteration is not making good progress, as
!                  measured by the improvement from the last
!                  five jacobian evaluations.
!  
!       info = 5   iteration is not making good progress, as
!                  measured by the improvement from the last
!                  ten iterations.
!  
!     nfev is an integer output variable set to the number of
!       calls to fcn with iflag = 1.
!  
!     njev is an integer output variable set to the number of
!       calls to fcn with iflag = 2.
!  
!     r is an output array of length lr which contains the
!       upper triangular matrix produced by the qr factorization
!       of the final approximate jacobian, stored rowwise.
!  
!     lr is a positive integer input variable not less than
!       (n*(n+1))/2.
!  
!     qtf is an output array of length n which contains
!       the vector (q transpose)*fvec.
!  
!     wa1, wa2, wa3, and wa4 are work arrays of length n.
!  
!   subprograms called
!  
!     user-supplied ...... fcn
!  
!     minpack-supplied ... dogleg,enorm,
!                          qform,qrfac,r1mpyq,r1updt
!  
!     fortran-supplied ... abs,dmax1,dmin1,mod
!  
!   argonne national laboratory. minpack project. march 1980.
!   burton s. garbow, kenneth e. hillstrom, jorge j. more
!  
!   **********
!   local variables
    integer :: i, iflag, iter, j, jm1, l, ncfail, ncsuc, nslow1, nslow2
    integer, dimension(1) :: iwa
    logical :: jeval, sing
    real(wp_) :: actred, delta, fnorm, fnorm1, pnorm, prered, &
            ratio, summ, temp, xnorm
!   parameters
    real(wp_), parameter :: p1 = 1.0e-1_wp_, p5 = 5.0e-1_wp_,   &
               p001 = 1.0e-3_wp_, p0001 = 1.0e-4_wp_
    
    interface
      subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
        use const_and_precisions, only : wp_
        implicit none
        integer, intent(in) :: n,ldfjac,iflag
        real(wp_), intent(in) :: x(n),f0(n)
        real(wp_), intent(inout) :: fvec(n),fjac(ldfjac,n)
      end subroutine fcn
    end interface
!  
    info = 0
    iflag = 0
    nfev = 0
    njev = 0
!  
!   check the input parameters for errors.
!  
    if (n <= 0 .or. ldfjac < n .or. xtol < zero &
        .or. maxfev <= 0 .or. factor <= zero &
        .or. lr < (n*(n + 1))/2) go to 300
    if (mode == 2) then
      do j = 1, n
        if (diag(j) <= zero) go to 300
      end do
    end if
!  
!   evaluate the function at the starting point
!   and calculate its norm.
!  
    iflag = 1
    call fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
    nfev = 1
    if (iflag < 0) go to 300
    fnorm = enorm(n,fvec)
!  
!   initialize iteration counter and monitors.
!  
    iter = 1
    ncsuc = 0
    ncfail = 0
    nslow1 = 0
    nslow2 = 0
!  
!   beginning of the outer loop.
!  
    do
      jeval = .true.
!  
!     calculate the jacobian matrix.
!  
      iflag = 2
      call fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
      njev = njev + 1
      if (iflag < 0) go to 300
!  
!     compute the qr factorization of the jacobian.
!  
      call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
!  
!     on the first iteration and if mode is 1, scale according
!     to the norms of the columns of the initial jacobian.
!  
      if (iter == 1) then
        if (mode /= 2) then
          do j = 1, n
            diag(j) = wa2(j)
            if (wa2(j) == zero) diag(j) = one
          end do
        end if
!  
!       on the first iteration, calculate the norm of the scaled x
!       and initialize the step bound delta.
!  
        do j = 1, n
          wa3(j) = diag(j)*x(j)
        end do
        xnorm = enorm(n,wa3)
        delta = factor*xnorm
        if (delta == zero) delta = factor
      end if
!  
!     form (q transpose)*fvec and store in qtf.
!  
      do i = 1, n
        qtf(i) = fvec(i)
      end do
      do j = 1, n
        if (fjac(j,j) /= zero) then
          summ = zero
          do i = j, n
            summ = summ + fjac(i,j)*qtf(i)
          end do
          temp = -summ/fjac(j,j)
          do i = j, n
            qtf(i) = qtf(i) + fjac(i,j)*temp
          end do
        end if
      end do
!  
!     copy the triangular factor of the qr factorization into r.
!  
      sing = .false.
      do j = 1, n
        l = j
        jm1 = j - 1
        do i = 1, jm1
          r(l) = fjac(i,j)
          l = l + n - i
        end do
        r(l) = wa1(j)
        if (wa1(j) == zero) sing = .true.
      end do
!  
!     accumulate the orthogonal factor in fjac.
!  
      call qform(n,n,fjac,ldfjac,wa1)
!  
!     rescale if necessary.
!  
      if (mode /= 2) then
        do j = 1, n
          diag(j) = dmax1(diag(j),wa2(j))
        end do
      end if
!  
!     beginning of the inner loop.
!  
      do
!  
!       if requested, call fcn to enable printing of iterates.
!  
        if (nprint > 0) then
          iflag = 0
          if (mod(iter-1,nprint) == 0) call fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
          if (iflag < 0) go to 300
        end if
!  
!       determine the direction p.
!  
        call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
!  
!       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 = dmin1(delta,pnorm)
!  
!       evaluate the function at x + p and calculate its norm.
!  
        iflag = 1
        call fcn(n,wa2,f0,wa4,fjac,ldfjac,iflag)
        nfev = nfev + 1
        if (iflag < 0) go to 300
        fnorm1 = enorm(n,wa4)
!  
!       compute the scaled actual reduction.
!  
        actred = -one
        if (fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
!  
!       compute the scaled predicted reduction.
!  
        l = 1
        do i = 1, n
          summ = zero
          do j = i, n
            summ = summ + r(l)*wa1(j)
            l = l + 1
          end do
          wa3(i) = qtf(i) + summ
        end do
        temp = enorm(n,wa3)
        prered = zero
        if (temp < fnorm) prered = one - (temp/fnorm)**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 < p1) then
          ncsuc = 0
          ncfail = ncfail + 1
          delta = p5*delta
        else
          ncfail = 0
          ncsuc = ncsuc + 1
          if (ratio >= p5 .or. ncsuc > 1) delta = dmax1(delta,pnorm/p5)
          if (abs(ratio-one) <= p1) delta = pnorm/p5
        end if
!  
!       test for successful iteration.
!  
        if (ratio >= p0001) then
!  
!         successful iteration. update x, fvec, and their norms.
!  
          do j = 1, n
            x(j) = wa2(j)
            wa2(j) = diag(j)*x(j)
            fvec(j) = wa4(j)
          end do
          xnorm = enorm(n,wa2)
          fnorm = fnorm1
          iter = iter + 1
        end if
!  
!       determine the progress of the iteration.
!  
        nslow1 = nslow1 + 1
        if (actred >= p001) nslow1 = 0
        if (jeval) nslow2 = nslow2 + 1
        if (actred >= p1) nslow2 = 0
!  
!       test for convergence.
!  
        if (delta <= xtol*xnorm .or. fnorm == zero) info = 1
        if (info /= 0) go to 300
!  
!       tests for termination and stringent tolerances.
!  
        if (nfev >= maxfev) info = 2
        if (p1*dmax1(p1*delta,pnorm) <= epsmch*xnorm) info = 3
        if (nslow2 == 5) info = 4
        if (nslow1 == 10) info = 5
        if (info /= 0) go to 300
!  
!       criterion for recalculating jacobian.
!  
        if (ncfail == 2) exit
!  
!       calculate the rank one modification to the jacobian
!       and update qtf if necessary.
!  
        do j = 1, n
          summ = zero
          do i = 1, n
            summ = summ + fjac(i,j)*wa4(i)
          end do
          wa2(j) = (summ - wa3(j))/pnorm
          wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
          if (ratio >= p0001) qtf(j) = summ
        end do
!  
!       compute the qr factorization of the updated jacobian.
!  
        call r1updt(n,n,r,lr,wa1,wa2,wa3,sing)
        call r1mpyq(n,n,fjac,ldfjac,wa2,wa3)
        call r1mpyq(1,n,qtf,1,wa2,wa3)
!  
!       end of the inner loop.
!  
        jeval = .false.
      end do
!  
!      end of the outer loop.
!  
    end do
    300 continue
!  
!   termination, either normal or user imposed.
!  
    if (iflag < 0) info = iflag
    iflag = 0
    if (nprint > 0) call fcn(n,x,f0,fvec,fjac,ldfjac,iflag)
  end subroutine hybrjmv

  subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
    use const_and_precisions, only : zero, one, epsmch=>comp_eps
    implicit none
!   arguments
    integer, intent(in) :: n, lr
    real(wp_), intent(in) :: delta, r(lr), diag(n), qtb(n)
    real(wp_), intent(out) :: x(n), wa1(n), wa2(n)
!   **********
!  
!   subroutine dogleg
!  
!   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 the convex combination x of the
!   gauss-newton and scaled gradient directions that minimizes
!   (a*x - b) in the least squares sense, subject to the
!   restriction that the euclidean norm of d*x be at most delta.
!  
!   this subroutine completes the solution of the problem
!   if it is provided with the necessary information from the
!   qr factorization of a. that is, if a = q*r, where q has
!   orthogonal columns and r is an upper triangular matrix,
!   then dogleg expects the full upper triangle of r and
!   the first n components of (q transpose)*b.
!  
!   the subroutine statement is
!  
!     subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
!  
!   where
!  
!     n is a positive integer input variable set to the order of r.
!  
!     r is an input array of length lr which must contain the upper
!       triangular matrix r stored by rows.
!  
!     lr is a positive integer input variable not less than
!       (n*(n+1))/2.
!  
!     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.
!  
!     x is an output array of length n which contains the desired
!       convex combination of the gauss-newton direction and the
!       scaled gradient direction.
!  
!     wa1 and wa2 are work arrays of length n.
!  
!   subprograms called
!  
!     minpack-supplied ... enorm
!  
!     fortran-supplied ... abs,dmax1,dmin1,sqrt
!  
!   argonne national laboratory. minpack project. march 1980.
!   burton s. garbow, kenneth e. hillstrom, jorge j. more
!  
!   **********
!   local variables
    integer       :: i, j, jj, jp1, k, l
    real(wp_) :: alpha, bnorm, gnorm, qnorm, sgnorm, summ, temp
!  
!   first, calculate the gauss-newton direction.
!  
    jj = (n*(n + 1))/2 + 1
    do k = 1, n
      j = n - k + 1
      jp1 = j + 1
      jj = jj - k
      l = jj + 1
      summ = zero
      do i = jp1, n
        summ = summ + r(l)*x(i)
        l = l + 1
      end do
      temp = r(jj)
      if (temp == zero) then
        l = j
        do i = 1, j
          temp = dmax1(temp,abs(r(l)))
          l = l + n - i
        end do
        temp = epsmch*temp
        if (temp == zero) temp = epsmch
      end if
      x(j) = (qtb(j) - summ)/temp
    end do
!  
!   test whether the gauss-newton direction is acceptable.
!  
    do j = 1, n
      wa1(j) = zero
      wa2(j) = diag(j)*x(j)
    end do
    qnorm = enorm(n,wa2)
    if (qnorm <= delta) return
!  
!   the gauss-newton direction is not acceptable.
!   next, calculate the scaled gradient direction.
!  
    l = 1
    do j = 1, n
      temp = qtb(j)
      do i = j, n
        wa1(i) = wa1(i) + r(l)*temp
        l = l + 1
      end do
      wa1(j) = wa1(j)/diag(j)
    end do
!  
!   calculate the norm of the scaled gradient and test for
!   the special case in which the scaled gradient is zero.
!  
    gnorm = enorm(n,wa1)
    sgnorm = zero
    alpha = delta/qnorm
    if (gnorm /= zero) then
!  
!   calculate the point along the scaled gradient
!   at which the quadratic is minimized.
!  
      do j = 1, n
        wa1(j) = (wa1(j)/gnorm)/diag(j)
      end do
      l = 1
      do j = 1, n
        summ = zero
        do i = j, n
          summ = summ + r(l)*wa1(i)
          l = l + 1
        end do
        wa2(j) = summ
      end do
      temp = enorm(n,wa2)
      sgnorm = (gnorm/temp)/temp
!  
!   test whether the scaled gradient direction is acceptable.
!  
      alpha = zero
      if (sgnorm < delta) then
!  
!   the scaled gradient direction is not acceptable.
!   finally, calculate the point along the dogleg
!   at which the quadratic is minimized.
!  
        bnorm = enorm(n,qtb)
        temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
        temp = temp - (delta/qnorm)*(sgnorm/delta)**2 &
               + sqrt((temp-(delta/qnorm))**2 &
                       +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2))
        alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp
      end if
    end if
!  
!   form appropriate convex combination of the gauss-newton
!   direction and the scaled gradient direction.
!  
    temp = (one - alpha)*dmin1(sgnorm,delta)
    do j = 1, n
      x(j) = temp*wa1(j) + alpha*x(j)
    end do
  end subroutine dogleg

  function enorm(n,x)
    use const_and_precisions, only : zero, one
    implicit none
    real(wp_) :: enorm
    integer, intent(in) :: n
    real(wp_), dimension(n), intent(in) :: x
!   **********
!  
!   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(8) 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(wp_) :: agiant,floatn,s1,s2,s3,xabs,x1max,x3max
    real(wp_), parameter :: rdwarf=3.834e-20_wp_,rgiant=1.304e19_wp_
    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 .or. xabs >= agiant) then
        if (xabs > rdwarf) then
!  
!         sum for large components.
!  
          if (xabs > x1max) then
            s1 = one + s1*(x1max/xabs)**2
            x1max = xabs
          else
            s1 = s1 + (xabs/x1max)**2
          end if
        else
!  
!         sum for small components.
!  
          if (xabs > x3max) then
            s3 = one + s3*(x3max/xabs)**2
            x3max = xabs
          else
            if (xabs /= zero) s3 = s3 + (xabs/x3max)**2
          end if
        end if
      else
!  
!       sum for intermediate components.
!  
        s2 = s2 + xabs**2
      end if
    end do
!  
!   calculation of norm.
!  
    if (s1 /= zero) then
      enorm = x1max*sqrt(s1+(s2/x1max)/x1max)
    else
      if (s2 /= zero) then
        if (s2 >= x3max) enorm = sqrt(s2*(one+(x3max/s2)*(x3max*s3)))
        if (s2 < x3max) enorm = sqrt(x3max*((s2/x3max)+(x3max*s3)))
      else
        enorm = x3max*sqrt(s3)
      end if
    end if
  end function enorm

  subroutine qform(m,n,q,ldq,wa)
    use const_and_precisions, only : zero, one
    implicit none
!   arguments
    integer, intent(in) :: m,n,ldq
    real(wp_), intent(out) :: wa(m)
    real(wp_), intent(inout) :: q(ldq,m)
!   **********
!  
!   subroutine qform
!  
!   this subroutine proceeds from the computed qr factorization of
!   an m by n matrix a to accumulate the m by m orthogonal matrix
!   q from its factored form.
!  
!   the subroutine statement is
!  
!     subroutine qform(m,n,q,ldq,wa)
!  
!   where
!  
!     m is a positive integer input variable set to the number
!       of rows of a and the order of q.
!  
!     n is a positive integer input variable set to the number
!       of columns of a.
!  
!     q is an m by m array. on input the full lower trapezoid in
!       the first min(m,n) columns of q contains the factored form.
!       on output q has been accumulated into a square matrix.
!  
!     ldq is a positive integer input variable not less than m
!       which specifies the leading dimension of the array q.
!  
!     wa is a work array of length m.
!  
!   subprograms called
!  
!     fortran-supplied ... min0
!  
!   argonne national laboratory. minpack project. march 1980.
!   burton s. garbow, kenneth e. hillstrom, jorge j. more
!  
!   **********
!   local variables
    integer :: i, j, jm1, k, l, minmn, np1
    real(wp_) :: summ, temp
!  
!   zero out upper triangle of q in the first min(m,n) columns.
!  
    minmn = min0(m,n)
    do j = 2, minmn
      jm1 = j - 1
      do i = 1, jm1
        q(i,j) = zero
      end do
    end do
!  
!   initialize remaining columns to those of the identity matrix.
!  
    np1 = n + 1
    do j = np1, m
      do i = 1, m
        q(i,j) = zero
      end do
      q(j,j) = one
    end do
!  
!   accumulate q from its factored form.
!  
    do l = 1, minmn
      k = minmn - l + 1
      do i = k, m
        wa(i) = q(i,k)
        q(i,k) = zero
      end do
      q(k,k) = one
      if (wa(k) /= zero) then
        do j = k, m
          summ = zero
          do i = k, m
            summ = summ + q(i,j)*wa(i)
          end do
          temp = summ/wa(k)
          do i = k, m
            q(i,j) = q(i,j) - temp*wa(i)
          end do
        end do
      end if
    end do
  end subroutine qform

  subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
    use const_and_precisions, only : zero, one, epsmch=>comp_eps
    implicit none
!   arguments
    integer, intent(in) :: m, n, lda, lipvt
    integer, intent(out) :: ipvt(lipvt)
    logical, intent(in) :: pivot
    real(wp_), intent(out) :: rdiag(n), acnorm(n), wa(n)
    real(wp_), intent(inout) :: a(lda,n)
!   **********
!  
!   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)
!  
!   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 ... enorm
!  
!     fortran-supplied ... dmax1,sqrt,min0
!  
!   argonne national laboratory. minpack project. march 1980.
!   burton s. garbow, kenneth e. hillstrom, jorge j. more
!  
!   **********
!   local variables
    integer       :: i, j, jp1, k, kmax, minmn
    real(wp_) :: ajnorm, summ, temp
!   parameters
    real(wp_), parameter :: p05=5.0e-2_wp_
!  
!   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 = min0(m,n)
    do j = 1, minmn
      if (pivot) then
!  
!       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) then
          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
        end if
      end if
!  
!     compute the householder transformation to reduce the
!     j-th column of a to a multiple of the j-th unit vector.
!  
      ajnorm = enorm(m-j+1,a(j,j))
      if (ajnorm /= zero) then
        if (a(j,j) < zero) ajnorm = -ajnorm
        do i = j, m
          a(i,j) = a(i,j)/ajnorm
        end do
        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
          summ = zero
          do i = j, m
            summ = summ + a(i,j)*a(i,k)
          end do
          temp = summ/a(j,j)
          do i = j, m
            a(i,k) = a(i,k) - temp*a(i,j)
          end do
          if (pivot .and. rdiag(k) /= zero) then
            temp = a(j,k)/rdiag(k)
            rdiag(k) = rdiag(k)*sqrt(dmax1(zero,one-temp**2))
            if (p05*(rdiag(k)/wa(k))**2 <= epsmch) then
              rdiag(k) = enorm(m-j,a(jp1,k))
              wa(k) = rdiag(k)
            end if
          end if
        end do
      end if
      rdiag(j) = -ajnorm
    end do
  end subroutine qrfac

  subroutine r1mpyq(m,n,a,lda,v,w)
    use const_and_precisions, only : one
    implicit none
!   arguments
    integer, intent(in) :: m, n, lda
    real(wp_), intent(in) :: v(n),w(n)
    real(wp_), intent(inout) :: a(lda,n)
!   **********
!  
!   subroutine r1mpyq
!  
!   given an m by n matrix a, this subroutine computes a*q where
!   q is the product of 2*(n - 1) transformations
!  
!         gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
!  
!   and gv(i), gw(i) are givens rotations in the (i,n) plane which
!   eliminate elements in the i-th and n-th planes, respectively.
!   q itself is not given, rather the information to recover the
!   gv, gw rotations is supplied.
!  
!   the subroutine statement is
!  
!     subroutine r1mpyq(m,n,a,lda,v,w)
!  
!   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 must contain the matrix
!       to be postmultiplied by the orthogonal matrix q
!       described above. on output a*q has replaced a.
!  
!     lda is a positive integer input variable not less than m
!       which specifies the leading dimension of the array a.
!  
!     v is an input array of length n. v(i) must contain the
!       information necessary to recover the givens rotation gv(i)
!       described above.
!  
!     w is an input array of length n. w(i) must contain the
!       information necessary to recover the givens rotation gw(i)
!       described above.
!  
!   subroutines called
!  
!     fortran-supplied ... abs,sqrt
!  
!   argonne national laboratory. minpack project. march 1980.
!   burton s. garbow, kenneth e. hillstrom, jorge j. more
!  
!   **********
!   local variables
    integer :: i, j, nmj, nm1
    real(wp_) :: cs, sn, temp
!  
!   apply the first set of givens rotations to a.
!  
    nm1 = n - 1
    if (nm1 < 1) return
    do nmj = 1, nm1
      j = n - nmj
      if (abs(v(j)) > one) cs = one/v(j)
      if (abs(v(j)) > one) sn = sqrt(one-cs**2)
      if (abs(v(j)) <= one) sn = v(j)
      if (abs(v(j)) <= one) cs = sqrt(one-sn**2)
      do i = 1, m
        temp = cs*a(i,j) - sn*a(i,n)
        a(i,n) = sn*a(i,j) + cs*a(i,n)
        a(i,j) = temp
      end do
    end do
!  
!   apply the second set of givens rotations to a.
!  
    do j = 1, nm1
      if (abs(w(j)) > one) cs = one/w(j)
      if (abs(w(j)) > one) sn = sqrt(one-cs**2)
      if (abs(w(j)) <= one) sn = w(j)
      if (abs(w(j)) <= one) cs = sqrt(one-sn**2)
      do i = 1, m
        temp = cs*a(i,j) + sn*a(i,n)
        a(i,n) = -sn*a(i,j) + cs*a(i,n)
        a(i,j) = temp
      end do
    end do
  end subroutine r1mpyq

  subroutine r1updt(m,n,s,ls,u,v,w,sing)
    use const_and_precisions, only : zero, one, giant=>comp_huge
    implicit none
!   arguments
    integer, intent(in) :: m, n, ls
    logical, intent(out) :: sing
    real(wp_), intent(in) :: u(m)
    real(wp_), intent(out) :: w(m)
    real(wp_), intent(inout) :: s(ls), v(n)
!   **********
!  
!   subroutine r1updt
!  
!   given an m by n lower trapezoidal matrix s, an m-vector u,
!   and an n-vector v, the problem is to determine an
!   orthogonal matrix q such that
!  
!                 t
!         (s + u*v )*q
!  
!   is again lower trapezoidal.
!  
!   this subroutine determines q as the product of 2*(n - 1)
!   transformations
!  
!         gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
!  
!   where gv(i), gw(i) are givens rotations in the (i,n) plane
!   which eliminate elements in the i-th and n-th planes,
!   respectively. q itself is not accumulated, rather the
!   information to recover the gv, gw rotations is returned.
!  
!   the subroutine statement is
!  
!     subroutine r1updt(m,n,s,ls,u,v,w,sing)
!  
!   where
!  
!     m is a positive integer input variable set to the number
!       of rows of s.
!  
!     n is a positive integer input variable set to the number
!       of columns of s. n must not exceed m.
!  
!     s is an array of length ls. on input s must contain the lower
!       trapezoidal matrix s stored by columns. on output s contains
!       the lower trapezoidal matrix produced as described above.
!  
!     ls is a positive integer input variable not less than
!       (n*(2*m-n+1))/2.
!  
!     u is an input array of length m which must contain the
!       vector u.
!  
!     v is an array of length n. on input v must contain the vector
!       v. on output v(i) contains the information necessary to
!       recover the givens rotation gv(i) described above.
!  
!     w is an output array of length m. w(i) contains information
!       necessary to recover the givens rotation gw(i) described
!       above.
!  
!     sing is a logical output variable. sing is set true if any
!       of the diagonal elements of the output s are zero. otherwise
!       sing is set false.
!  
!   subprograms called
!  
!     fortran-supplied ... abs,sqrt
!  
!   argonne national laboratory. minpack project. march 1980.
!   burton s. garbow, kenneth e. hillstrom, jorge j. more,
!   john l. nazareth
!  
!   **********
!   local variables
    integer :: i, j, jj, l, nmj, nm1
    real(wp_) :: cs, cotan, sn, tn, tau, temp
!   parameters
    real(wp_), parameter :: p5=5.0e-1_wp_, p25=2.5e-1_wp_
!  
!   initialize the diagonal element pointer.
!  
    jj = (n*(2*m - n + 1))/2 - (m - n)
!  
!   move the nontrivial part of the last column of s into w.
!  
    l = jj
    do i = n, m
      w(i) = s(l)
      l = l + 1
    end do
!  
!   rotate the vector v into a multiple of the n-th unit vector
!   in such a way that a spike is introduced into w.
!  
    nm1 = n - 1
    do nmj = 1, nm1
      j = n - nmj
      jj = jj - (m - j + 1)
      w(j) = zero
      if (v(j) /= zero) then
!  
!       determine a givens rotation which eliminates the
!       j-th element of v.
!  
        if (abs(v(n)) < abs(v(j))) then
          cotan = v(n)/v(j)
          sn = p5/sqrt(p25+p25*cotan**2)
          cs = sn*cotan
          tau = one
          if (abs(cs)*giant > one) tau = one/cs
        else
          tn = v(j)/v(n)
          cs = p5/sqrt(p25+p25*tn**2)
          sn = cs*tn
          tau = sn
        end if
!  
!       apply the transformation to v and store the information
!       necessary to recover the givens rotation.
!  
        v(n) = sn*v(j) + cs*v(n)
        v(j) = tau
!  
!       apply the transformation to s and extend the spike in w.
!  
        l = jj
        do i = j, m
          temp = cs*s(l) - sn*w(i)
          w(i) = sn*s(l) + cs*w(i)
          s(l) = temp
          l = l + 1
        end do
      end if
    end do
!  
!   add the spike from the rank 1 update to w.
!  
    do i = 1, m
      w(i) = w(i) + v(n)*u(i)
    end do
!  
!   eliminate the spike.
!  
    sing = .false.
    do j = 1, nm1
      if (w(j) /= zero) then
!  
!       determine a givens rotation which eliminates the
!       j-th element of the spike.
!  
        if (abs(s(jj)) < abs(w(j))) then
          cotan = s(jj)/w(j)
          sn = p5/sqrt(p25+p25*cotan**2)
          cs = sn*cotan
          tau = one
          if (abs(cs)*giant > one) tau = one/cs
        else
          tn = w(j)/s(jj)
          cs = p5/sqrt(p25+p25*tn**2)
          sn = cs*tn
          tau = sn
        end if
!  
!       apply the transformation to s and reduce the spike in w.
!  
        l = jj
        do i = j, m
          temp = cs*s(l) + sn*w(i)
          w(i) = -sn*s(l) + cs*w(i)
          s(l) = temp
          l = l + 1
        end do
!  
!       store the information necessary to recover the
!       givens rotation.
!  
        w(j) = tau
      end if
!  
!     test for zero diagonal elements in the output s.
!  
      if (s(jj) == zero) sing = .true.
      jj = jj + (m - j + 1)
    end do
!  
!   move w back into the last column of the output s.
!  
    l = jj
    do i = n, m
      s(l) = w(i)
      l = l + 1
    end do
    if (s(jj) == zero) sing = .true.
!  
  end subroutine r1updt

end module minpack