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



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, 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



SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag)
 
! N.B. Arguments LDR & WA have been removed.

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(:)