subroutine snsqe(fcn,jac,iopt,n,x,fvec,tol,nprint,info,wa,lwa
     1 ,maxfev)
c***begin prologue  snsqe
c***revision date  850627   (yymmdd)
c***latest revision by M. Sagen
c        changed to double precision
c        allows maximum number of iterations to be read in
c
c***category no.  c5a,c5b
c***keywords easy-to-use,nonlinear square system, zero, powell
c       hybrid method
c***date written  march 1980
c***author  hiebert k.l. (sla)
c***purpose
c       snsqe is the easy-to-use version of snsq which finds a zero
c       of a system of n nonlinear functions in n variables by a modi-
c       fication of powell hybrid method.  this code is the combination
c       of the minpack codes (argonne) hybrd1 and hybrj1.
c***description
c
c
c 1. purpose.
c
c
c       the purpose of snsqe is to find a zero of a system of n non-
c       linear functions in n variables by a modification of the powell
c       hybrid method.  this is done by using the more general nonlinear
c       equation solver snsq.  the user must provide a subroutine which
c       calculates the functions.  the user has the option of either to
c       provide a subroutine which calculates the jacobian or to let the
c       code calculate it by a forward-difference approximation.  this
c       code is the combination of the minpack codes (argonne) hybrd1
c       and hybrj1.
c
c
c 2. subroutine and type statements.
c
c       subroutine snsqe(fcn,jac,iopt,n,x,fvec,tol,nprint,info,
c      *                  wa,lwa,maxfev)
c       integer iopt,n,nprint,info,lwa,maxfev
c       real tol
c       real x(n),fvec(n),wa(lwa)
c       external fcn,jac
c
c
c 3. parameters.
c
c       parameters designated as input parameters must be specified on
c       entry to snsqe and are not changed on exit, while parameters
c       designated as output parameters need not be specified on entry
c       and are set to appropriate values on exit from snsqe.
c
c       fcn is the name of the user-supplied subroutine which calculates
c         the functions.  fcn must be declared in an external statement
c         in the user calling program, and should be written as follows.
c
c         subroutine fcn(n,x,fvec,iflag)
c         integer n,iflag
c         real x(n),fvec(n)
c         ----------
c         calculate the functions at x and
c         return this vector in fvec.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by fcn unless the
c         user wants to terminate execution of snsqe.  in this case set
c         iflag to a negative integer.
c
c       jac is the name of the user-supplied subroutine which calculates
c         the jacobian.  if iopt=1, then jac must be declared in an
c         external statement in the user calling program, and should be
c         written as follows.
c
c         subroutine jac(n,x,fvec,fjac,ldfjac,iflag)
c         integer n,ldfjac,iflag
c         real x(n),fvec(n),fjac(ldfjac,n)
c         ----------
c         calculate the jacobian at x and return this
c         matrix in fjac.  fvec contains the function
c         values at x and should not be altered.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by jac unless the
c         user wants to terminate execution of snsqe. in this case set
c         iflag to a negative integer.
c
c         if iopt=2, jac can be ignored (treat it as a dummy argument).
c
c       iopt is an input variable which specifies how the jacobian will
c         be calculated.  if iopt=1, then the user must supply the
c         jacobian through the subroutine jac.  if iopt=2, then the
c         code will approximate the jacobian by forward-differencing.
c
c       n is a positive integer input variable set to the number of
c         functions and variables.
c
c       x is an array of length n.  on input x must contain an initial
c         estimate of the solution vector.  on output x contains the
c         final estimate of the solution vector.
c
c       fvec is an output array of length n which contains the functions
c         evaluated at the output x.
c
c       tol is a nonnegative input variable.  termination occurs when
c         the algorithm estimates that the relative error between x and
c         the solution is at most tol.  section 4 contains more details
c         about tol.
c
c       nprint is an integer input variable that enables controlled
c         printing of iterates if it is positive.  in this case, fcn is
c         called with iflag = 0 at the beginning of the first iteration
c         and every nprint iterations thereafter and immediately prior
c         to return, with x and fvec available for printing. appropriate
c         print statements must be added to fcn(see example).  if nprint
c         is not positive, no special calls of fcn with iflag = 0 are
c         made.
c
c       info is an integer output variable.  if the user has terminated
c         execution, info is set to the (negative) value of iflag.  see
c         description of fcn and jac. otherwise, info is set as follows.
c
c         info = 0  improper input parameters.
c
c         info = 1  algorithm estimates that the relative error between
c                   x and the solution is at most tol.
c
c         info = 2  number of calls to fcn has reached or exceeded
c                   maxfev if maxfev > 0.  If maxfev .le. 0 then
c                   number of calls is 100*(n+1) for iopt = 1 and
c                   200*(n+1) for iopt = 2. 
c
c         info = 3  tol is too small.  no further improvement in the
c                   approximate solution x is possible.
c
c         info = 4  iteration is not making good progress.
c
c         sections 4 and 5 contain more details about info.
c
c       wa is a work array of length lwa.
c
c       lwa is a positive integer input variable not less than
c         (3*n**2+13*n))/2.
c
c       maxfev is the maximum number of iterations
c
c 4. successful completion.
c
c       the accuracy of snsqe is controlled by the convergence parame-
c       ter tol.  this parameter is used in a test which makes a compar-
c       ison between the approximation x and a solution xsol.  snsqe
c       terminates when the test is satisfied.  if tol is less than the
c       machine precision (as defined by the  function r1mach(4)), then
c       snsqe only attempts to satisfy the test defined by the machine
c       precision.  further progress is not usually possible.  unless
c       high precision solutions are required, the recommended value
c       for tol is the square root of the machine precision.
c
c       the test assumes that the functions are reasonably well behaved,
c       and, if the jacobian is supplied by the user, that the functions
c       and the jacobian are coded consistently. if these conditions are
c       not satisfied, then snsqe may incorrectly indicate convergence.
c       the coding of the jacobian can be checked by the subroutine
c       chkder.  if the jacobian is coded correctly or iopt=2, then
c       the validity of the answer can be checked, for example, by
c       rerunning snsqe with a tighter tolerance.
c
c       convergence test.  if enorm(z) denotes the euclidean norm of a
c         vector z, then this test attempts to guarantee that
c
c               enorm(x-xsol) .le. tol*enorm(xsol).
c
c         if this condition is satisfied with tol = 10**(-k), then the
c         larger components of x have k significant decimal digits and
c         info is set to 1.  there is a danger that the smaller compo-
c         nents of x may have large relative errors, but the fast rate
c         of convergence of snsqe usually avoids this possibility.
c
c
c 5. unsuccessful completion.
c
c       unsuccessful termination of snsqe can be due to improper input
c       parameters, arithmetic interrupts, an excessive number of func-
c       tion evaluations, errors in the functions, or lack of good prog-
c       ress.
c
c       improper input parameters.  info is set to 0 if iopt .lt. 1, or
c         iopt .gt. 2, or n .le. 0, or tol .lt. 0.e0, or
c         lwa .lt. (3*n**2+13*n)/2.
c
c       arithmetic interrupts.  if these interrupts occur in the fcn
c         subroutine during an early stage of the computation, they may
c         be caused by an unacceptable choice of x by snsqe.  in this
c         case, it may be possible to remedy the situation by not evalu-
c         ating the functions here, but instead setting the components
c         of fvec to numbers that exceed those in the initial fvec.
c
c       excessive number of function evaluations.  if the number of
c         calls to fcn reaches maxfev for maxfev .gt. 0 or 100*(n+1) for
c         for iopt=1 or 200*(n+1) for iopt=2 when maxfev .le. 0,
c         then this indicates that the routine is converging
c         very slowly as measured by the progress of fvec, and info is
c         set to 2.  this situation should be unusual because, as
c         indicated below, lack of good progress is usually diagnosed
c         earlier by snsqe, causing termination with info = 4.
c
c       errors in the functions.  when iopt=2, the choice of step length
c         in the forward-difference approximation to the jacobian
c         assumes that the relative errors in the functions are of the
c         order of the machine precision.  if this is not the case,
c         snsqe may fail (usually with info = 4).  the user should
c         then either use snsq and set the step length or use iopt=1
c         and supply the jacobian.
c
c       lack of good progress.  snsqe searches for a zero of the system
c         by minimizing the sum of the squares of the functions.  in so
c         doing, it can become trapped in a region where the minimum
c         does not correspond to a zero of the system and, in this situ-
c         ation, the iteration eventually fails to make good progress.
c         in particular, this will happen if the system does not have a
c         zero.  if the system has a zero, rerunning snsqe from a dif-
c         ferent starting point may be helpful.
c
c
c 6. characteristics of the algorithm.
c
c       snsqe is a modification of the powell hybrid method.  two of
c       its main characteristics involve the choice of the correction as
c       a convex combination of the newton and scaled gradient direc-
c       tions, and the updating of the jacobian by the rank-1 method of
c       broyden.  the choice of the correction guarantees (under reason-
c       able conditions) global convergence for starting points far from
c       the solution and a fast rate of convergence.  the jacobian is
c       calculated at the starting point by either the user-supplied
c       subroutine or a forward-difference approximation, but it is not
c       recalculated until the rank-1 method fails to produce satis-
c       factory progress.
c
c       timing.  the time required by snsqe to solve a given problem
c         depends on n, the behavior of the functions, the accuracy
c         requested, and the starting point.  the number of arithmetic
c         operations needed by snsqe is about 11.5*(n**2) to process
c         each evaluation of the functions (call to fcn) and 1.3*(n**3)
c         to process each evaluation of the jacobian (call to jac,
c         if iopt = 1).  unless fcn and jac can be evaluated quickly,
c         the timing of snsqe will be strongly influenced by the time
c         spent in fcn and jac.
c
c       storage.  snsqe requires (3*n**2 + 17*n)/2 single precision
c         storage locations, in addition to the storage required by the
c         program.  there are no internally declared storage arrays.
c
c
c 7. example.
c
c       the problem is to determine the values of x(1), x(2), ..., x(9),
c       which solve the system of tridiagonal equations
c
c       (3-2*x(1))*x(1)           -2*x(2)                   = -1
c               -x(i-1) + (3-2*x(i))*x(i)         -2*x(i+1) = -1, i=2-8
c                                   -x(8) + (3-2*x(9))*x(9) = -1
c
c       **********
c
c       program test(input,output,tape6=output)
c c
c c     driver for snsqe example.
c c
c       integer j,n,iopt,nprint,info,lwa,nwrite,maxfev
c       real tol,fnorm
c       real x(9),fvec(9),wa(180)
c       real enorm,r1mach
c       external fcn
c       data nwrite /6/
c c
c       iopt = 2
c       n = 9
c c
c c     the following starting values provide a rough solution.
c c
c       do 10 j = 1, 9
c          x(j) = -1.e0
c    10    continue
c
c       lwa = 180
c       nprint = 0
c       maxfev = 20
c c
c c     set tol to the square root of the machine precision.
c c     unless high precision solutions are required,
c c     this is the recommended setting.
c c
c       tol = sqrt(r1mach(4))
c c
c       call snsqe(fcn,jac,iopt,n,x,fvec,tol,nprint,info,wa,lwa,maxfev)
c       fnorm = enorm(n,fvec)
c       write (nwrite,1000) fnorm,info,(x(j),j=1,n)
c       stop
c  1000 format (5x,31h final l2 norm of the residuals,e15.7 //
c      *        5x,15h exit parameter,16x,i10 //
c      *        5x,27h final approximate solution // (5x,3e15.7))
c       end
c       subroutine fcn(n,x,fvec,iflag)
c       integer n,iflag
c       real x(n),fvec(n)
c       integer k
c       real one,temp,temp1,temp2,three,two,zero
c       data zero,one,two,three /0.e0,1.e0,2.e0,3.e0/
c c
c       do 10 k = 1, n
c          temp = (three - two*x(k))*x(k)
c          temp1 = zero
c          if (k .ne. 1) temp1 = x(k-1)
c          temp2 = zero
c          if (k .ne. n) temp2 = x(k+1)
c          fvec(k) = temp - temp1 - two*temp2 + one
c    10    continue
c       return
c       end
c
c       results obtained with different compilers or machines
c       may be slightly different.
c
c       final l2 norm of the residuals  0.1192636e-07
c
c       exit parameter                         1
c
c       final approximate solution
c
c       -0.5706545e+00 -0.6816283e+00 -0.7017325e+00
c       -0.7042129e+00 -0.7013690e+00 -0.6918656e+00
c       -0.13257920e+00 -0.5960342e+00 -0.4164121e+00
c
c
c***references
c       powell, m. j. d.
c       a hybrid method for nonlinear equations.
c       numerical methods for nonlinear algebraic equations,
c       p. rabinowitz, editor.  gordon and breach, 1970.
c***routines called  snsq,xerror
c***end prologue  snsqe
      implicit real*8 (a-h,o-z)
      integer iopt,n,nprint,info,lwa
      real*8 tol,xtol
      real*8 x(n),fvec(n),wa(lwa)
      external fcn,jac
      integer index,j,lr,maxfev,ml,mode,mu,nfev,njev
      real*8 epsfcn,factor,one,zero
      data factor,one,zero /1.0d2,1.0d0,0.0d0/
c***first executable statement  snsqe
      info = 0
c
c     check the input parameters for errors.
c
      if (iopt .lt. 1 .or. iopt .gt. 2 .or. n .le. 0
     *    .or. tol .lt. zero .or. lwa .lt. (3*n**2 +13*n)/2)
     *   go to 20
c
c     call snsq.
c
      if (maxfev .le. 0) then
		 maxfev = 100*(n + 1)
		 if (iopt .eq. 2) maxfev = 2*maxfev
      endif
      xtol = tol
      ml = n - 1
      mu = n - 1
      epsfcn = zero
      mode = 2
      do 10 j = 1, n
	 wa(j) = one
   10    continue
      lr = (n*(n + 1))/2
      index=6*n+lr
      call snsq(fcn,jac,iopt,n,x,fvec,wa(index+1),n,xtol,maxfev,ml,mu,
     *           epsfcn,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 .eq. 5) info = 4
   20 continue
      if (info .eq. 0) call xerror(34hsnsqe  -- invalid input parameter.
     1,34,2,1)
      return
c
c     last card of subroutine snsqe.
c
      end
c======================================================================================		
	subroutine caxpy(n,ca,cx,incx,cy,incy)
c***begin prologue  caxpy
c***revision date  811015   (yymmdd)
c***category no.  f1a
c***keywords  blas,complex
c***date written  october 1979
c
c***changed to double precision june 1985 by m.s.
c
c***author lawson c. (jpl),hanson r. (sla),
c                            kincaid d. (u texas), krogh f. (jpl)
c***purpose
c  complex computation y = a*x + y
c***description
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       ca  complex scalar multiplier
c       cx  complex vector with n elements
c     incx  storage spacing between elements of cx
c       cy  complex vector with n elements
c     incy  storage spacing between elements of cy
c
c     --output--
c       cy  complex result (unchanged if n.le.0)
c
c     overwrite complex cy with complex  ca*cx + cy.
c     for i = 0 to n-1, replace  cy(ly+i*incy) with ca*cx(lx+i*incx) +
c       cy(ly+i*incy), where lx = 1 if incx .ge. 0, else lx = (-incx)*n
c       and ly is defined in a similar way using incy.
c
c
c***references
c  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c   *basic linear algebra subprograms for fortran usage*,
c  algorithm no. 539, transactions on mathematical software,
c  volume 5, number 3, september 1979, 308-323
c***routines called    (none)
c***end prologue  caxpy
c
      implicit real*8 (a-h,o-z)
      complex*16 cx(1),cy(1),ca
c***first executable statement  caxpy
      canorm = abs(dble(ca)) + abs(dimag(ca))
      if(n.le.0.or.canorm.eq.0.e0) return
      if(incx.eq.incy.and.incx.gt.0) go to 20
      kx = 1
      ky = 1
      if(incx.lt.0) kx = 1+(1-n)*incx
      if(incy.lt.0) ky = 1+(1-n)*incy
	  do 10 i = 1,n
	  cy(ky) = cy(ky) + ca*cx(kx)
	  kx = kx + incx
	  ky = ky + incy
   10 continue
      return
   20 continue
      ns = n*incx
	  do 30 i=1,ns,incx
	  cy(i) = ca*cx(i) + cy(i)
   30     continue
      return
      end
      subroutine cgedi(a,lda,n,ipvt,det,work,job)
c***begin prologue  cgedi
c***revision date  811015   (yymmdd)
c***category no.  f4f, f3
c***keywords  linpack,matrix,linear equations,complex,inverse,
c             determinant
c***date written  august 14, 1978
c
c***changed to double precision june 1985 by m.s.
c
c***author  moler c.b. (unm)
c***purpose
c  computes the determinant and inverse of a complex matrix
c   using the factors computed by cgeco or cgefa.
c***description
c     cgedi computes the determinant and inverse of a matrix
c     using the factors computed by cgeco or cgefa.
c
c     on entry
c
c        a       complex(lda, n)
c                the output from cgeco or cgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from cgeco or cgefa.
c
c        work    complex(n)
c                work vector.  contents destroyed.
c
c        job     integer
c                = 11   both determinant and inverse.
c                = 01   inverse only.
c                = 10   determinant only.
c
c     on return
c
c        a       inverse of original matrix if requested.
c                otherwise unchanged.
c
c        det     complex(2)
c                determinant of original matrix if requested.
c                otherwise not referenced.
c                determinant = det(1) * 10.0**det(2)
c                with  1.0 .le. cabs1(det(1)) .lt. 10.0
c                or  det(1) .eq. 0.0 .
c
c     error condition
c
c        a division by zero will occur if the input factor contains
c        a zero on the diagonal and the inverse is requested.
c        it will not occur if the subroutines are called correctly
c        and if cgeco has set rcond .gt. 0.0 or cgefa has set
c        info .eq. 0 .
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas caxpy,cscal,cswap
c     fortran abs,aimag,cmplx,mod,real
c
c***references
c  dongarra j.j., bunch j.r., moler c.b., stewart g.w.,
c   *linpack users  guide*, siam, 1979.
c***routines called  cswap,caxpy,cscal
c***end prologue  cgedi
      implicit real*8 (a-h,o-z)
      integer lda,n,ipvt(1),job
      complex*16 a(lda,1),det(2),work(1)
c
      complex*16 t
      real*8 ten
      integer i,j,k,kb,kp1,l,nm1
      complex*16 zdum
      real*8 cabs1
      cabs1(zdum) = abs(dble(zdum)) + abs(dimag(zdum))
c
c     compute determinant
c
c***first executable statement  cgedi
      if (job/10 .eq. 0) go to 70
	 det(1) = (1.0e0,0.0e0)
	 det(2) = (0.0e0,0.0e0)
	 ten = 10.0e0
	 do 50 i = 1, n
	    if (ipvt(i) .ne. i) det(1) = -det(1)
	    det(1) = a(i,i)*det(1)
c        ...exit
	    if (cabs1(det(1)) .eq. 0.0e0) go to 60
   10       if (cabs1(det(1)) .ge. 1.0e0) go to 20
	       det(1) = cmplx(ten,0.0e0)*det(1)
	       det(2) = det(2) - (1.0e0,0.0e0)
	    go to 10
   20       continue
   30       if (cabs1(det(1)) .lt. ten) go to 40
	       det(1) = det(1)/cmplx(ten,0.0e0)
	       det(2) = det(2) + (1.0e0,0.0e0)
	    go to 30
   40       continue
   50    continue
   60    continue
   70 continue
c
c     compute inverse(u)
c
      if (mod(job,10) .eq. 0) go to 150
	 do 100 k = 1, n
	    a(k,k) = (1.0e0,0.0e0)/a(k,k)
	    t = -a(k,k)
	    call cscal(k-1,t,a(1,k),1)
	    kp1 = k + 1
	    if (n .lt. kp1) go to 90
	    do 80 j = kp1, n
	       t = a(k,j)
	       a(k,j) = (0.0e0,0.0e0)
	       call caxpy(k,t,a(1,k),1,a(1,j),1)
   80       continue
   90       continue
  100    continue
c
c        form inverse(u)*inverse(l)
c
	 nm1 = n - 1
	 if (nm1 .lt. 1) go to 140
	 do 130 kb = 1, nm1
	    k = n - kb
	    kp1 = k + 1
	    do 110 i = kp1, n
	       work(i) = a(i,k)
	       a(i,k) = (0.0e0,0.0e0)
  110       continue
	    do 120 j = kp1, n
	       t = work(j)
	       call caxpy(n,t,a(1,j),1,a(1,k),1)
  120       continue
	    l = ipvt(k)
	    if (l .ne. k) call cswap(n,a(1,k),1,a(1,l),1)
  130    continue
  140    continue
  150 continue
      return
      end
      subroutine cgefa(a,lda,n,ipvt,info)
c***begin prologue  cgefa
c***revision date  811015   (yymmdd)
c***category no.  f4f
c***keywords  linpack,matrix,linear equations,complex,factor,
c             gauss elimination
c***date written  august 14, 1978
c
c***changed to double precision june 1985 by m.s.
c
c***author  moler c.b. (unm)
c***purpose
c   factors a complex matrix by gaussian elimination.
c***description
c     cgefa factors a complex matrix by gaussian elimination.
c
c     cgefa is usually called by cgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for cgeco) = (1 + 9/n)*(time for cgefa) .
c
c     on entry
c
c        a       complex(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that cgesl or cgedi will divide by zero
c                     if called.  use  rcond  in cgeco for a reliable
c                     indication of singularity.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas caxpy,cscal,icamax
c     fortran abs,aimag,real
c
c***references
c  dongarra j.j., bunch j.r., moler c.b., stewart g.w.,
c   *linpack users  guide*, siam, 1979.
c***routines called  caxpy,cscal,icamax
c***end prologue  cgefa
      implicit real*8 (a-h,o-z)
      integer lda,n,ipvt(1),info
      complex*16 a(lda,1)
c
      complex*16 t
      integer icamax,j,k,kp1,l,nm1
      complex*16 zdum
      real*8 cabs1
      cabs1(zdum) = abs(dble(zdum)) + abs(dimag(zdum))
c
c     gaussian elimination with partial pivoting
c
c***first executable statement  cgefa
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
	 kp1 = k + 1
c
c        find l = pivot index
c
	 l = icamax(n-k+1,a(k,k),1) + k - 1
	 ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
	 if (cabs1(a(l,k)) .eq. 0.0e0) go to 40
c
c           interchange if necessary
c
	    if (l .eq. k) go to 10
	       t = a(l,k)
	       a(l,k) = a(k,k)
	       a(k,k) = t
   10       continue
c
c           compute multipliers
c
	    t = -(1.0e0,0.0e0)/a(k,k)
	    call cscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
	    do 30 j = kp1, n
	       t = a(l,j)
	       if (l .eq. k) go to 20
		  a(l,j) = a(k,j)
		  a(k,j) = t
   20          continue
	       call caxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
	 go to 50
   40    continue
	    info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (cabs1(a(n,n)) .eq. 0.0e0) info = n
      return
      end
C=================================================================
      subroutine cmatin(a,b,indexa,indexb,ipivot,n,ndim,detr,deti)
c
c
c     finds the inverse of complex*16 matrix h=(rea,aima)=(a,b)
c
c     code provided by norman bardsley to kwon+tabakin s prog bopit
c      n=no of elements in a assigned(and in the sub=dimension)
c      matrix a,b(real+imag parts) destroyed by call
c      the other arguments are dummy just to save space
      implicit real*8 (a-h,o-z)
      dimension  indexa(n),indexb(n),ipivot(n)
      dimension a(ndim,ndim),b(ndim,ndim)
      equivalence (irow,jrow),(icolum,jcolum)
c
      detr=1.0
c      rhl add to define tempi for real*8 potential
      tempi=0.
      deti=0.0
      do 10 j=1,n
  10  ipivot(j)=0
      do 70 i=1,n
      xmax=0.0
      do 35 j=1,n
      if(ipivot(j)-1) 15, 35, 15
  15  do 30 k=1,n
      if (ipivot(k)-1) 20, 30, 90
  20  if(abs(a(j,k)).lt.1.0e-20) go to 30
      xjk=abs(a(j,k))+abs(b(j,k))
      if(xmax-xjk)  25,30,30
  25  irow=j
      icolum=k
      xmax=xjk
  30  continue
  35  continue
      ipivot(icolum)=ipivot(icolum)+1
      if(irow.eq.icolum) go to 50
      detr=-deti
      deti=-deti
      do 45 l=1,n
      swap=a(irow,l)
      swapi=b(irow,l)
      a(irow,l)=a(icolum,l)
      b(irow,l)=b(icolum,l)
      b(icolum,l)=swapi
  45  a(icolum,l)=swap
  50  indexa(i)=irow
      indexb(i)=icolum
      pivotr=a(icolum,icolum)
      pivoti=b(icolum,icolum)
      temp=detr*pivotr-deti*pivoti
      deti=detr*pivoti+deti*pivotr
      detr=temp
      a(icolum,icolum)=1.0
      b(icolum,icolum)=0.0
      if(pivoti.eq.0.0) tempr=1.0/pivotr
      if(pivoti.ne.0.0) tempr=pivotr/(pivotr*pivotr+pivoti*pivoti)
      if(pivoti.ne.0.0) tempi=-pivoti/(pivotr*pivotr+pivoti*pivoti)
      do 55 l=1,n
      temp=a(icolum,l)*tempr-b(icolum,l)*tempi
      b(icolum,l)=a(icolum,l)*tempi+b(icolum,l)*tempr
  55  a(icolum,l)=temp
      do 70 l1=1,n
      if(l1-icolum) 60, 70, 60
  60  tempa=a(l1,icolum)
      tempb=b(l1,icolum)
      a(l1,icolum)=0.0
      b(l1,icolum)=0.0
      do 65 l=1,n
      b(l1,l)=b(l1,l)-a(icolum,l)*tempb-b(icolum,l)*tempa
      a(l1,l)=a(l1,l)-a(icolum,l)*tempa+b(icolum,l)*tempb
  65  continue
  70  continue
      do 85 i=1,n
      l=n+1-i
      if(indexa(l)-indexb(l)) 75, 85, 75
  75  jrow=indexa(l)
      jcolum=indexb(l)
      do 80 k=1,n
      swap=a(k,jrow)
      swapi=b(k,jrow)
      a(k,jrow)=a(k,jcolum)
      b(k,jrow)=b(k,jcolum)
      a(k,jcolum)=swap
      b(k,jcolum)=swapi
  80  continue
  85  continue
  90  return
      end
      subroutine cscal(n,ca,cx,incx)
c***begin prologue  cscal
c***revision date  811015   (yymmdd)
c***category no.  f1a, m2
c***keywords  complex,blas,vector,scale
c***date written  october 1979
c
c***changed to double precision june 1985 by m.s.
c
c***author lawson c. (jpl),hanson r. (sla),
c                            kincaid d. (u texas), krogh f. (jpl)
c***purpose
c    complex vector scale x = a*x
c***description
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       ca  complex scale factor
c       cx  complex vector with n elements
c     incx  storage spacing between elements of cx
c
c     --output--
c    cscal  complex result (unchanged if n.le.0)
c
c     replace complex cx by complex ca*cx.
c     for i = 0 to n-1, replace cx(1+i*incx) with  ca * cx(1+i*incx)
c
c
c***references
c  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c   *basic linear algebra subprograms for fortran usage*,
c  algorithm no. 539, transactions on mathematical software,
c  volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  cscal
c
      implicit real*8 (a-h,o-z)
      complex*16 ca,cx(1)
c***first executable statement  cscal
      if(n .le. 0) return
      ns = n*incx
	  do 10 i = 1,ns,incx
	  cx(i) = ca*cx(i)
   10     continue
      return
      end
      subroutine cswap(n,cx,incx,cy,incy)
c***begin prologue  cswap
c***revision date  811015   (yymmdd)
c***category no.  f1a
c***keywords  complex,blas,vector,interchange
c***date written  october 1979
c
c***changed to double precision june 1985 by m.s.
c
c***author lawson c. (jpl),hanson r. (sla),
c                            kincaid d. (u texas), krogh f. (jpl)
c***purpose
c     interchange complex vectors
c***description
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       cx  complex vector with n elements
c     incx  storage spacing between elements of cx
c       cy  complex vector with n elements
c     incy  storage spacing between elements of cy
c
c     --output--
c       cx  input vector cy (unchanged if n.le.0)
c       cy  input vector cx (unchanged if n.le.0)
c
c     interchange complex cx and complex cy
c     for i = 0 to n-1, interchange  cx(lx+i*incx) and cy(ly+i*incy),
c     where lx = 1 if incx .gt. 0, else lx = (-incx)*n, and ly is
c     defined in a similar way using incy.
c
c
c***references
c  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c   *basic linear algebra subprograms for fortran usage*,
c  algorithm no. 539, transactions on mathematical software,
c  volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  cswap
c
      implicit real*8 (a-h,o-z)
      complex*16 cx(1),cy(1),ctemp
c***first executable statement  cswap
      if(n .le. 0)return
      if(incx.eq.incy.and.incx.gt.0) go to 20
      kx = 1
      ky = 1
      if(incx.lt.0) kx = 1+(1-n)*incx
      if(incy.lt.0) ky = 1+(1-n)*incy
	  do 10 i = 1,n
	  ctemp = cx(kx)
	  cx(kx) = cy(ky)
	  cy(ky) = ctemp
	  kx = kx + incx
	  ky = ky + incy
   10 continue
      return
   20 continue
      ns = n*incx
	  do 30 i=1,ns,incx
	  ctemp = cx(i)
	  cx(i) = cy(i)
	  cy(i) = ctemp
   30     continue
      return
      end
      subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
c***begin prologue  dogleg
c
c***changed to double precision
c
c***refer to  snsqe,snsq
c     **********
c
c     subroutine dogleg
c
c     given an m by n matrix a, an n by n nonsingular diagonal
c     matrix d, an m-vector b, and a positive number delta, the
c     problem is to determine the convex combination x of the
c     gauss-newton and scaled gradient directions that minimizes
c     (a*x - b) in the least squares sense, subject to the
c     restriction that the euclidean norm of d*x be at most delta.
c
c     this subroutine completes the solution of the problem
c     if it is provided with the necessary information from the
c     qr factorization of a. that is, if a = q*r, where q has
c     orthogonal columns and r is an upper triangular matrix,
c     then dogleg expects the full upper triangle of r and
c     the first n components of (q transpose)*b.
c
c     the subroutine statement is
c
c       subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
c
c     where
c
c       n is a positive integer input variable set to the order of r.
c
c       r is an input array of length lr which must contain the upper
c         triangular matrix r stored by rows.
c
c       lr is a positive integer input variable not less than
c         (n*(n+1))/2.
c
c       diag is an input array of length n which must contain the
c         diagonal elements of the matrix d.
c
c       qtb is an input array of length n which must contain the first
c         n elements of the vector (q transpose)*b.
c
c       delta is a positive input variable which specifies an upper
c         bound on the euclidean norm of d*x.
c
c       x is an output array of length n which contains the desired
c         convex combination of the gauss-newton direction and the
c         scaled gradient direction.
c
c       wa1 and wa2 are work arrays of length n.
c
c     subprograms called
c
c       minpack-supplied ... r1mach,enorm
c
c       fortran-supplied ... abs,dmax1,dmin1,sqrt
c
c     minpack. version of july 1978.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     epsmch is the machine precision.
c
c***routines called  r1mach,enorm
c***end prologue  dogleg
      implicit real*8 (a-h,o-z)
      integer n,lr
      real*8 delta
      real*8 r(lr),diag(n),qtb(n),x(n),wa1(n),wa2(n)
      integer i,j,jj,jp1,k,l
      real*8 alpha,bnorm,epsmch,gnorm,one,qnorm,sgnorm,sum,temp,zero
      real*8 r1mach,enorm
      data one,zero /1.0d0,0.0d0/
c***first executable statement  dogleg
      epsmch = r1mach(4)
c
c     first, calculate the gauss-newton direction.
c
      jj = (n*(n + 1))/2 + 1
      do 50 k = 1, n
	 j = n - k + 1
	 jp1 = j + 1
	 jj = jj - k
	 l = jj + 1
	 sum = zero
	 if (n .lt. jp1) go to 20
	 do 10 i = jp1, n
	    sum = sum + r(l)*x(i)
	    l = l + 1
   10       continue
   20    continue
	 temp = r(jj)
	 if (temp .ne. zero) go to 40
	 l = j
	 do 30 i = 1, j
	    temp = dmax1(temp,abs(r(l)))
	    l = l + n - i
   30       continue
	 temp = epsmch*temp
	 if (temp .eq. zero) temp = epsmch
   40    continue
	 x(j) = (qtb(j) - sum)/temp
   50    continue
c
c     test whether the gauss-newton direction is acceptable.
c
      do 60 j = 1, n
	 wa1(j) = zero
	 wa2(j) = diag(j)*x(j)
   60    continue
      qnorm = enorm(n,wa2)
      if (qnorm .le. delta) go to 140
c
c     the gauss-newton direction is not acceptable.
c     next, calculate the scaled gradient direction.
c
      l = 1
      do 80 j = 1, n
	 temp = qtb(j)
	 do 70 i = j, n
	    wa1(i) = wa1(i) + r(l)*temp
	    l = l + 1
   70       continue
	 wa1(j) = wa1(j)/diag(j)
   80    continue
c
c     calculate the norm of the scaled gradient direction,
c     normalize, and rescale the gradient.
c
      gnorm = enorm(n,wa1)
      sgnorm = zero
      alpha = delta/qnorm
      if (gnorm .eq. zero) go to 120
      do 90 j = 1, n
	 wa1(j) = (wa1(j)/gnorm)/diag(j)
   90    continue
c
c     calculate the point along the scaled gradient
c     at which the quadratic is minimized.
c
      l = 1
      do 110 j = 1, n
	 sum = zero
	 do 100 i = j, n
	    sum = sum + r(l)*wa1(i)
	    l = l + 1
  100       continue
	 wa2(j) = sum
  110    continue
      temp = enorm(n,wa2)
      sgnorm = (gnorm/temp)/temp
c
c     test whether the scaled gradient direction is acceptable.
c
      alpha = zero
      if (sgnorm .ge. delta) go to 120
c
c     the scaled gradient direction is not acceptable.
c     finally, calculate the point along the dogleg
c     at which the quadratic is minimized.
c
      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
  120 continue
c
c     form appropriate convex combination of the gauss-newton
c     direction and the scaled gradient direction.
c
      temp = (one - alpha)*dmin1(sgnorm,delta)
      do 130 j = 1, n
	 x(j) = temp*wa1(j) + alpha*x(j)
  130    continue
  140 continue
      return
c
c     last card of subroutine dogleg.
c
      end
C=================================================================
      subroutine eigenc(n,hr,hi,imag,iw,p0,q0,ia,ib,ip,uu,vv,ndim)
c
c    rhl slight revision of
c     inverse iteration method
c     code provided by norman bardsley to kwon+tabakin s code bopit
c      n=no pts used (and in sub dimension of matrix)
c      (p0,q0) = intial complex*16 energy guess
c      this prog does only 1 eigenvalue per call
c      imag = 0 no imag potent used     =1  both real*8 + imag used
c      iw = -1 for intermediate print out    = 0 else
c      uu = complex*16 wavefuntion in p space (eigenvector) vv=dummy
c      other variables are just dummy passed to save space
c      nb hr hi matrices are destroyed by call to eigenc, need redefine
c
      implicit real*8 (a-h,o-z)
      complex*16 uu,vv,zero,one,uniti,temp,e0,ee
      dimension  hr(ndim,ndim),hi(ndim,ndim),ia(n),ib(n),ip(n),
     1 uu(ndim),vv(ndim)
      data  ivec,iterx,conv/1,20,1.e-6/
      data zero,one,uniti/(0.0,0.0),(1.0,0.0),(0.0,1.0)/
      itern=0
      temp=zero
      e0=p0*one+q0*uniti*imag
c
c      rhl fix of problem with nused(n)=ndim
      if(n.gt.ndim)write(6,900)
      if(n.gt.ndim)stop
 900  format(' *****n gt ndim in eigenc **********')
c      set up intial values for eigenvectors and h-eguess matrix
      do 40 i=1,n
      hr(i,i)=hr(i,i)-p0
       uu(i)=zero
      if(imag.ne.0) hi(i,i)=hi(i,i)-q0
 40    continue
      do 39 j=1,n
      if(imag.eq.0)uu(j)=one
      if(imag.ne.0)uu(j)=(1.0,1.0)
 39   continue
	if(imag.ne.0.and.imag.ne.1)return
c
c     finds the inverse of h matrix.
c
      call cmatin(hr,hi,ia,ib,ip,n,ndim)
c
c      mult psi by inverse of h-eguess
  50  imax=0
      amax=0.0
      do 70 i=1,n
      vv(i)=zero
      do 60 j=1,n
  60  vv(i)=vv(i)+(hr(i,j)*one+hi(i,j)*imag*uniti)*uu(j)
      if (abs(vv(i)).lt.amax) go to 70
      imax=i
      amax=abs(vv(i))
  70  continue
c      renorm psi, determine new e = ee (complex*16 guess e=e0)
      do 80 i=1,n
  80  uu(i)=vv(i)/vv(imax)
c
      ee=e0+one/vv(imax)
      if(iw.eq.-1) write(6,105) itern,ee,uu(ivec)
      change=abs((ee-temp)/ee)
      if(itern.ge.5.and.change.lt.conv) go to 85
      if(itern.eq.iterx) write(6,103)
      if(itern.eq.iterx) go to 84
      temp=ee
      itern=itern+1
      go to 50
c
c85   write(6,520)
85    continue
c520   format( ' momentum space wave function')
c     write(6,530)
c530   format( ' mom. in inv. fermis',' real*8 part   ',' imagi
c  a  nary part')
c     do 500 ij=1,n
c     aq=uu(ij)
c     bq=aimag(uu(ij))
c     write(6,510) p(ij),aq,bq
c510   format(3e15.5)
 500   continue
  84  p0=ee
      q0=-uniti*ee
  103 format(/15h non-convergent)
  105 format(19h itern,ee,uu(ivec)=,i4,5x,4(1pe20.10))
      return
      end
      double precision function enorm(n,x)
c***begin prologue  enorm
c
c***changed to double precision
c
c***refer to  snsqe,snsq,snlse,snls
c     **********
c
c     function enorm
c
c     given an n-vector x, this function calculates the
c     euclidean norm of x.
c
c     the euclidean norm is computed by accumulating the sum of
c     squares in three different sums. the sums of squares for the
c     small and large components are scaled so that no overflows
c     occur. non-destructive underflows are permitted. underflows
c     and overflows do not occur in the computation of the unscaled
c     sum of squares for the intermediate components.
c     the definitions of small, intermediate and large components
c     depend on two constants, rdwarf and rgiant. the main
c     restrictions on these constants are that rdwarf**2 not
c     underflow and rgiant**2 not overflow. the constants
c     given here are suitable for every known computer.
c
c     the function statement is
c
c       real function enorm(n,x)
c
c     where
c
c       n is a positive integer input variable.
c
c       x is an input array of length n.
c
c     subprograms called
c
c       fortran-supplied ... abs,sqrt
c
c     minpack. version of october 1979.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c***routines called  (none)
c***end prologue  enorm
      implicit real*8 (a-h,o-z)
      integer n
      real*8 x(n)
      integer i
      real*8 agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs,x1max,x3max,
     *     zero
      data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/
c***first executable statement  enorm
      s1 = zero
      s2 = zero
      s3 = zero
      x1max = zero
      x3max = zero
      floatn = n
      agiant = rgiant/floatn
      do 90 i = 1, n
	 xabs = abs(x(i))
	 if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70
	    if (xabs .le. rdwarf) go to 30
c
c              sum for large components.
c
	       if (xabs .le. x1max) go to 10
		  s1 = one + s1*(x1max/xabs)**2
		  x1max = xabs
		  go to 20
   10          continue
		  s1 = s1 + (xabs/x1max)**2
   20          continue
	       go to 60
   30       continue
c
c              sum for small components.
c
	       if (xabs .le. x3max) go to 40
		  s3 = one + s3*(x3max/xabs)**2
		  x3max = xabs
		  go to 50
   40          continue
		  if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2
   50          continue
   60       continue
	    go to 80
   70    continue
c
c           sum for intermediate components.
c
	    s2 = s2 + xabs**2
   80    continue
   90    continue
c
c     calculation of norm.
c
      if (s1 .eq. zero) go to 100
	 enorm = x1max*sqrt(s1+(s2/x1max)/x1max)
	 go to 130
  100 continue
	 if (s2 .eq. zero) go to 110
	    if (s2 .ge. x3max)
     *         enorm = sqrt(s2*(one+(x3max/s2)*(x3max*s3)))
	    if (s2 .lt. x3max)
     *         enorm = sqrt(x3max*((s2/x3max)+(x3max*s3)))
	    go to 120
  110    continue
	    enorm = x3max*sqrt(s3)
  120    continue
  130 continue
      return
c
c     last card of function enorm.
c
      end
      subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
     *                  wa1,wa2)
c***begin prologue  fdjac1
c
c***changed to double precision
c
c***refer to snsqe,snsq
c
c     subroutine fdjac1
c
c     this subroutine computes a forward-difference approximation
c     to the n by n jacobian matrix associated with a specified
c     problem of n functions in n variables. if the jacobian has
c     a banded form, then function evaluations are saved by only
c     approximating the nonzero terms.
c
c     the subroutine statement is
c
c       subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
c                         wa1,wa2)
c
c     where
c
c       fcn is the name of the user-supplied subroutine which
c         calculates the functions. fcn must be declared
c         in an external statement in the user calling
c         program, and should be written as follows.
c
c         subroutine fcn(n,x,fvec,iflag)
c         integer n,iflag
c         real x(n),fvec(n)
c         ----------
c         calculate the functions at x and
c         return this vector in fvec.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by fcn unless
c         the user wants to terminate execution of fdjac1.
c         in this case set iflag to a negative integer.
c
c       n is a positive integer input variable set to the number
c         of functions and variables.
c
c       x is an input array of length n.
c
c       fvec is an input array of length n which must contain the
c         functions evaluated at x.
c
c       fjac is an output n by n array which contains the
c         approximation to the jacobian matrix evaluated at x.
c
c       ldfjac is a positive integer input variable not less than n
c         which specifies the leading dimension of the array fjac.
c
c       iflag is an integer variable which can be used to terminate
c         the execution of fdjac1. see description of fcn.
c
c       ml is a nonnegative integer input variable which specifies
c         the number of subdiagonals within the band of the
c         jacobian matrix. if the jacobian is not banded, set
c         ml to at least n - 1.
c
c       epsfcn is an input variable used in determining a suitable
c         step length for the forward-difference approximation. this
c         approximation assumes that the relative errors in the
c         functions are of the order of epsfcn. if epsfcn is less
c         than the machine precision, it is assumed that the relative
c         errors in the functions are of the order of the machine
c         precision.
c
c       mu is a nonnegative integer input variable which specifies
c         the number of superdiagonals within the band of the
c         jacobian matrix. if the jacobian is not banded, set
c         mu to at least n - 1.
c
c       wa1 and wa2 are work arrays of length n. if ml + mu + 1 is at
c         least n, then the jacobian is considered dense, and wa2 is
c         not referenced.
c
c     subprograms called
c
c       minpack-supplied ... r1mach
c
c       fortran-supplied ... abs,dmax1,sqrt
c
c     minpack. version of june 1979.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     epsmch is the machine precision.
c
c***routines called  r1mach
c***end prologue  fdjac1
      implicit real*8 (a-h,o-z)
      integer n,ldfjac,iflag,ml,mu
      real*8 epsfcn
      real*8 x(n),fvec(n),fjac(ldfjac,n),wa1(n),wa2(n)
      integer i,j,k,msum
      real*8 eps,epsmch,h,temp,zero
      real*8 r1mach
      data zero /0.0d0/
c***first executable statement  fdjac1
      epsmch = r1mach(4)
c
      eps = sqrt(dmax1(epsfcn,epsmch))
      msum = ml + mu + 1
      if (msum .lt. n) go to 40
c
c        computation of dense approximate jacobian.
c
	 do 20 j = 1, n
	    temp = x(j)
	    h = eps*abs(temp)
	    if (h .eq. zero) h = eps
	    x(j) = temp + h
	    call fcn(n,x,wa1,iflag)
	    if (iflag .lt. 0) go to 30
	    x(j) = temp
	    do 10 i = 1, n
	       fjac(i,j) = (wa1(i) - fvec(i))/h
   10          continue
   20       continue
   30    continue
	 go to 110
   40 continue
c
c        computation of banded approximate jacobian.
c
	 do 90 k = 1, msum
	    do 60 j = k, n, msum
	       wa2(j) = x(j)
	       h = eps*abs(wa2(j))
	       if (h .eq. zero) h = eps
	       x(j) = wa2(j) + h
   60          continue
	    call fcn(n,x,wa1,iflag)
	    if (iflag .lt. 0) go to 100
	    do 80 j = k, n, msum
	       x(j) = wa2(j)
	       h = eps*abs(wa2(j))
	       if (h .eq. zero) h = eps
	       do 70 i = 1, n
		  fjac(i,j) = zero
		  if (i .ge. j - mu .and. i .le. j + ml)
     *               fjac(i,j) = (wa1(i) - fvec(i))/h
   70             continue
   80          continue
   90       continue
  100    continue
  110 continue
      return
c
c     last card of subroutine fdjac1.
c
      end
      subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
c***begin prologue  fdjac2
c***refer to snlse,snls
c     **********
c
c     subroutine fdjac2
c
c     this subroutine computes a forward-difference approximation
c     to the m by n jacobian matrix associated with a specified
c     problem of m functions in n variables.
c
c     the subroutine statement is
c
c       subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
c
c     where
c
c       fcn is the name of the user-supplied subroutine which
c         calculates the functions. fcn must be declared
c         in an external statement in the user calling
c         program, and should be written as follows.
c
c         subroutine fcn(m,n,x,fvec,iflag)
c         integer m,n,iflag
c         double precision x(n),fvec(m)
c         ----------
c         calculate the functions at x and
c         return this vector in fvec.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by fcn unless
c         the user wants to terminate execution of fdjac2.
c         in this case set iflag to a negative integer.
c
c       m is a positive integer input variable set to the number
c         of functions.
c
c       n is a positive integer input variable set to the number
c         of variables. n must not exceed m.
c
c       x is an input array of length n.
c
c       fvec is an input array of length m which must contain the
c         functions evaluated at x.
c
c       fjac is an output m by n array which contains the
c         approximation to the jacobian matrix evaluated at x.
c
c       ldfjac is a positive integer input variable not less than m
c         which specifies the leading dimension of the array fjac.
c
c       iflag is an integer variable which can be used to terminate
c         the execution of fdjac2. see description of fcn.
c
c       epsfcn is an input variable used in determining a suitable
c         step length for the forward-difference approximation. this
c         approximation assumes that the relative errors in the
c         functions are of the order of epsfcn. if epsfcn is less
c         than the machine precision, it is assumed that the relative
c         errors in the functions are of the order of the machine
c         precision.
c
c       wa is a work array of length m.
c
c     subprograms called
c
c       user-supplied ...... fcn
c
c       minpack-supplied ... r1mach
c
c       fortran-supplied ... abs,dmax1,sqrt
c
c     minpack. version of june 1979.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     epsmch is the machine precision.
c
c***routines called  r1mach
c***end prologue  fdjac2
      integer m,n,ldfjac,iflag
      double precision epsfcn
      double precision x(n),fvec(m),fjac(ldfjac,n),wa(m)
      integer i,j
      double precision eps,epsmch,h,temp,zero
      double precision r1mach
      data zero /0.0d0/
c***first executable statement  fdjac2
      epsmch = r1mach(4)
c
      eps = sqrt(dmax1(epsfcn,epsmch))
      do 20 j = 1, n
	 temp = x(j)
	 h = eps*abs(temp)
	 if (h .eq. zero) h = eps
	 x(j) = temp + h
	 call fcn(m,n,x,wa,iflag)
	 if (iflag .lt. 0) go to 30
	 x(j) = temp
	 do 10 i = 1, m
	    fjac(i,j) = (wa(i) - fvec(i))/h
   10       continue
   20    continue
   30 continue
      return
c
c     last card of subroutine fdjac2.
c
      end
      subroutine fdump
c***begin prologue  fdump
c
c***changed to double precision
c
c***revision date  811015   (yymmdd)
c***category no.  q2,n2
c***keywords  symbolic dump,dump
c***date written  8/79
c***author  jones r.e. (sla)
c***purpose
c  symbolic dump (should be locally written)
c***description
c        ***note*** machine dependent routine
c        fdump is intended to be replaced by a locally written
c        version which produces a symbolic dump.  failing this,
c        it should be replaced by a version which prints the
c        subprogram nesting list.  note that this dump must be
c        printed on each of up to five files, as indicated by the
c        xgetua routine.  see xsetua and xgetua for details.
c
c     written by ron jones, with slatec common math library subcommittee
c     latest revision ---  23 may 1979
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  (none)
c***end prologue  fdump
c
c***first executable statement  fdump
      return
      end
C=================================================================
      integer function i1mach(i)
c***begin prologue  i1mach
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  q
c***keywords  machine constants,integer
c***date written   1975
c***author fox p.a.,hall a.d.,schryer n.l. (bell labs)
c***purpose
c returns integer machine dependent constants
c***description
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c   these machine constant routines must be activated for
c   a particular environment.
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     i1mach can be used to obtain machine-dependent parameters
c     for the local machine environment.  it is a function
c     subroutine with one (input) argument, and can be called
c     as follows, for example
c
c          k = i1mach(i)
c
c     where i=1,...,16.  the (output) value of k above is
c     determined by the (input) value of i.  the results for
c     various values of i are discussed below.
c
c  i/o unit numbers.
c    i1mach( 1) = the standard input unit.
c    i1mach( 2) = the standard output unit.
c    i1mach( 3) = the standard punch unit.
c    i1mach( 4) = the standard error message unit.
c
c  words.
c    i1mach( 5) = the number of bits per integer storage unit.
c    i1mach( 6) = the number of characters per integer storage unit.
c
c  integers.
c    assume integers are represented in the s-digit, base-a form
c
c               sign ( x(s-1)*a**(s-1) + ... + x(1)*a + x(0) )
c
c               where 0 .le. x(i) .lt. a for i=0,...,s-1.
c    i1mach( 7) = a, the base.
c    i1mach( 8) = s, the number of base-a digits.
c    i1mach( 9) = a**s - 1, the largest magnitude.
c
c  floating-point numbers.
c    assume floating-point numbers are represented in the t-digit,
c    base-b form
c               sign (b**e)*( (x(1)/b) + ... + (x(t)/b**t) )
c
c               where 0 .le. x(i) .lt. b for i=1,...,t,
c               0 .lt. x(1), and emin .le. e .le. emax.
c    i1mach(10) = b, the base.
c
c  single-precision
c    i1mach(11) = t, the number of base-b digits.
c    i1mach(12) = emin, the smallest exponent e.
c    i1mach(13) = emax, the largest exponent e.
c
c  double-precision
c    i1mach(14) = t, the number of base-b digits.
c    i1mach(15) = emin, the smallest exponent e.
c    i1mach(16) = emax, the largest exponent e.
c
c  to alter this function for a particular environment,
c  the desired set of data statements should be activated by
c  removing the c from column 1.  also, the values of
c  i1mach(1) - i1mach(4) should be checked for consistency
c  with the local operating system.
c
c***references
c  fox p.a., hall a.d., schryer n.l.,*framework for a portable library*,
c  acm transaction on mathematical software, vol. 4, no. 2,
c  june 1978, pp. 177-188.
c***routines called  (none)
c***end prologue  i1mach
c
      implicit real*8 (a-h,o-z)
      integer imach(16),output
c
c     equivalence (imach(4),output)
c
c     machine constants for the burroughs 1700 system.
c
c     data imach( 1) /    7 /
c     data imach( 2) /    2 /
c     data imach( 3) /    2 /
c     data imach( 4) /    2 /
c     data imach( 5) /   36 /
c     data imach( 6) /    4 /
c     data imach( 7) /    2 /
c     data imach( 8) /   33 /
c     data imach( 9) / z1ffffffff /
c     data imach(10) /    2 /
c     data imach(11) /   24 /
c     data imach(12) / -256 /
c     data imach(13) /  255 /
c     data imach(14) /   60 /
c     data imach(15) / -256 /
c     data imach(16) /  255 /
c
c     machine constants for the burroughs 5700 system.
c
c     data imach( 1) /   5 /
c     data imach( 2) /   6 /
c     data imach( 3) /   7 /
c     data imach( 4) /   6 /
c     data imach( 5) /  48 /
c     data imach( 6) /   6 /
c     data imach( 7) /   2 /
c     data imach( 8) /  39 /
c     data imach( 9) / o0007777777777777 /
c     data imach(10) /   8 /
c     data imach(11) /  13 /
c     data imach(12) / -50 /
c     data imach(13) /  76 /
c     data imach(14) /  26 /
c     data imach(15) / -50 /
c     data imach(16) /  76 /
c
c     machine constants for the burroughs 6700/7700 systems.
c
c     data imach( 1) /   5 /
c     data imach( 2) /   6 /
c     data imach( 3) /   7 /
c     data imach( 4) /   6 /
c     data imach( 5) /  48 /
c     data imach( 6) /   6 /
c     data imach( 7) /   2 /
c     data imach( 8) /  39 /
c     data imach( 9) / o0007777777777777 /
c     data imach(10) /   8 /
c     data imach(11) /  13 /
c     data imach(12) / -50 /
c     data imach(13) /  76 /
c     data imach(14) /  26 /
c     data imach(15) / -32754 /
c     data imach(16) /  32780 /
c
c     machine constants for the cdc 6000/7000 series.
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    7 /
c     data imach( 4) /6loutput/
c     data imach( 5) /   60 /
c     data imach( 6) /   10 /
c     data imach( 7) /    2 /
c     data imach( 8) /   48 /
c     data imach( 9) / 00007777777777777777b /
c     data imach(10) /    2 /
c     data imach(11) /   48 /
c     data imach(12) / -974 /
c     data imach(13) / 1070 /
c     data imach(14) /   96 /
c     data imach(15) / -927 /
c     data imach(16) / 1070 /
c
c     machine constants for the cray 1
c
c     data imach( 1) /   100 /
c     data imach( 2) /   101 /
c     data imach( 3) /   102 /
c     data imach( 4) /   101 /
c     data imach( 5) /    64 /
c     data imach( 6) /     8 /
c     data imach( 7) /     2 /
c     data imach( 8) /    63 /
c     data imach( 9) /  777777777777777777777b /
c     data imach(10) /     2 /
c     data imach(11) /    48 /
c     data imach(12) / -8192 /
c     data imach(13) /  8191 /
c     data imach(14) /    96 /
c     data imach(15) / -8192 /
c     data imach(16) /  8191 /
c
c     machine constants for the data general eclipse s/200
c
c     data imach( 1) /   11 /
c     data imach( 2) /   12 /
c     data imach( 3) /    8 /
c     data imach( 4) /   10 /
c     data imach( 5) /   16 /
c     data imach( 6) /    2 /
c     data imach( 7) /    2 /
c     data imach( 8) /   15 /
c     data imach( 9) /32767 /
c     data imach(10) /   16 /
c     data imach(11) /    6 /
c     data imach(12) /  -64 /
c     data imach(13) /   63 /
c     data imach(14) /   14 /
c     data imach(15) /  -64 /
c     data imach(16) /   63 /
c
c     machine constants for the harris 220
c
c     data imach( 1) /       5 /
c     data imach( 2) /       6 /
c     data imach( 3) /       0 /
c     data imach( 4) /       6 /
c     data imach( 5) /      24 /
c     data imach( 6) /       3 /
c     data imach( 7) /       2 /
c     data imach( 8) /      23 /
c     data imach( 9) / 8388607 /
c     data imach(10) /       2 /
c     data imach(11) /      23 /
c     data imach(12) /    -127 /
c     data imach(13) /     127 /
c     data imach(14) /      38 /
c     data imach(15) /    -127 /
c     data imach(16) /     127 /
c
c     machine constants for the honeywell 600/6000 series.
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /   43 /
c     data imach( 4) /    6 /
c     data imach( 5) /   36 /
c     data imach( 6) /    6 /
c     data imach( 7) /    2 /
c     data imach( 8) /   35 /
c     data imach( 9) / o377777777777 /
c     data imach(10) /    2 /
c     data imach(11) /   27 /
c     data imach(12) / -127 /
c     data imach(13) /  127 /
c     data imach(14) /   63 /
c     data imach(15) / -127 /
c     data imach(16) /  127 /
c
c     machine constants for the hp 2100
c     3 word double precision option with ftn4
c
c     data imach(1) /      5/
c     data imach(2) /      6 /
c     data imach(3) /      4 /
c     data imach(4) /      1 /
c     data imach(5) /     16 /
c     data imach(6) /      2 /
c     data imach(7) /      2 /
c     data imach(8) /     15 /
c     data imach(9) /  32767 /
c     data imach(10)/      2 /
c     data imach(11)/     23 /
c     data imach(12)/   -128 /
c     data imach(13)/    127 /
c     data imach(14)/     39 /
c     data imach(15)/   -128 /
c     data imach(16)/    127 /
c
c     machine constants for the hp 2100
c     4 word double precision option with ftn4
c
c     data imach(1) /      5 /
c     data imach(2) /      6 /
c     data imach(3) /      4 /
c     data imach(4) /      1 /
c     data imach(5) /     16 /
c     data imach(6) /      2 /
c     data imach(7) /      2 /
c     data imach(8) /     15 /
c     data imach(9) /  32767 /
c     data imach(10)/      2 /
c     data imach(11)/     23 /
c     data imach(12)/   -128 /
c     data imach(13)/    127 /
c     data imach(14)/     55 /
c     data imach(15)/   -128 /
c     data imach(16)/    127 /
c
c     machine constants for the ibm 360/370 series,
c     the xerox sigma 5/7/9, the sel systems 85/86, and
c     the perkin elmer (interdata) 7/32.
c
c     data imach( 1) /   5 /
c     data imach( 2) /   6 /
c     data imach( 3) /   7 /
c     data imach( 4) /   6 /
c     data imach( 5) /  32 /
c     data imach( 6) /   4 /
c     data imach( 7) /   2 /
c     data imach( 8) /  31 /
c     data imach( 9) / z7fffffff /
c     data imach(10) /  16 /
c     data imach(11) /   6 /
c     data imach(12) / -64 /
c     data imach(13) /  63 /
c     data imach(14) /  14 /
c     data imach(15) / -64 /
c     data imach(16) /  63 /
c
c     machine constants for the pdp-10 (ka processor).
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    5 /
c     data imach( 4) /    6 /
c     data imach( 5) /   36 /
c     data imach( 6) /    5 /
c     data imach( 7) /    2 /
c     data imach( 8) /   35 /
c     data imach( 9) / "377777777777 /
c     data imach(10) /    2 /
c     data imach(11) /   27 /
c     data imach(12) / -128 /
c     data imach(13) /  127 /
c     data imach(14) /   54 /
c     data imach(15) / -101 /
c     data imach(16) /  127 /
c
c     machine constants for the pdp-10 (ki processor).
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    5 /
c     data imach( 4) /    6 /
c     data imach( 5) /   36 /
c     data imach( 6) /    5 /
c     data imach( 7) /    2 /
c     data imach( 8) /   35 /
c     data imach( 9) / "377777777777 /
c     data imach(10) /    2 /
c     data imach(11) /   27 /
c     data imach(12) / -128 /
c     data imach(13) /  127 /
c     data imach(14) /   62 /
c     data imach(15) / -128 /
c     data imach(16) /  127 /
c
c     machine constants for pdp-11 fortran[s supporting
c     32-bit integer arithmetic.
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    5 /
c     data imach( 4) /    6 /
c     data imach( 5) /   32 /
c     data imach( 6) /    4 /
c     data imach( 7) /    2 /
c     data imach( 8) /   31 /
c     data imach( 9) / 2147483647 /
c     data imach(10) /    2 /
c     data imach(11) /   24 /
c     data imach(12) / -127 /
c     data imach(13) /  127 /
c     data imach(14) /   56 /
c     data imach(15) / -127 /
c     data imach(16) /  127 /
c
c     machine constants for pdp-11 fortran[s supporting
c     16-bit integer arithmetic.
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    5 /
c     data imach( 4) /    6 /
c     data imach( 5) /   16 /
c     data imach( 6) /    2 /
c     data imach( 7) /    2 /
c     data imach( 8) /   15 /
c     data imach( 9) / 32767 /
c     data imach(10) /    2 /
c     data imach(11) /   24 /
c     data imach(12) / -127 /
c     data imach(13) /  127 /
c     data imach(14) /   56 /
c     data imach(15) / -127 /
c     data imach(16) /  127 /
c
c     machine constants for the univac 1100 series. ftn compiler
c
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    1 /
c     data imach( 4) /    6 /
c     data imach( 5) /   36 /
c     data imach( 6) /    4 /
c     data imach( 7) /    2 /
c     data imach( 8) /   35 /
c     data imach( 9) / o377777777777 /
c     data imach(10) /    2 /
c     data imach(11) /   27 /
c     data imach(12) / -128 /
c     data imach(13) /  127 /
c     data imach(14) /   60 /
c     data imach(15) /-1024 /
c     data imach(16) / 1023 /
c
c     machine constants for the univac 1100 series. for compiler
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    7 /
c     data imach( 4) /    6 /
c     data imach( 5) /   36 /
c     data imach( 6) /    6 /
c     data imach( 7) /    2 /
c     data imach( 8) /   35 /
c     data imach( 9) / o377777777777 /
c     data imach(10) /    2 /
c     data imach(11) /   27 /
c     data imach(12) / -128 /
c     data imach(13) /  127 /
c     data imach(14) /   60 /
c     data imach(15) /-1024/
c     data imach(16) / 1023 /
c
c
c     machine constants for the vax 11/780
c
c     data imach(1) /    5 /
c     data imach(2) /    6 /
c     data imach(3) /    5 /
c     data imach(4) /    6 /
c     data imach(5) /   32 /
c     data imach(6) /    4 /
c     data imach(7) /    2 /
c     data imach(8) /   31 /
c     data imach(9) /2147483647 /
c     data imach(10)/    2 /
c     data imach(11)/   24 /
c     data imach(12)/ -127 /
c     data imach(13)/  127 /
c     data imach(14)/   56 /
c     data imach(15)/ -127 /
c     data imach(16)/  127 /
c
c     machine constants for the ridge 32

      data imach(1)  /      5 /
      data imach(2)  /      6 /
      data imach(3)  /      6 /
      data imach(4)  /      6 /
      data imach(5)  /     32 /
      data imach(6)  /      4 /
      data imach(7)  /      2 /
      data imach(8)  /     31 /
      data imach(9)  / 2147483647 /
      data imach(10) /      2 /
      data imach(11) /     23 /
      data imach(12) /   -127 /
      data imach(13) /    128 /
      data imach(14) /     53 /
      data imach(15) /  -1023 /
      data imach(16) /   1024 /
c
c***first executable statement  i1mach
c
      if (i .lt. 1  .or.  i .gt. 16) go to 10
c
      i1mach=imach(i)
      return
c
   10 continue
      write(output,9000)
9000  format(39h1error    1 in i1mach - i out of bounds )
c
c     call fdump
c
c
      stop
      end
      integer function icamax(n,cx,incx)
c***begin prologue  icamax
c***revision date  811015   (yymmdd)
c***category no.  f1a
c***keywords  blas,vector,complex,component,largest component
c***date written  october 1979
c
c***changed to double precision june 1985 by m.s.
c
c***author lawson c. (jpl),hanson r. (sla),
c                            kincaid d. (u texas), krogh f. (jpl)
c***purpose
c    find largest component of complex vector
c***description
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       cx  complex vector with n elements
c     incx  storage spacing between elements of cx
c
c     --output--
c   icamax  smallest index (zero if n.le.0)
c
c      returns the index of the component of cx having the
c      largest sum of magnitudes of real and imaginary parts.
c     icamax = first i, i = 1 to n, to minimize
c        abs(real(cx(1-incx+i*incx))) + abs(imag(cx(1-incx+i*incx)))
c
c
c***references
c  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c   *basic linear algebra subprograms for fortran usage*,
c  algorithm no. 539, transactions on mathematical software,
c  volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  icamax
c
      implicit real*8 (a-h,o-z)
      complex*16 cx(1)
c***first executable statement  icamax
      icamax = 0
      if(n.le.0) return
      icamax = 1
      if(n .le. 1) return
      ns = n*incx
      ii = 1
      summax = abs(dble(cx(1))) + abs(dimag(cx(1)))
	  do 20 i=1,ns,incx
	  sumri = abs(dble(cx(i))) + abs(dimag(cx(i)))
	  if(summax.ge.sumri) go to 10
	  summax = sumri
	  icamax = ii
   10     ii = ii + 1
   20     continue
      return
      end
      function j4save(iwhich,ivalue,iset)
c***begin prologue  j4save
c
c***changed to double precision
c
c***refer to  xerror
c     abstract
c        j4save saves and recalls several global variables needed
c        by the library error handling routines.
c
c     description of parameters
c      --input--
c        iwhich - index of item desired.
c                 = 1 refers to current error number.
c                 = 2 refers to current error control flag.
c                 = 3 refers to current unit number to which error
c                     messages are to be sent.  (0 means use standard.)
c                 = 4 refers to the maximum number of times any
c                     message is to be printed (as set by xermax).
c                 = 5 refers to the total number of units to which
c                     each error message is to be written.
c                 = 6 refers to the 2nd unit for error messages
c                 = 7 refers to the 3rd unit for error messages
c                 = 8 refers to the 4th unit for error messages
c                 = 9 refers to the 5th unit for error messages
c        ivalue - the value to be set for the iwhich-th parameter,
c                 if iset is .true. .
c        iset   - if iset=.true., the iwhich-th parameter will be
c                 given the value, ivalue.  if iset=.false., the
c                 iwhich-th parameter will be unchanged, and ivalue
c                 is a dummy parameter.
c      --output--
c        the (old) value of the iwhich-th parameter will be returned
c        in the function value, j4save.
c
c     written by ron jones, with slatec common math library subcommittee
c     adapted from bell laboratories port library error handler
c     latest revision ---  23 may 1979
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  (none)
c***end prologue  j4save
c
      logical iset
      integer iparam(9)
      data iparam(1),iparam(2),iparam(3),iparam(4)/0,1,0,10/
      data iparam(5)/1/
      data iparam(6),iparam(7),iparam(8),iparam(9)/0,0,0,0/
c***first executable statement  j4save
      j4save = iparam(iwhich)
      if (iset) iparam(iwhich) = ivalue
      return
      end
      subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sigma,wa1,
     *                 wa2)
c***begin prologue  lmpar
c***refer to snlse,snls
c     **********
c
c     subroutine lmpar
c
c     given an m by n matrix a, an n by n nonsingular diagonal
c     matrix d, an m-vector b, and a positive number delta,
c     the problem is to determine a value for the parameter
c     par such that if x solves the system
c
c           a*x = b ,     sqrt(par)*d*x = 0 ,
c
c     in the least squares sense, and dxnorm is the euclidean
c     norm of d*x, then either par is zero and
c
c           (dxnorm-delta) .le. 0.1*delta ,
c
c     or par is positive and
c
c           abs(dxnorm-delta) .le. 0.1*delta .
c
c     this subroutine completes the solution of the problem
c     if it is provided with the necessary information from the
c     qr factorization, with column pivoting, of a. that is, if
c     a*p = q*r, where p is a permutation matrix, q has orthogonal
c     columns, and r is an upper triangular matrix with diagonal
c     elements of nonincreasing magnitude, then lmpar expects
c     the full upper triangle of r, the permutation matrix p,
c     and the first n components of (q transpose)*b. on output
c     lmpar also provides an upper triangular matrix s such that
c
c            t   t                   t
c           p *(a *a + par*d*d)*p = s *s .
c
c     s is employed within lmpar and may be of separate interest.
c
c     only a few iterations are generally needed for convergence
c     of the algorithm. if, however, the limit of 10 iterations
c     is reached, then the output par will contain the best
c     value obtained so far.
c
c     the subroutine statement is
c
c       subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sigma,
c                        wa1,wa2)
c
c     where
c
c       n is a positive integer input variable set to the order of r.
c
c       r is an n by n array. on input the full upper triangle
c         must contain the full upper triangle of the matrix r.
c         on output the full upper triangle is unaltered, and the
c         strict lower triangle contains the strict upper triangle
c         (transposed) of the upper triangular matrix s.
c
c       ldr is a positive integer input variable not less than n
c         which specifies the leading dimension of the array r.
c
c       ipvt is an integer input array of length n which defines the
c         permutation matrix p such that a*p = q*r. column j of p
c         is column ipvt(j) of the identity matrix.
c
c       diag is an input array of length n which must contain the
c         diagonal elements of the matrix d.
c
c       qtb is an input array of length n which must contain the first
c         n elements of the vector (q transpose)*b.
c
c       delta is a positive input variable which specifies an upper
c         bound on the euclidean norm of d*x.
c
c       par is a nonnegative variable. on input par contains an
c         initial estimate of the levenberg-marquardt parameter.
c         on output par contains the final estimate.
c
c       x is an output array of length n which contains the least
c         squares solution of the system a*x = b, sqrt(par)*d*x = 0,
c         for the output par.
c
c       sigma is an output array of length n which contains the
c         diagonal elements of the upper triangular matrix s.
c
c       wa1 and wa2 are work arrays of length n.
c
c     subprograms called
c
c       minpack-supplied ... r1mach,enorm,qrsolv
c
c       fortran-supplied ... abs,dmax1,dmin1,sqrt
c
c     minpack. version of june 1979.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     dwarf is the smallest positive magnitude.
c
c***routines called  r1mach,enorm,qrsolv
c***end prologue  lmpar
      integer n,ldr
      integer ipvt(n)
      double precision delta,par
      double precision r(ldr,n),diag(n),qtb(n),x(n),
     &                 sigma(n),wa1(n),wa2(n)
      integer i,iter,j,jm1,jp1,k,l,nsing
      double precision dxnorm,dwarf,fp,gnorm,parc,parl,
     &                 paru,p1,p001,sum,temp,zero
      double precision r1mach,enorm
      data p1,p001,zero /1.0d-1,1.0d-3,0.0d0/
c***first executable statement  lmpar
      dwarf = r1mach(1)
c
c     compute and store in x the gauss-newton direction. if the
c     jacobian is rank-deficient, obtain a least squares solution.
c
      nsing = n
      do 10 j = 1, n
	 wa1(j) = qtb(j)
	 if (r(j,j) .eq. zero .and. nsing .eq. n) nsing = j - 1
	 if (nsing .lt. n) wa1(j) = zero
   10    continue
      if (nsing .lt. 1) go to 50
      do 40 k = 1, nsing
	 j = nsing - k + 1
	 wa1(j) = wa1(j)/r(j,j)
	 temp = wa1(j)
	 jm1 = j - 1
	 if (jm1 .lt. 1) go to 30
	 do 20 i = 1, jm1
	    wa1(i) = wa1(i) - r(i,j)*temp
   20       continue
   30    continue
   40    continue
   50 continue
      do 60 j = 1, n
	 l = ipvt(j)
	 x(l) = wa1(j)
   60    continue
c
c     initialize the iteration counter.
c     evaluate the function at the origin, and test
c     for acceptance of the gauss-newton direction.
c
      iter = 0
      do 70 j = 1, n
	 wa2(j) = diag(j)*x(j)
   70    continue
      dxnorm = enorm(n,wa2)
      fp = dxnorm - delta
      if (fp .le. p1*delta) go to 220
c
c     if the jacobian is not rank deficient, the newton
c     step provides a lower bound, parl, for the zero of
c     the function. otherwise set this bound to zero.
c
      parl = zero
      if (nsing .lt. n) go to 120
      do 80 j = 1, n
	 l = ipvt(j)
	 wa1(j) = diag(l)*(wa2(l)/dxnorm)
   80    continue
      do 110 j = 1, n
	 sum = zero
	 jm1 = j - 1
	 if (jm1 .lt. 1) go to 100
	 do 90 i = 1, jm1
	    sum = sum + r(i,j)*wa1(i)
   90       continue
  100    continue
	 wa1(j) = (wa1(j) - sum)/r(j,j)
  110    continue
      temp = enorm(n,wa1)
      parl = ((fp/delta)/temp)/temp
  120 continue
c
c     calculate an upper bound, paru, for the zero of the function.
c
      do 140 j = 1, n
	 sum = zero
	 do 130 i = 1, j
	    sum = sum + r(i,j)*qtb(i)
  130       continue
	 l = ipvt(j)
	 wa1(j) = sum/diag(l)
  140    continue
      gnorm = enorm(n,wa1)
      paru = gnorm/delta
      if (paru .eq. zero) paru = dwarf/dmin1(delta,p1)
c
c     if the input par lies outside of the interval (parl,paru),
c     set par to the closer endpoint.
c
      par = dmax1(par,parl)
      par = dmin1(par,paru)
      if (par .eq. zero) par = gnorm/dxnorm
c
c     beginning of an iteration.
c
  150 continue
	 iter = iter + 1
c
c        evaluate the function at the current value of par.
c
	 if (par .eq. zero) par = dmax1(dwarf,p001*paru)
	 temp = sqrt(par)
	 do 160 j = 1, n
	    wa1(j) = temp*diag(j)
  160       continue
	 call qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sigma,wa2)
	 do 170 j = 1, n
	    wa2(j) = diag(j)*x(j)
  170       continue
	 dxnorm = enorm(n,wa2)
	 temp = fp
	 fp = dxnorm - delta
c
c        if the function is small enough, accept the current value
c        of par. also test for the exceptional cases where parl
c        is zero or the number of iterations has reached 10.
c
	 if (abs(fp) .le. p1*delta
     *       .or. parl .eq. zero .and. fp .le. temp
     *            .and. temp .lt. zero .or. iter .eq. 10) go to 220
c
c        compute the newton correction.
c
	 do 180 j = 1, n
	    l = ipvt(j)
	    wa1(j) = diag(l)*(wa2(l)/dxnorm)
  180       continue
	 do 210 j = 1, n
	    wa1(j) = wa1(j)/sigma(j)
	    temp = wa1(j)
	    jp1 = j + 1
	    if (n .lt. jp1) go to 200
	    do 190 i = jp1, n
	       wa1(i) = wa1(i) - r(i,j)*temp
  190          continue
  200       continue
  210       continue
	 temp = enorm(n,wa1)
	 parc = ((fp/delta)/temp)/temp
c
c        depending on the sign of the function, update parl or paru.
c
	 if (fp .gt. zero) parl = dmax1(parl,par)
	 if (fp .lt. zero) paru = dmin1(paru,par)
c
c        compute an improved estimate for par.
c
	 par = dmax1(parl,par+parc)
c
c        end of an iteration.
c
	 go to 150
  220 continue
c
c     termination.
c
      if (iter .eq. 0) par = zero
      return
c
c     last card of subroutine lmpar.
c
      end
C=================================================================
      subroutine psibs(zsip,npt,nn,ldum,irmax,deltar,ptptw,pt,wt)
c      normalize p space wf and calculate r space wf via kwon-tabakin
      implicit real*8 (a-h,o-y)
      complex*16 zsip,psir,sump,sumr,cnorm
      dimension ptptw(npt),sbl(50),pt(npt),wt(npt),zsip(nn)
      data pi,hbarc/3.14159265,197.3289/
c     most of this sub follows kwon-tabakin 's sub wavfn
c
c      this portion normalizes the p-space w.f.
c
      sump =0.
c      nb psi ... are complex*16 numbers in this sub
      do 10 i=1,npt
      ptptw(i)=wt(i)*(pt(i)**2)
      sump=sump+zsip(i)*zsip(i)*ptptw(i)
c      nb this renom mixes real*8 and im parts of wf and can give nodes in
 10   continue
c     cnorm=1./csqrt(sump)
      cnorm=1./sqrt(sump)
      write(6,100) cnorm
      do 20 i=1,npt
      zsip(i)=cnorm*zsip(i)
 20   continue
      write(6,110)
      write(6,120) (pt(i),zsip(i),i=1,npt)
      write(6,130)
c
c     r space w.f.
c     irmax= no rvalues wanted     deltar= increment in r
      fac=sqrt(2./pi)
      do 40 i=1,irmax
      sumr=0.
      rval=deltar*(i-1)
      if(rval.eq.0.)rval=1.e-8
      do 30 j=1,npt
      pr=pt(j)*rval/hbarc
c     call spbesl(ldum,pr,sbl)
      sbl(ldum)=sin(pr)/pr
c    ********special l=0 only
      sumr=sumr+sbl(ldum)*zsip(j)*ptptw(j)
 30   continue
      psir=fac*sumr
      write(6,140) rval,psir
 40   continue
 100  format(' cnorm <complex*16> =',1pe14.6,e14.6,/)
 110  format(6x,1hp,10x,8hre(zsip),12x,8him(psip))
 120  format(6(1pe14.6,2x))
 130  format(/,' r (fermi)     re(psir)         im(psir)')
 140  format(0pf8.2,2x,1pe14.6,2x,e14.6)
      return
      end
      subroutine qform(m,n,q,ldq,wa)
c***begin prologue  qform
c
c***changed to double precision
c
c***refer to  snsqe,snsq
c     **********
c
c     subroutine qform
c
c     this subroutine proceeds from the computed qr factorization of
c     an m by n matrix a to accumulate the m by m orthogonal matrix
c     q from its factored form.
c
c     the subroutine statement is
c
c       subroutine qform(m,n,q,ldq,wa)
c
c     where
c
c       m is a positive integer input variable set to the number
c         of rows of a and the order of q.
c
c       n is a positive integer input variable set to the number
c         of columns of a.
c
c       q is an m by m array. on input the full lower trapezoid in
c         the first min(m,n) columns of q contains the factored form.
c         on output q has been accumulated into a square matrix.
c
c       ldq is a positive integer input variable not less than m
c         which specifies the leading dimension of the array q.
c
c       wa is a work array of length m.
c
c     subprograms called
c
c       fortran-supplied ... min0
c
c     minpack. version of january 1979.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     zero out upper triangle of q in the first min(m,n) columns.
c
c***routines called  (none)
c***end prologue  qform
      implicit real*8 (a-h,o-z)
      integer m,n,ldq
      real*8 q(ldq,m),wa(m)
      integer i,j,jm1,k,l,minmn,np1
      real*8 one,sum,temp,zero
      data one,zero /1.0d0,0.0d0/
c***first executable statement  qform
      minmn = min0(m,n)
      if (minmn .lt. 2) go to 30
      do 20 j = 2, minmn
	 jm1 = j - 1
	 do 10 i = 1, jm1
	    q(i,j) = zero
   10       continue
   20    continue
   30 continue
c
c     initialize remaining columns to those of the identity matrix.
c
      np1 = n + 1
      if (m .lt. np1) go to 60
      do 50 j = np1, m
	 do 40 i = 1, m
	    q(i,j) = zero
   40       continue
	 q(j,j) = one
   50    continue
   60 continue
c
c     accumulate q from its factored form.
c
      do 120 l = 1, minmn
	 k = minmn - l + 1
	 do 70 i = k, m
	    wa(i) = q(i,k)
	    q(i,k) = zero
   70       continue
	 q(k,k) = one
	 if (wa(k) .eq. zero) go to 110
	 do 100 j = k, m
	    sum = zero
	    do 80 i = k, m
	       sum = sum + q(i,j)*wa(i)
   80          continue
	    temp = sum/wa(k)
	    do 90 i = k, m
	       q(i,j) = q(i,j) - temp*wa(i)
   90          continue
  100       continue
  110    continue
  120    continue
      return
c
c     last card of subroutine qform.
c
      end
C=================================================================
      subroutine qlfun(l,z,ql)
c
c     legrendre functn of second kind for all l values
c     from tabakin s bopit
      implicit real*8 (a-h,o-z)
      dimension fac(15),ql(15)
      if(z.gt.1) go to 30
c     rhl mod use series for z gt 1 vs 10 as tabakin
      if(z.eq.1)write(6,900)
      if(z.eq.1)return
 900  format('qlfun called for z=1')
      ql(1)=0.5*dlog(abs((z+1)/(z-1)))
      ql(2)=z*ql(1)-1.
      if(l.le.2) go to 20
      do 10 n=2,l
      ql(n+1)=(2*n-1)*z*ql(n)-(n-1)*ql(n-1)
 10   ql(n+1)=ql(n+1)/n
 20   return
c
c     series expansion
c
 30   fac(1)=1.
      nmax=l+2
      zsqinv=1/z/z
      do 50 n=1,nmax
 50   fac(n+1)=fac(n)*n/(2.*n+1)
       do 70 lp1=1,l
      ll=lp1-1
      prefac=fac(ll+1)/(z**(ll+1))
      fk=1.
      alfa=0.5*ll
      beta=alfa-.5
      gama=ll+.5
      sum=1.
      do 60 k=1,100
      fk=fk*zsqinv*(alfa+k)*(beta+k)/(gama+k)/k
      if(fk.lt.1.e-17) go to 70
 60   sum=sum+fk
 70   ql(lp1)= prefac*sum
      return
      end
      subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,sigma,acnorm,wa)
c***begin prologue  qrfac
c
c***changed to double precision
c
c***refer to  snlse,snls,snsqe,snsq
c     **********
c
c     subroutine qrfac
c
c     this subroutine uses householder transformations with column
c     pivoting (optional) to compute a qr factorization of the
c     m by n matrix a. that is, qrfac determines an orthogonal
c     matrix q, a permutation matrix p, and an upper trapezoidal
c     matrix r with diagonal elements of nonincreasing magnitude,
c     such that a*p = q*r. the householder transformation for
c     column k, k = 1,2,...,min(m,n), is of the form
c
c                           t
c           i - (1/u(k))*u*u
c
c     where u has zeros in the first k-1 positions. the form of
c     this transformation and the method of pivoting first
c     appeared in the corresponding linpack subroutine.
c
c     the subroutine statement is
c
c       subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,sigma,acnorm,wa)
c
c     where
c
c       m is a positive integer input variable set to the number
c         of rows of a.
c
c       n is a positive integer input variable set to the number
c         of columns of a.
c
c       a is an m by n array. on input a contains the matrix for
c         which the qr factorization is to be computed. on output
c         the strict upper trapezoidal part of a contains the strict
c         upper trapezoidal part of r, and the lower trapezoidal
c         part of a contains a factored form of q (the non-trivial
c         elements of the u vectors described above).
c
c       lda is a positive integer input variable not less than m
c         which specifies the leading dimension of the array a.
c
c       pivot is a logical input variable. if pivot is set true,
c         then column pivoting is enforced. if pivot is set false,
c         then no column pivoting is done.
c
c       ipvt is an integer output array of length lipvt. ipvt
c         defines the permutation matrix p such that a*p = q*r.
c         column j of p is column ipvt(j) of the identity matrix.
c         if pivot is false, ipvt is not referenced.
c
c       lipvt is a positive integer input variable. if pivot is false,
c         then lipvt may be as small as 1. if pivot is true, then
c         lipvt must be at least n.
c
c       sigma is an output array of length n which contains the
c         diagonal elements of r.
c
c       acnorm is an output array of length n which contains the
c         norms of the corresponding columns of the input matrix a.
c         if this information is not needed, then acnorm can coincide
c         with sigma.
c
c       wa is a work array of length n. if pivot is false, then wa
c         can coincide with sigma.
c
c     subprograms called
c
c       minpack-supplied ... r1mach,enorm
c
c       fortran-supplied ... dmax1,sqrt,min0
c
c     minpack. version of december 1978.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     epsmch is the machine precision.
c
c***routines called  r1mach,enorm
c***end prologue  qrfac
      implicit real*8 (a-h,o-z)
      integer m,n,lda,lipvt
      integer ipvt(lipvt)
      logical pivot
      real*8 a(lda,n),sigma(n),acnorm(n),wa(n)
      integer i,j,jp1,k,kmax,minmn
      real*8 ajnorm,epsmch,one,p05,sum,temp,zero
      real*8 r1mach,enorm
      data one,p05,zero /1.0d0,5.0d-2,0.0d0/
c***first executable statement  qrfac
      epsmch = r1mach(4)
c
c     compute the initial column norms and initialize several arrays.
c
      do 10 j = 1, n
	 acnorm(j) = enorm(m,a(1,j))
	 sigma(j) = acnorm(j)
	 wa(j) = sigma(j)
	 if (pivot) ipvt(j) = j
   10    continue
c
c     reduce a to r with householder transformations.
c
      minmn = min0(m,n)
      do 110 j = 1, minmn
	 if (.not.pivot) go to 40
c
c        bring the column of largest norm into the pivot position.
c
	 kmax = j
	 do 20 k = j, n
	    if (sigma(k) .gt. sigma(kmax)) kmax = k
   20       continue
	 if (kmax .eq. j) go to 40
	 do 30 i = 1, m
	    temp = a(i,j)
	    a(i,j) = a(i,kmax)
	    a(i,kmax) = temp
   30       continue
	 sigma(kmax) = sigma(j)
	 wa(kmax) = wa(j)
	 k = ipvt(j)
	 ipvt(j) = ipvt(kmax)
	 ipvt(kmax) = k
   40    continue
c
c        compute the householder transformation to reduce the
c        j-th column of a to a multiple of the j-th unit vector.
c
	 ajnorm = enorm(m-j+1,a(j,j))
	 if (ajnorm .eq. zero) go to 100
	 if (a(j,j) .lt. zero) ajnorm = -ajnorm
	 do 50 i = j, m
	    a(i,j) = a(i,j)/ajnorm
   50       continue
	 a(j,j) = a(j,j) + one
c
c        apply the transformation to the remaining columns
c        and update the norms.
c
	 jp1 = j + 1
	 if (n .lt. jp1) go to 100
	 do 90 k = jp1, n
	    sum = zero
	    do 60 i = j, m
	       sum = sum + a(i,j)*a(i,k)
   60          continue
	    temp = sum/a(j,j)
	    do 70 i = j, m
	       a(i,k) = a(i,k) - temp*a(i,j)
   70          continue
	    if (.not.pivot .or. sigma(k) .eq. zero) go to 80
	    temp = a(j,k)/sigma(k)
	    sigma(k) = sigma(k)*sqrt(dmax1(zero,one-temp**2))
	    if (p05*(sigma(k)/wa(k))**2 .gt. epsmch) go to 80
	    sigma(k) = enorm(m-j,a(jp1,k))
	    wa(k) = sigma(k)
   80       continue
   90       continue
  100    continue
	 sigma(j) = -ajnorm
  110    continue
      return
c
c     last card of subroutine qrfac.
c
      end
      subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sigma,wa)
c***begin prologue  qrsolv
c***refer to  snlse,snls
c     **********
c
c     subroutine qrsolv
c
c     given an m by n matrix a, an n by n diagonal matrix d,
c     and an m-vector b, the problem is to determine an x which
c     solves the system
c
c           a*x = b ,     d*x = 0 ,
c
c     in the least squares sense.
c
c     this subroutine completes the solution of the problem
c     if it is provided with the necessary information from the
c     qr factorization, with column pivoting, of a. that is, if
c     a*p = q*r, where p is a permutation matrix, q has orthogonal
c     columns, and r is an upper triangular matrix with diagonal
c     elements of nonincreasing magnitude, then qrsolv expects
c     the full upper triangle of r, the permutation matrix p,
c     and the first n components of (q transpose)*b. the system
c     a*x = b, d*x = 0, is then equivalent to
c
c                  t       t
c           r*z = q *b ,  p *d*p*z = 0 ,
c
c     where x = p*z. if this system does not have full rank,
c     then a least squares solution is obtained. on output qrsolv
c     also provides an upper triangular matrix s such that
c
c            t   t               t
c           p *(a *a + d*d)*p = s *s .
c
c     s is computed within qrsolv and may be of separate interest.
c
c     the subroutine statement is
c
c       subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sigma,wa)
c
c     where
c
c       n is a positive integer input variable set to the order of r.
c
c       r is an n by n array. on input the full upper triangle
c         must contain the full upper triangle of the matrix r.
c         on output the full upper triangle is unaltered, and the
c         strict lower triangle contains the strict upper triangle
c         (transposed) of the upper triangular matrix s.
c
c       ldr is a positive integer input variable not less than n
c         which specifies the leading dimension of the array r.
c
c       ipvt is an integer input array of length n which defines the
c         permutation matrix p such that a*p = q*r. column j of p
c         is column ipvt(j) of the identity matrix.
c
c       diag is an input array of length n which must contain the
c         diagonal elements of the matrix d.
c
c       qtb is an input array of length n which must contain the first
c         n elements of the vector (q transpose)*b.
c
c       x is an output array of length n which contains the least
c         squares solution of the system a*x = b, d*x = 0.
c
c       sigma is an output array of length n which contains the
c         diagonal elements of the upper triangular matrix s.
c
c       wa is a work array of length n.
c
c     subprograms called
c
c       fortran-supplied ... abs,sqrt
c
c     minpack. version of december 1978.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     copy r and (q transpose)*b to preserve input and initialize s.
c     in particular, save the diagonal elements of r in x.
c
c***routines called  (none)
c***end prologue  qrsolv
      integer n,ldr
      integer ipvt(n)
      double precision r(ldr,n),diag(n),qtb(n),x(n),sigma(n),wa(n)
      integer i,j,jp1,k,kp1,l,nsing
      double precision cos,cotan,p5,p25,qtbpj,sin,sum,tan,temp,zero
      data p5,p25,zero /5.0d-1,2.5d-1,0.0d0/
c***first executable statement  qrsolv
      do 20 j = 1, n
	 do 10 i = j, n
	    r(i,j) = r(j,i)
   10       continue
	 x(j) = r(j,j)
	 wa(j) = qtb(j)
   20    continue
c
c     eliminate the diagonal matrix d using a givens rotation.
c
      do 100 j = 1, n
c
c        prepare the row of d to be eliminated, locating the
c        diagonal element using p from the qr factorization.
c
	 l = ipvt(j)
	 if (diag(l) .eq. zero) go to 90
	 do 30 k = j, n
	    sigma(k) = zero
   30       continue
	 sigma(j) = diag(l)
c
c        the transformations to eliminate the row of d
c        modify only a single element of (q transpose)*b
c        beyond the first n, which is initially zero.
c
	 qtbpj = zero
	 do 80 k = j, n
c
c           determine a givens rotation which eliminates the
c           appropriate element in the current row of d.
c
	    if (sigma(k) .eq. zero) go to 70
	    if (abs(r(k,k)) .ge. abs(sigma(k))) go to 40
	       cotan = r(k,k)/sigma(k)
	       sin = p5/sqrt(p25+p25*cotan**2)
	       cos = sin*cotan
	       go to 50
   40       continue
	       tan = sigma(k)/r(k,k)
	       cos = p5/sqrt(p25+p25*tan**2)
	       sin = cos*tan
   50       continue
c
c           compute the modified diagonal element of r and
c           the modified element of ((q transpose)*b,0).
c
	    r(k,k) = cos*r(k,k) + sin*sigma(k)
	    temp = cos*wa(k) + sin*qtbpj
	    qtbpj = -sin*wa(k) + cos*qtbpj
	    wa(k) = temp
c
c           accumulate the tranformation in the row of s.
c
	    kp1 = k + 1
	    if (n .lt. kp1) go to 70
	    do 60 i = kp1, n
	       temp = cos*r(i,k) + sin*sigma(i)
	       sigma(i) = -sin*r(i,k) + cos*sigma(i)
	       r(i,k) = temp
   60          continue
   70       continue
   80       continue
   90    continue
c
c        store the diagonal element of s and restore
c        the corresponding diagonal element of r.
c
	 sigma(j) = r(j,j)
	 r(j,j) = x(j)
  100    continue
c
c     solve the triangular system for z. if the system is
c     singular, then obtain a least squares solution.
c
      nsing = n
      do 110 j = 1, n
	 if (sigma(j) .eq. zero .and. nsing .eq. n) nsing = j - 1
	 if (nsing .lt. n) wa(j) = zero
  110    continue
      if (nsing .lt. 1) go to 150
      do 140 k = 1, nsing
	 j = nsing - k + 1
	 sum = zero
	 jp1 = j + 1
	 if (nsing .lt. jp1) go to 130
	 do 120 i = jp1, nsing
	    sum = sum + r(i,j)*wa(i)
  120       continue
  130    continue
	 wa(j) = (wa(j) - sum)/sigma(j)
  140    continue
  150 continue
c
c     permute the components of z back to components of x.
c
      do 160 j = 1, n
	 l = ipvt(j)
	 x(l) = wa(j)
  160    continue
      return
c
c     last card of subroutine qrsolv.
c
      end
      double precision function r1mach(i)
c***begin prologue  r1mach
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  q
c***keywords  machine constants
c***date written  1979
c***author  fox p.a.,.hall a.d., schryer n.l. (bell labs)
c***purpose
c returns single precision machine dependent constants
c***description
c
c     r1mach can be used to obtain machine-dependent parameters
c     for the local machine environment.  it is a function
c     subroutine with one (input) argument, and can be called
c     as follows, for example
c
c          a = r1mach(i)
c
c     where i=1,...,5.  the (output) value of a above is
c     determined by the (input) value of i.  the results for
c     various values of i are discussed below.
c
c  single-precision machine constants
c  r1mach(1) = b**(emin-1), the smallest positive magnitude.
c  r1mach(2) = b**emax*(1 - b**(-t)), the largest magnitude.
c  r1mach(3) = b**(-t), the smallest relative spacing.
c  r1mach(4) = b**(1-t), the largest relative spacing.
c  r1mach(5) = log10(b)
c
c  to alter this function for a particular environment,
c  the desired set of data statements should be activated by
c  removing the c from column 1.
c
c  where possible, octal or hexadecimal constants have been used
c  to specify the constants exactly which has in some cases
c  required the use of equivalent integer arrays.
c
c***references
c  fox p.a., hall a.d., schryer n.l.,*framework for a portable library*,
c  acm transaction on mathematical software, vol. 4, no. 2,
c  june 1978, pp. 177-188.
c***routines called  xerror
c***end prologue  r1mach
c
      implicit real*8 (a-h,o-z)
c      integer small(2)
c      integer large(2)
c      integer right(2)
c      integer diver(2)
c      integer log10(2)
c
      real*8 rmach(5)
c
c      equivalence (rmach(1),small(1))
c      equivalence (rmach(2),large(1))
c      equivalence (rmach(3),right(1))
c      equivalence (rmach(4),diver(1))
c      equivalence (rmach(5),log10(1))
c
c     machine constants for the burroughs 1700 system.
c
c     data rmach(1) / z400800000 /
c     data rmach(2) / z5ffffffff /
c     data rmach(3) / z4e9800000 /
c     data rmach(4) / z4ea800000 /
c     data rmach(5) / z500e730e8 /
c
c     machine constants for the burroughs 5700/6700/7700 systems.
c
c     data rmach(1) / o1771000000000000 /
c     data rmach(2) / o0777777777777777 /
c     data rmach(3) / o1311000000000000 /
c     data rmach(4) / o1301000000000000 /
c     data rmach(5) / o1157163034761675 /
c
c     machine constants for the cdc 6000/7000 series.
c
c     data rmach(1) / 00014000000000000000b /
c     data rmach(2) / 37767777777777777777b /
c     data rmach(3) / 16404000000000000000b /
c     data rmach(4) / 16414000000000000000b /
c     data rmach(5) / 17164642023241175720b /
c
c     machine constants for the cray 1
c
c     data rmach(1) / 200004000000000000000b /
c     data rmach(2) / 577777777777777777777b /
c     data rmach(3) / 377214000000000000000b /
c     data rmach(4) / 377224000000000000000b /
c     data rmach(5) / 377774642023241175720b /
c
c     machine constants for the data general eclipse s/200
c
c     note - it may be appropriate to include the following card -
c     static rmach(5)
c
c     data small/20k,0/,large/77777k,177777k/
c     data right/35420k,0/,diver/36020k,0/
c     data log10/40423k,42023k/
c
c     machine constants for the harris 220
c
c     data small(1),small(2) / [20000000, [00000201 /
c     data large(1),large(2) / [37777777, [00000177 /
c     data right(1),right(2) / [20000000, [00000352 /
c     data diver(1),diver(2) / [20000000, [00000353 /
c     data log10(1),log10(2) / [23210115, [00000377 /
c
c     machine constants for the honeywell 600/6000 series.
c
c     data rmach(1) / o402400000000 /
c     data rmach(2) / o376777777777 /
c     data rmach(3) / o714400000000 /
c     data rmach(4) / o716400000000 /
c     data rmach(5) / o776464202324 /
c
c     machine constants for the hp 2100
c
c     3 word double precision with ftn4
c
c     data small(1), small(2) / 40000b,       1 /
c     data large(1), large(2) / 77777b, 177776b /
c     data right(1), right(2) / 40000b,    325b /
c     data diver(1), diver(2) / 40000b,    327b /
c     data log10(1), log10(2) / 46420b,  46777b /
c
c     machine constants for the hp 2100
c     4 word double precision with ftn4
c
c     data small(1), small(2) / 40000b,       1 /
c     data large91), large(2) / 77777b, 177776b /
c     data right(1), right(2) / 40000b,    325b /
c     data diver(1), diver(2) / 40000b,    327b /
c     data log10(1), log10(2) / 46420b,  46777b /
c
c     machine constants for the ibm 360/370 series,
c     the xerox sigma 5/7/9, the sel systems 85/86  and
c     the perkin elmer (interdata) 7/32.
c
c     data rmach(1) / z00100000 /
c     data rmach(2) / z7fffffff /
c     data rmach(3) / z3b100000 /
c     data rmach(4) / z3c100000 /
c     data rmach(5) / z41134413 /
c
c     machine constants for the pdp-10 (ka or ki processor).
c
c     data rmach(1) / "000400000000 /
c     data rmach(2) / "377777777777 /
c     data rmach(3) / "146400000000 /
c     data rmach(4) / "147400000000 /
c     data rmach(5) / "177464202324 /
c
c     machine constants for pdp-11 fortran[s supporting
c     32-bit integers (expressed in integer and octal).
c
c     data small(1) /    8388608 /
c     data large(1) / 2147483647 /
c     data right(1) /  880803840 /
c     data diver(1) /  889192448 /
c     data log10(1) / 1067065499 /
c
c     data rmach(1) / o00040000000 /
c     data rmach(2) / o17777777777 /
c     data rmach(3) / o06440000000 /
c     data rmach(4) / o06500000000 /
c     data rmach(5) / o07746420233 /
c
c     machine constants for pdp-11 fortran[s supporting
c     16-bit integers  (expressed in integer and octal).
c
c     data small(1),small(2) /   128,     0 /
c     data large(1),large(2) / 32767,    -1 /
c     data right(1),right(2) / 13440,     0 /
c     data diver(1),diver(2) / 13568,     0 /
c     data log10(1),log10(2) / 16282,  8347 /
c
c     data small(1),small(2) / o000200, o000000 /
c     data large(1),large(2) / o077777, o177777 /
c     data right(1),right(2) / o032200, o000000 /
c     data diver(1),diver(2) / o032400, o000000 /
c     data log10(1),log10(2) / o037632, o020233 /
c
c     machine constants for the univac 1100 series.
c
c     data rmach(1) / o000400000000 /
c     data rmach(2) / o377777777777 /
c     data rmach(3) / o146400000000 /
c     data rmach(4) / o147400000000 /
c     data rmach(5) / o177464202324 /
c
c     machine constants for the vax 11/780
c
c     data small(1) / z00000080 /
c     data large(2) / zffff7fff /
c     data right(3) / z00003480 /
c     data diver(4) / z00003500 /
c     data log10(5) / z20983f9a /
c
c    machine constants for the ridge double precision

      data rmach(1) / 2.2250738585072027d-308/
      data rmach(2) / 1.7976931348623148d308/
      data rmach(3) / 1.110223d-16/
      data rmach(4) / 2.220446d-16/
      data rmach(5) / 3.010299d-01/
c
c***first executable statement  r1mach
c
      if (i .lt. 1  .or.  i .gt. 5)
     1   call xerror (25hr1mach -- i out of bounds,25,1,2)
c
      r1mach = rmach(i)
      return
c
      end
      subroutine r1mpyq(m,n,a,lda,v,w)
c***begin prologue  r1mpyq
c
c***changed to double precision
c
c***refer to snsq,snsqe
c     **********
c
c     subroutine r1mpyq
c
c     given an m by n matrix a, this subroutine computes a*q where
c     q is the product of 2*(n - 1) transformations
c
c           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
c
c     and gv(i), gw(i) are givens rotations in the (i,n) plane which
c     eliminate elements in the i-th and n-th planes, respectively.
c     q itself is not given, rather the information to recover the
c     gv, gw rotations is supplied.
c
c     the subroutine statement is
c
c       subroutine r1mpyq(m,n,a,lda,v,w)
c
c     where
c
c       m is a positive integer input variable set to the number
c         of rows of a.
c
c       n is a positive integer input variable set to the number
c         of columns of a.
c
c       a is an m by n array. on input a must contain the matrix
c         to be postmultiplied by the orthogonal matrix q
c         described above. on output a*q has replaced a.
c
c       lda is a positive integer input variable not less than m
c         which specifies the leading dimension of the array a.
c
c       v is an input array of length n. v(i) must contain the
c         information necessary to recover the givens rotation gv(i)
c         described above.
c
c       w is an input array of length n. w(i) must contain the
c         information necessary to recover the givens rotation gw(i)
c         described above.
c
c     subroutines called
c
c       fortran-supplied ... abs,sqrt
c
c     minpack. version of december 1978.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     apply the first set of givens rotations to a.
c
c***routines called  (none)
c***end prologue  r1mpyq
      implicit real*8 (a-h,o-z)
      integer m,n,lda
      real*8 a(lda,n),v(n),w(n)
      integer i,j,nmj,nm1
      real*8 cos,one,sin,temp
      data one /1.0d0/
c***first executable statement  r1mpyq
      nm1 = n - 1
      if (nm1 .lt. 1) go to 50
      do 20 nmj = 1, nm1
	 j = n - nmj
	 if (abs(v(j)) .gt. one) cos = one/v(j)
	 if (abs(v(j)) .gt. one) sin = sqrt(one-cos**2)
	 if (abs(v(j)) .le. one) sin = v(j)
	 if (abs(v(j)) .le. one) cos = sqrt(one-sin**2)
	 do 10 i = 1, m
	    temp = cos*a(i,j) - sin*a(i,n)
	    a(i,n) = sin*a(i,j) + cos*a(i,n)
	    a(i,j) = temp
   10       continue
   20    continue
c
c     apply the second set of givens rotations to a.
c
      do 40 j = 1, nm1
	 if (abs(w(j)) .gt. one) cos = one/w(j)
	 if (abs(w(j)) .gt. one) sin = sqrt(one-cos**2)
	 if (abs(w(j)) .le. one) sin = w(j)
	 if (abs(w(j)) .le. one) cos = sqrt(one-sin**2)
	 do 30 i = 1, m
	    temp = cos*a(i,j) + sin*a(i,n)
	    a(i,n) = -sin*a(i,j) + cos*a(i,n)
	    a(i,j) = temp
   30       continue
   40    continue
   50 continue
      return
c
c     last card of subroutine r1mpyq.
c
      end
      subroutine r1updt(m,n,s,ls,u,v,w,sing)
c***begin prologue  r1updt
c
c***changed to double precision
c
c***refer to  snsqe,snsq
c     **********
c
c     subroutine r1updt
c
c     given an m by n lower trapezoidal matrix s, an m-vector u,
c     and an n-vector v, the problem is to determine an
c     orthogonal matrix q such that
c
c                   t
c           (s + u*v )*q
c
c     is again lower trapezoidal.
c
c     this subroutine determines q as the product of 2*(n - 1)
c     transformations
c
c           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
c
c     where gv(i), gw(i) are givens rotations in the (i,n) plane
c     which eliminate elements in the i-th and n-th planes,
c     respectively. q itself is not accumulated, rather the
c     information to recover the gv, gw rotations is returned.
c
c     the subroutine statement is
c
c       subroutine r1updt(m,n,s,ls,u,v,w,sing)
c
c     where
c
c       m is a positive integer input variable set to the number
c         of rows of s.
c
c       n is a positive integer input variable set to the number
c         of columns of s. n must not exceed m.
c
c       s is an array of length ls. on input s must contain the lower
c         trapezoidal matrix s stored by columns. on output s contains
c         the lower trapezoidal matrix produced as described above.
c
c       ls is a positive integer input variable not less than
c         (n*(2*m-n+1))/2.
c
c       u is an input array of length m which must contain the
c         vector u.
c
c       v is an array of length n. on input v must contain the vector
c         v. on output v(i) contains the information necessary to
c         recover the givens rotation gv(i) described above.
c
c       w is an output array of length m. w(i) contains information
c         necessary to recover the givens rotation gw(i) described
c         above.
c
c       sing is a logical output variable. sing is set true if any
c         of the diagonal elements of the output s are zero. otherwise
c         v is set false.
c
c     subprograms called
c
c       minpack-supplied ... r1mach
c
c       fortran-supplied ... abs,sqrt
c
c     minpack. version of december 1978.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more,
c     john l. nazareth
c
c     **********
c
c     giant is the largest magnitude.
c
c***routines called  r1mach
c***end prologue  r1updt
      implicit real*8 (a-h,o-z)
      integer m,n,ls
      logical sing
      real*8 s(ls),u(m),v(n),w(m)
      integer i,j,jj,l,nmj,nm1
      real*8 cos,cotan,giant,one,p5,p25,sin,tan,tau,temp,zero
      real*8 r1mach
      data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
c***first executable statement  r1updt
      giant = r1mach(2)
c
c     initialize the diagonal element pointer.
c
      jj = (n*(2*m - n + 1))/2 - (m - n)
c
c     move the nontrivial part of the last column of s into w.
c
      l = jj
      do 10 i = n, m
	 w(i) = s(l)
	 l = l + 1
   10    continue
c
c     rotate the vector v into a multiple of the n-th unit vector
c     in such a way that a spike is introduced into w.
c
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 nmj = 1, nm1
	 j = n - nmj
	 jj = jj - (m - j + 1)
	 w(j) = zero
	 if (v(j) .eq. zero) go to 50
c
c        determine a givens rotation which eliminates the
c        j-th element of v.
c
	 if (abs(v(n)) .ge. abs(v(j))) go to 20
	    cotan = v(n)/v(j)
	    sin = p5/sqrt(p25+p25*cotan**2)
	    cos = sin*cotan
	    tau = one
	    if (abs(cos)*giant .gt. one) tau = one/cos
	    go to 30
   20    continue
	    tan = v(j)/v(n)
	    cos = p5/sqrt(p25+p25*tan**2)
	    sin = cos*tan
	    tau = sin
   30    continue
c
c        apply the transformation to v and store the information
c        necessary to recover the givens rotation.
c
	 v(n) = sin*v(j) + cos*v(n)
	 v(j) = tau
c
c        apply the transformation to s and extend the spike in w.
c
	 l = jj
	 do 40 i = j, m
	    temp = cos*s(l) - sin*w(i)
	    w(i) = sin*s(l) + cos*w(i)
	    s(l) = temp
	    l = l + 1
   40       continue
   50    continue
   60    continue
   70 continue
c
c     add the spike from the rank 1 update to w.
c
      do 80 i = 1, m
	 w(i) = w(i) + v(n)*u(i)
   80    continue
c
c     eliminate the spike.
c
      sing = .false.
      if (nm1 .lt. 1) go to 140
      do 130 j = 1, nm1
	 if (w(j) .eq. zero) go to 120
c
c        determine a givens rotation which eliminates the
c        j-th element of the spike.
c
	 if (abs(s(jj)) .ge. abs(w(j))) go to 90
	    cotan = s(jj)/w(j)
	    sin = p5/sqrt(p25+p25*cotan**2)
	    cos = sin*cotan
	    tau = one
	    if (abs(cos)*giant .gt. one) tau = one/cos
	    go to 100
   90    continue
	    tan = w(j)/s(jj)
	    cos = p5/sqrt(p25+p25*tan**2)
	    sin = cos*tan
	    tau = sin
  100    continue
c
c        apply the transformation to s and reduce the spike in w.
c
	 l = jj
	 do 110 i = j, m
	    temp = cos*s(l) + sin*w(i)
	    w(i) = -sin*s(l) + cos*w(i)
	    s(l) = temp
	    l = l + 1
  110       continue
c
c        store the information necessary to recover the
c        givens rotation.
c
	 w(j) = tau
  120    continue
c
c        test for zero diagonal elements in the output s.
c
	 if (s(jj) .eq. zero) sing = .true.
	 jj = jj + (m - j + 1)
  130    continue
  140 continue
c
c     move w back into the last column of the output s.
c
      l = jj
      do 150 i = n, m
	 s(l) = w(i)
	 l = l + 1
  150    continue
      if (s(jj) .eq. zero) sing = .true.
      return
c
c     last card of subroutine r1updt.
c
      end
C=================================================================
      subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
c***begin prologue  rwupdt
c***refer to  snlse,snls
c     **********
c
c     subroutine rwupdt
c
c     given an n by n upper triangular matrix r, this subroutine
c     computes the qr decomposition of the matrix formed when a row
c     is added to r. if the row is specified by the vector w, then
c     rwupdt determines an orthogonal matrix q such that when the
c     n+1 by n matrix composed of r augmented by w is premultiplied
c     by (q transpose), the resulting matrix is upper trapezoidal.
c     the orthogonal matrix q is the product of n transformations
c
c           g(1)*g(2)* ... *g(n)
c
c     where g(i) is a givens rotation in the (i,n+1) plane which
c     eliminates elements in the i-th plane. rwupdt also
c     computes the product (q transpose)*c where c is the
c     (n+1)-vector (b,alpha). q itself is not accumulated, rather
c     the information to recover the g rotations is supplied.
c
c     the subroutine statement is
c
c       subroutine rwupdt(n,r,ldr,w,b,alpha,cos,sin)
c
c     where
c
c       n is a positive integer input variable set to the order of r.
c
c       r is an n by n array. on input the upper triangular part of
c         r must contain the matrix to be updated. on output r
c         contains the updated triangular matrix.
c
c       ldr is a positive integer input variable not less than n
c         which specifies the leading dimension of the array r.
c
c       w is an input array of length n which must contain the row
c         vector to be added to r.
c
c       b is an array of length n. on input b must contain the
c         first n elements of the vector c. on output b contains
c         the first n elements of the vector (q transpose)*c.
c
c       alpha is a variable. on input alpha must contain the
c         (n+1)-st element of the vector c. on output alpha contains
c         the (n+1)-st element of the vector (q transpose)*c.
c
c       cos is an output array of length n which contains the
c         cosines of the transforming givens rotations.
c
c       sin is an output array of length n which contains the
c         sines of the transforming givens rotations.
c
c     subprograms called
c
c       fortran-supplied ... abs,sqrt
c
c     minpack. version of december 1978.
c     burton s. garbow, dudley v. goetschel, kenneth e. hillstrom,
c     jorge j. more
c
c     **********
c
c***routines called  (none)
c***end prologue  rwupdt
      integer n,ldr
      double precision alpha
      double precision r(ldr,n),w(n),b(n),cos(n),sin(n)
      integer i,j,jm1
      double precision cotan,one,p5,p25,rowj,tan,temp,zero
      data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
c***first executable statement  rwupdt
      do 60 j = 1, n
	 rowj = w(j)
	 jm1 = j - 1
c
c        apply the previous transformations to
c        r(i,j), i=1,2,...,j-1, and to w(j).
c
	 if (jm1 .lt. 1) go to 20
	 do 10 i = 1, jm1
	    temp = cos(i)*r(i,j) + sin(i)*rowj
	    rowj = -sin(i)*r(i,j) + cos(i)*rowj
	    r(i,j) = temp
   10       continue
   20    continue
c
c        determine a givens rotation which eliminates w(j).
c
	 cos(j) = one
	 sin(j) = zero
	 if (rowj .eq. zero) go to 50
	 if (abs(r(j,j)) .ge. abs(rowj)) go to 30
	    cotan = r(j,j)/rowj
	    sin(j) = p5/sqrt(p25+p25*cotan**2)
	    cos(j) = sin(j)*cotan
	    go to 40
   30    continue
	    tan = rowj/r(j,j)
	    cos(j) = p5/sqrt(p25+p25*tan**2)
	    sin(j) = cos(j)*tan
   40    continue
c
c        apply the current transformation to r(j,j), b(j), and alpha.
c
	 r(j,j) = cos(j)*r(j,j) + sin(j)*rowj
	 temp = cos(j)*b(j) + sin(j)*alpha
	 alpha = -sin(j)*b(j) + cos(j)*alpha
	 b(j) = temp
   50    continue
   60    continue
      return
c
c     last card of subroutine rwupdt.
c
      end
      subroutine s88fmt(n,ivalue,ifmt)
c***begin prologue  s88fmt
c
c***changed to double precision
c
c***refer to  xerror
c     abstract
c        s88fmt replaces ifmt(1), ... ,ifmt(n) with the
c        characters corresponding to the n least significant
c        digits of ivalue.
c
c     taken from the bell laboratories port library error handler
c     latest revision ---  7 june 1978
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  (none)
c***end prologue  s88fmt
c
      implicit real*8 (a-h,o-z)
      dimension ifmt(n),idigit(10)
      data idigit(1),idigit(2),idigit(3),idigit(4),idigit(5),
     1     idigit(6),idigit(7),idigit(8),idigit(9),idigit(10)
     2     /1h0,1h1,1h2,1h3,1h4,1h5,1h6,1h7,1h8,1h9/
c***first executable statement  s88fmt
      nt = n
      it = ivalue
   10    if (nt .eq. 0) return
	 index = mod(it,10)
	 ifmt(nt) = idigit(index+1)
	 it = it/10
	 nt = nt - 1
	 go to 10
      end
      real*8 function scasum(n,cx,incx)
c***begin prologue  scasum
c***revision date  811015   (yymmdd)
c***category no.  f1a
c***keywords  blas,vector,complex,sum
c***date written  october 1979
c
c***changed to double precision june 1985 by m.s.
c
c***author lawson c. (jpl),hanson r. (sla),
c                            kincaid d. (u texas), krogh f. (jpl)
c***purpose
c  sum of magnitudes of real and imaginary components of complex vector
c***description
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       cx  complex vector with n elements
c     incx  storage spacing between elements of cx
c
c     --output--
c   scasum  single precision result (zero if n.le.0)
c
c     returns sums of magnitudes of real and imaginary parts of
c     components of cx.  note that this is not the l1 norm of cx.
c     casum = sum from 0 to n-1 of abs(real(cx(1+i*incx))) +
c             abs(imag(cx(1+i*incx)))
c
c
c***references
c  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c   *basic linear algebra subprograms for fortran usage*,
c  algorithm no. 539, transactions on mathematical software,
c  volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  scasum
      implicit real*8 (a-h,o-z)
      complex*16 cx(1)
c***first executable statement  scasum
      scasum=0.
      if(n .le. 0) return
      ns = n*incx
	  do 10 i=1,ns,incx
	  scasum = scasum + abs(dble(cx(i))) + abs(dimag(cx(i)))
   10     continue
      return
      end
      subroutine snls(fcn,jac,iopt,m,n,x,fvec,fjac,ldfjac,ftol,xtol,
     *                 gtol,maxfev,epsfcn,diag,mode,factor,nprint,
     *                 info,nfev,njev,ipvt,qtf,wa1,wa2,wa3,wa4)
c***begin prologue  snls
c***revision date  811015   (yymmdd)
c***category no.  e2g1b1,e2g1b2
c***keywords nonlinear least squares,nonlinear data fitting,
c            levenberg-marquardt
c***date written  march 1980
c***author  hiebert k.l. (sla)
c***purpose
c       snls minimizes the sum of the squares of m nonlinear functions
c       in n variables by a modification of the levenberg-marquardt
c       algorithm.  this code is the combination of the minpack codes
c       (argonne) lmder, lmdif, and lmstr.
c***description
c
c
c 1. purpose.
c
c       the purpose of snls is to minimize the sum of the squares of m
c       nonlinear functions in n variables by a modification of the
c       levenberg-marquardt algorithm.  the user must provide a subrou-
c       tine which calculates the functions.  the user has the option
c       of how the jacobian will be supplied.  the user can supply the
c       full jacobian, or the rows of the jacobian (to avoid storing
c       the full jacobian), or let the code approximate the jacobian by
c       forward-differencing.   this code is the combination of the
c       minpack codes (argonne) lmder, lmdif, and lmstr.
c
c
c 2. subroutine and type statements.
c
c       subroutine snls(fcn,jac,iopt,m,n,x,fvec,fjac,ldfjac,ftol,xtol,
c      *                 gtol,maxfev,epsfcn,diag,mode,factor,nprint,info
c      *                 ,nfev,njev,ipvt,qtf,wa1,wa2,wa3,wa4)
c       integer iopt,m,n,ldfjac,maxfev,mode,nprint,info,nfev,njev
c       integer ipvt(n)
c       double precision ftol,xtol,gtol,epsfcn,factor
c       double precision x(n),fvec(m),fjac(ldfjac,n),diag(n),qtf(n),
c      *     wa1(n),wa2(n),wa3(n),wa4(m)
c
c
c 3. parameters.
c
c       parameters designated as input parameters must be specified on
c       entry to snls and are not changed on exit, while parameters
c       designated as output parameters need not be specified on entry
c       and are set to appropriate values on exit from snls.
c
c       fcn is the name of the user-supplied subroutine which calculates
c         the functions.  fcn must be declared in an external statement
c         in the user calling program, and should be written as follows.
c
c         subroutine fcn(m,n,x,fvec,iflag)
c         integer m,n,iflag
c         double precision x(n),fvec(m)
c         ----------
c         calculate the functions at x and
c         return this vector in fvec.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by fcn unless the
c         user wants to terminate execution of snls.  in this case set
c         iflag to a negative integer.
c
c       jac is the name of the user-supplied subroutine which calculates
c         the jacobian. if iopt=1 or 3, then jac must be declared in an
c         external statement in the user calling program, and should be
c         written as follows.
c
c         for iopt = 1
c         subroutine jac(m,n,x,fvec,fjac,ldfjac,iflag)
c         integer m,n,ldfjac,iflag
c         double precision x(n),fvec(m),fjac(ldfjac,n)
c         ----------
c         calculate the jacobian at x and return this
c         matrix in fjac.  fvec contains the function
c         values at x and should not be altered.
c         ----------
c         return
c         end
c
c         for iopt = 3
c         subroutine jac(m,n,x,fvec,fjrow,iflag)
c         integer m,n,iflag
c         double precision x(n),fvec(m),fjrow(n)
c         ----------
c         if iflag = i, then calculate the i-th row of
c         jacobian at x and return this vector in fjrow.
c         fvec contains the function values at x and
c         should not be altered.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by jac unless the
c         user wants to terminate execution of snls.  in this case set
c         iflag to a negative integer.
c
c         if iopt=2, jac can be ignored (treat it as a dummy argument).
c
c       iopt is an input variable which specifies how the jacobian will
c         be calculated.  if iopt=1 or 3, then the user must supply the
c         jacobian through the subroutine jac.  if iopt=1, then the user
c         supplies the full jacobian with each call to jac.  if iopt=3,
c         then the user supplies one row of the jacobian with each call.
c         (in this manner, storage can be saved because the full jacobia
c         jacobian is not stored.)   if iopt=2, then the code will
c         approximate the jacobian by forward-differencing.
c
c       m is a positive integer input variable set to the number of
c         functions.
c
c       n is a positive integer input variable set to the number of
c         variables.  n must not exceed m.
c
c       x is an array of length n.  on input x must contain an initial
c         estimate of the solution vector.  on output x contains the
c         final estimate of the solution vector.
c
c       fvec is an output array of length m which contains the functions
c         evaluated at the output x.
c
c       fjac is an output array.  for iopt=1 and 2, fjac is an m by n
c         array. for iopt=3, fjac is an n by n array.  the upper n by n
c         submatrix of fjac contains an upper triangular matrix r with
c         diagonal elements of nonincreasing magnitude such that
c
c                t     t           t
c               p *(jac *jac)*p = r *r,
c
c         where p is a permutation matrix and jac is the final calcu-
c         lated jacobian.  column j of p is column ipvt(j) (see below)
c         of the identity matrix.  the lower part of fjac contains
c         information generated during the computation of r.
c
c       ldfjac is a positive integer input variable which specifies
c         the leading dimension of the array fjac.  for iopt=1 and 2,
c         ldfjac must not be less than m.  for iopt=3, ldfjac must not
c         be less than n.
c
c       ftol is a nonnegative input variable.  termination occurs when
c         both the actual and predicted relative reductions in the sum
c         of squares are at most ftol.  therefore, ftol measures the
c         relative error desired in the sum of squares.  section 4 con-
c         tains more details about ftol.
c
c       xtol is a nonnegative input variable.  termination occurs when
c         the relative error between two consecutive iterates is at most
c         xtol.  therefore, xtol measures the relative error desired in
c         the approximate solution.  section 4 contains more details
c         about xtol.
c
c       gtol is a nonnegative input variable.  termination occurs when
c         the cosine of the angle between fvec and any column of the
c         jacobian is at most gtol in absolute value.  therefore, gtol
c         measures the orthogonality desired between the function vector
c         and the columns of the jacobian.  section 4 contains more
c         details about gtol.
c
c       maxfev is a positive integer input variable.  termination occurs
c         when the number of calls to fcn has reached maxfev.
c
c       epsfcn is an input variable used in determining a suitable step
c         for the forward-difference approximation.  this approximation
c         assumes that the relative errors in the functions are of the
c         order of epsfcn.  if epsfcn is less than the machine preci-
c         sion, it is assumed that the relative errors in the functions
c         are of the order of the machine precision.  if iopt=1 or 3,
c         then epsfcn can be ignored (treat it as a dummy argument).
c
c       diag is an array of length n.  if mode = 1 (see below), diag is
c         internally set.  if mode = 2, diag must contain positive
c         entries that serve as implicit (multiplicative) scale factors
c         for the variables.
c
c       mode is an integer input variable.  if mode = 1, the variables
c         will be scaled internally.  if mode = 2, the scaling is speci-
c         fied by the input diag.  other values of mode are equivalent
c         to mode = 1.
c
c       factor is a positive input variable used in determining the ini-
c         tial step bound.  this bound is set to the product of factor
c         and the euclidean norm of diag*x if nonzero, or else to factor
c         itself.  in most cases factor should lie in the interval
c         (.1,100.).  100. is a generally recommended value.
c
c       nprint is an integer input variable that enables controlled
c         printing of iterates if it is positive.  in this case, fcn is
c         called with iflag = 0 at the beginning of the first iteration
c         and every nprint iterations thereafter and immediately prior
c         to return, with x and fvec available for printing. appropriate
c         print statements must be added to fcn (see example) and
c         fvec should not be altered.  if nprint is not positive, no
c         special calls of fcn with iflag = 0 are made.
c
c       info is an integer output variable.  if the user has terminated
c         execution, info is set to the (negative) value of iflag.  see
c         description of fcn and jac. otherwise, info is set as follows.
c
c         info = 0  improper input parameters.
c
c         info = 1  both actual and predicted relative reductions in the
c                   sum of squares are at most ftol.
c
c         info = 2  relative error between two consecutive iterates is
c                   at most xtol.
c
c         info = 3  conditions for info = 1 and info = 2 both hold.
c
c         info = 4  the cosine of the angle between fvec and any column
c                   of the jacobian is at most gtol in absolute value.
c
c         info = 5  number of calls to fcn has reached maxfev.
c
c         info = 6  ftol is too small.  no further reduction in the sum
c                   of squares is possible.
c
c         info = 7  xtol is too small.  no further improvement in the
c                   approximate solution x is possible.
c
c         info = 8  gtol is too small.  fvec is orthogonal to the
c                   columns of the jacobian to machine precision.
c
c         sections 4 and 5 contain more details about info.
c
c       nfev is an integer output variable set to the number of calls to
c         fcn.
c
c       njev is an integer output variable set to the number of
c         evaluation of the full jacobian.  if iopt=1, only one call to
c         jac is required.  if iopt=3, then m calls to jac are required.
c         if iopt=2, then njev is set to zero.
c
c       ipvt is an integer output array of length n.  ipvt defines a
c         permutation matrix p such that jac*p = q*r, where jac is the
c         final calculated jacobian, q is orthogonal (not stored), and r
c         is upper triangular with diagonal elements of nonincreasing
c         magnitude.  column j of p is column ipvt(j) of the identity
c         matrix.
c
c       qtf is an output array of length n which contains the first n
c         elements of the vector (q transpose)*fvec.
c
c       wa1, wa2, and wa3 are work arrays of length n.
c
c       wa4 is a work array of length m.
c
c
c 4. successful completion.
c
c       the accuracy of snls is controlled by the convergence parame-
c       ters ftol, xtol, and gtol.  these parameters are used in tests
c       which make three types of comparisons between the approximation
c       x and a solution xsol.  snls terminates when any of the tests
c       is satisfied.  if any of the convergence parameters is less than
c       the machine precision (as defined by the function r1mach(4))
c       then snls only attempts to satisfy the test defined by the
c       machine precision.  further progress is not usually possible.
c
c       the tests assume that the functions are reasonably well behaved,
c       and, if the jacobian is supplied by the user, that the functions
c       and the jacobian are coded consistently.  if these conditions
c       are not satisfied, then snls may incorrectly indicate conver-
c       gence.  the coding of the jacobian can be checked by the
c       subroutine chkder. if the jacobian is coded correctly or iopt=2,
c       then the validity of the answer can be checked, for example, by
c       rerunning snls with tighter tolerances.
c
c       first convergence test.  if enorm(z) denotes the euclidean norm
c         of a vector z, then this test attempts to guarantee that
c
c               enorm(fvec) .le. (1+ftol)*enorm(fvecs),
c
c         where fvecs denotes the functions evaluated at xsol.  if this
c         condition is satisfied with ftol = 10**(-k), then the final
c         residual norm enorm(fvec) has k significant decimal digits and
c         info is set to 1 (or to 3 if the second test is also satis-
c         fied).  unless high precision solutions are required, the
c         recommended value for ftol is the square root of the machine
c         precision.
c
c       second convergence test.  if d is the diagonal matrix whose
c         entries are defined by the array diag, then this test attempts
c         to guarantee that
c
c               enorm(d*(x-xsol)) .le. xtol*enorm(d*xsol).
c
c         if this condition is satisfied with xtol = 10**(-k), then the
c         larger components of d*x have k significant decimal digits and
c         info is set to 2 (or to 3 if the first test is also satis-
c         fied).  there is a danger that the smaller components of d*x
c         may have large relative errors, but if mode = 1, then the
c         accuracy of the components of x is usually related to their
c         sensitivity.  unless high precision solutions are required,
c         the recommended value for xtol is the square root of the
c         machine precision.
c
c       third convergence test.  this test is satisfied when the cosine
c         of the angle between fvec and any column of the jacobian at x
c         is at most gtol in absolute value.  there is no clear rela-
c         tionship between this test and the accuracy of snls, and
c         furthermore, the test is equally well satisfied at other crit-
c         ical points, namely maximizers and saddle points.  therefore,
c         termination caused by this test (info = 4) should be examined
c         carefully.  the recommended value for gtol is zero.
c
c
c 5. unsuccessful completion.
c
c       unsuccessful termination of snls can be due to improper input
c       parameters, arithmetic interrupts, or an excessive number of
c       function evaluations.
c
c       improper input parameters.  info is set to 0 if iopt .lt. 1
c         or iopt .gt. 3, or n .le. 0, or m .lt. n, or for iopt=1 or 2
c         ldfjac .lt. m, or for iopt=3 ldfjac .lt. n, or ftol .lt. 0.e0,
c         or xtol .lt. 0.e0, or gtol .lt. 0.e0, or maxfev .le. 0, or
c         factor .le. 0.e0.
c
c       arithmetic interrupts.  if these interrupts occur in the fcn
c         subroutine during an early stage of the computation, they may
c         be caused by an unacceptable choice of x by snls.  in this
c         case, it may be possible to remedy the situation by rerunning
c         snls with a smaller value of factor.
c
c       excessive number of function evaluations.  a reasonable value
c         for maxfev is 100*(n+1) for iopt=1 or 3 and 200*(n+1) for
c         iopt=2.  if the number of calls to fcn reaches maxfev, then
c         this indicates that the routine is converging very slowly
c         as measured by the progress of fvec, and info is set to 5.
c         in this case, it may be helpful to restart snls with mode
c         set to 1.
c
c
c 6. characteristics of the algorithm.
c
c       snls is a modification of the levenberg-marquardt algorithm.
c       two of its main characteristics involve the proper use of
c       implicitly scaled variables (if mode = 1) and an optimal choice
c       for the correction.  the use of implicitly scaled variables
c       achieves scale invariance of snls and limits the size of the
c       correction in any direction where the functions are changing
c       rapidly.  the optimal choice of the correction guarantees (under
c       reasonable conditions) global convergence from starting points
c       far from the solution and a fast rate of convergence for
c       problems with small residuals.
c
c       timing.  the time required by snls to solve a given problem
c         depends on m and n, the behavior of the functions, the accu-
c         racy requested, and the starting point.  the number of arith-
c         metic operations needed by snls is about n**3 to process each
c         evaluation of the functions (call to fcn) and to process each
c         evaluation of the jacobian it takes m*n**2 for iopt=1 (one
c         call to jac), m*n**2 for iopt=2 (n calls to fcn) and
c         1.5*m*n**2 for iopt=3 (m calls to jac).  unless fcn and jac
c         can be evaluated quickly, the timing of snls will be
c         strongly influenced by the time spent in fcn and jac.
c
c       storage.  snls requires (m*n + 2*m + 6*n) for iopt=1 or 2 and
c         (n**2 + 2*m + 6*n) for iopt=3 single precision storage
c         locations and n integer storage locations, in addition to
c         the storage required by the program.  there are no internally
c         declared storage arrays.
c
c
c 7. example.
c
c       the problem is to determine the values of x(1), x(2), and x(3)
c       which provide the best fit (in the least squares sense) of
c
c             x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)),  i = 1, 15
c
c       to the data
c
c             y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,
c                  0.37,0.58,0.73,0.96,1.34,2.10,4.39),
c
c       where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)).  the
c       i-th component of fvec is thus defined by
c
c             y(i) - (x(1) + u(i)/(v(i)*x(2) + w(i)*x(3))).
c
c       **********
c
c       program test(input,output,tape6=output)
c c
c c     driver for snls example.
c c
c       integer j,iopt,m,n,ldfjac,maxfev,mode,nprint,info,nfev,njev,
c      *        nwrite
c       integer ipvt(3)
c       double precision ftol,xtol,gtol,factor,fnorm,epsfcn
c       double precision x(3),fvec(15),fjac(15,3),diag(3),qtf(3),
c      *     wa1(3),wa2(3),wa3(3),wa4(15)
c       double precision enorm,r1mach
c       external fcn,jac
c       data nwrite /6/
c c
c       iopt = 1
c       m = 15
c       n = 3
c c
c c     the following starting values provide a rough fit.
c c
c       x(1) = 1.e0
c       x(2) = 1.e0
c       x(3) = 1.e0
c c
c       ldfjac = 15
c c
c c     set ftol and xtol to the square root of the machine precision
c c     and gtol to zero. unless high precision solutions are
c c     required, these are the recommended settings.
c c
c       ftol = sqrt(r1mach(4))
c       xtol = sqrt(r1mach(4))
c       gtol = 0.e0
c c
c       maxfev = 400
c       epsfcn = 0.0
c       mode = 1
c       factor = 1.e2
c       nprint = 0
c c
c       call snls(fcn,jac,iopt,m,n,x,fvec,fjac,ldfjac,ftol,xtol,
c      *           gtol,maxfev,epsfcn,diag,mode,factor,nprint,
c      *           info,nfev,njev,ipvt,qtf,wa1,wa2,wa3,wa4)
c       fnorm = enorm(m,fvec)
c       write (nwrite,1000) fnorm,nfev,njev,info,(x(j),j=1,n)
c       stop
c  1000 format (5x,31h final l2 norm of the residuals,e15.7 //
c      *        5x,31h number of function evaluations,i10 //
c      *        5x,31h number of jacobian evaluations,i10 //
c      *        5x,15h exit parameter,16x,i10 //
c      *        5x,27h final approximate solution // 5x,3e15.7)
c       end
c       subroutine fcn(m,n,x,fvec,iflag)
c       integer m,n,iflag
c       double precision x(n),fvec(m)
c       integer i
c       double precision tmp1,tmp2,tmp3,tmp4
c       double precision y(15)
c       data y(1),y(2),y(3),y(4),y(5),y(6),y(7),y(8),
c      *     y(9),y(10),y(11),y(12),y(13),y(14),y(15)
c      *     /1.4e-1,1.8e-1,2.2e-1,2.5e-1,2.9e-1,3.2e-1,3.5e-1,3.9e-1,
c      *      3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0/
c c
c       if (iflag .ne. 0) go to 5
c c
c c     insert print statements here when nprint is positive.
c c
c       return
c     5 continue
c       do 10 i = 1, 15
c          tmp1 = i
c          tmp2 = 16 - i
c          tmp3 = tmp1
c          if (i .gt. 8) tmp3 = tmp2
c          fvec(i) = y(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3))
c    10    continue
c       return
c       end
c       subroutine jac(m,n,x,fvec,fjac,ldfjac,iflag)
c       integer m,n,ldfjac,iflag
c       double precision x(n),fvec(m),fjac(ldfjac,n)
c       integer i
c       double precision tmp1,tmp2,tmp3,tmp4
c c
c       do 10 i = 1, 15
c          tmp1 = i
c          tmp2 = 16 - i
c          tmp3 = tmp1
c          if (i .gt. 8) tmp3 = tmp2
c          tmp4 = (x(2)*tmp2 + x(3)*tmp3)**2
c          fjac(i,1) = -1.e0
c          fjac(i,2) = tmp1*tmp2/tmp4
c          fjac(i,3) = tmp1*tmp3/tmp4
c    10    continue
c       return
c       end
c
c
c       results obtained with different compilers or machines
c       may be slightly different.
c
c       final l2 norm of the residuals  0.9063596e-01
c
c       number of function evaluations         6
c
c       number of jacobian evaluations         5
c
c       exit parameter                         1
c
c       final approximate solution
c
c        0.8241058e-01  0.1133037e+01  0.2343695e+01
c
c
c       for iopt = 2, fcn would be the same and jac would not be
c         supplied.
c
c       for iopt = 3, fcn would be the same, fjac wolud be dimensioned
c         as fjac(3,3), ldfjac would be set to 3, and jac would be
c         written as follows.
c
c       subroutine jac(m,n,x,fvec,fjrow,iflag)
c       integer m,n,iflag
c       double precision x(n),fvec(m),fjrow(n)
c       integer i
c       double precision tmp1,tmp2,tmp3,tmp4
c
c       i = iflag
c          tmp1 = i
c          tmp2 = 16 - i
c          tmp3 = tmp1
c          if (i .gt. 8) tmp3 = tmp2
c          tmp4 = (x(2)*tmp2 + x(3)*tmp3)**2
c          fjrow(1) = -1.e0
c          fjrow(2) = tmp1*tmp2/tmp4
c          fjrow(3) = tmp1*tmp3/tmp4
c       return
c       end
c
c
c
c***references
c       more, jorge j.
c       the levenberg-marquardt algorithm, implementation and theory.
c       numerical analysis, g. a. watson, editor.
c       lecture notes in mathematics 630, springer-verlag, 1977.
c***routines called  enorm,fdjac2,lmpar,qrfac,rwupdt,r1mach,xerror
c***end prologue  snls
      integer iopt,m,n,ldfjac,maxfev,mode,nprint,info,nfev,njev
      integer ipvt(n)
      double precision ftol,xtol,gtol,factor,epsfcn
      double precision x(n),fvec(m),fjac(ldfjac,n),diag(n),
     &     qtf(n),wa1(n),wa2(n),
     *     wa3(n),wa4(m)
      logical sing
      external fcn
      integer i,iflag,iter,j,l
      double precision actred,delta,dirder,epsmch,fnorm,fnorm1,
     &     gnorm,one,par,
     *     pnorm,prered,p1,p5,p25,p75,p0001,ratio,sum,temp,temp1,
     *     temp2,xnorm,zero
      double precision r1mach,enorm
      data one,p1,p5,p25,p75,p0001,zero
     *     /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,1.0d-4,0.0d0/
c
c***first executable statement  snls
      epsmch = r1mach(4)
c
      info = 0
      iflag = 0
      nfev = 0
      njev = 0
c
c     check the input parameters for errors.
c
      if (iopt .lt. 1 .or. iopt .gt. 3 .or. n .le. 0 .or.
     *    m .lt. n .or. ldfjac .lt. n .or. ftol .lt. zero
     *    .or. xtol .lt. zero .or. gtol .lt. zero
     *    .or. maxfev .le. 0 .or. factor .le. zero) go to 300
      if (iopt .lt. 3 .and. ldfjac .lt. m) go to 300
      if (mode .ne. 2) go to 20
      do 10 j = 1, n
	 if (diag(j) .le. zero) go to 300
   10    continue
   20 continue
c
c     evaluate the function at the starting point
c     and calculate its norm.
c
      iflag = 1
      call fcn(m,n,x,fvec,iflag)
      nfev = 1
      if (iflag .lt. 0) go to 300
      fnorm = enorm(m,fvec)
c
c     initialize levenberg-marquardt parameter and iteration counter.
c
      par = zero
      iter = 1
c
c     beginning of the outer loop.
c
   30 continue
c
c        if requested, call fcn to enable printing of iterates.
c
	 if (nprint .le. 0) go to 40
	 iflag = 0
	 if (mod(iter-1,nprint) .eq. 0)
     *      call fcn(m,n,x,fvec,iflag)
	 if (iflag .lt. 0) go to 300
   40    continue
c
c        calculate the jacobian matrix.
c
      if (iopt .eq. 3) go to 475
c
c     store the full jacobian using m*n storage
c
      iflag = 2
      if (iopt .eq. 2) go to 410
c
c     the user supplies the jacobian
c
	 call jac(m,n,x,fvec,fjac,ldfjac,iflag)
	 njev = njev + 1
	 go to 420
c
c     the code approximates the jacobian
c
  410    call fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4)
	 nfev = nfev + n
  420    if (iflag .lt. 0) go to 300
c
c        compute the qr factorization of the jacobian.
c
	 call qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
c
c        form (q transpose)*fvec and store the first n components in
c        qtf.
c
	 do 430 i = 1, m
	    wa4(i) = fvec(i)
  430         continue
	 do 470 j = 1, n
	    if (fjac(j,j) .eq. zero) go to 460
	    sum = zero
	    do 440 i = j, m
	       sum = sum + fjac(i,j)*wa4(i)
  440          continue
	    temp = -sum/fjac(j,j)
	    do 450 i = j, m
	       wa4(i) = wa4(i) + fjac(i,j)*temp
  450          continue
  460       continue
	    fjac(j,j) = wa1(j)
	    qtf(j) = wa4(j)
  470       continue
	 go to 560
c
c        accumulate the jacobian by rows in order to save storage.
c        compute the qr factorization of the jacobian matrix
c        calculated one row at a time, while simultaneously
c        forming (q transpose)*fvec and storing the first
c        n components in qtf.
c
  475    do 490 j = 1, n
	    qtf(j) = zero
	    do 480 i = 1, n
	       fjac(i,j) = zero
  480          continue
  490        continue
	 iflag = 1
	 do 500 i = 1, m
	    call jac(m,n,x,fvec,wa3,iflag)
	    if (iflag .lt. 0) go to 300
	    temp = fvec(i)
	    call rwupdt(n,fjac,ldfjac,wa3,qtf,temp,wa1,wa2)
	    iflag = iflag + 1
  500       continue
	 njev = njev + 1
c
c        if the jacobian is rank deficient, call qrfac to
c        reorder its columns and update the components of qtf.
c
	 sing = .false.
	 do 510 j = 1, n
	    if (fjac(j,j) .eq. zero) sing = .true.
	    ipvt(j) = j
	    wa2(j) = enorm(j,fjac(1,j))
  510       continue
	 if (.not.sing) go to 560
	 call qrfac(n,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
	 do 550 j = 1, n
	    if (fjac(j,j) .eq. zero) go to 540
	    sum = zero
	    do 520 i = j, n
	       sum = sum + fjac(i,j)*qtf(i)
  520         continue
	    temp = -sum/fjac(j,j)
	    do 530 i = j, n
	       qtf(i) = qtf(i) + fjac(i,j)*temp
  530          continue
  540       continue
	    fjac(j,j) = wa1(j)
  550       continue
  560    continue
c
c        on the first iteration and if mode is 1, scale according
c        to the norms of the columns of the initial jacobian.
c
	 if (iter .ne. 1) go to 80
	 if (mode .eq. 2) go to 60
	 do 50 j = 1, n
	    diag(j) = wa2(j)
	    if (wa2(j) .eq. zero) diag(j) = one
   50       continue
   60    continue
c
c        on the first iteration, calculate the norm of the scaled x
c        and initialize the step bound delta.
c
	 do 70 j = 1, n
	    wa3(j) = diag(j)*x(j)
   70       continue
	 xnorm = enorm(n,wa3)
	 delta = factor*xnorm
	 if (delta .eq. zero) delta = factor
   80    continue
c
c        compute the norm of the scaled gradient.
c
	 gnorm = zero
	 if (fnorm .eq. zero) go to 170
	 do 160 j = 1, n
	    l = ipvt(j)
	    if (wa2(l) .eq. zero) go to 150
	    sum = zero
	    do 140 i = 1, j
	       sum = sum + fjac(i,j)*(qtf(i)/fnorm)
  140          continue
	    gnorm = dmax1(gnorm,abs(sum/wa2(l)))
  150       continue
  160       continue
  170    continue
c
c        test for convergence of the gradient norm.
c
	 if (gnorm .le. gtol) info = 4
	 if (info .ne. 0) go to 300
c
c        rescale if necessary.
c
	 if (mode .eq. 2) go to 190
	 do 180 j = 1, n
	    diag(j) = dmax1(diag(j),wa2(j))
  180       continue
  190    continue
c
c        beginning of the inner loop.
c
  200    continue
c
c           determine the levenberg-marquardt parameter.
c
	    call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,
     *                 wa3,wa4)
c
c           store the direction p and x + p. calculate the norm of p.
c
	    do 210 j = 1, n
	       wa1(j) = -wa1(j)
	       wa2(j) = x(j) + wa1(j)
	       wa3(j) = diag(j)*wa1(j)
  210          continue
	    pnorm = enorm(n,wa3)
c
c           on the first iteration, adjust the initial step bound.
c
	    if (iter .eq. 1) delta = dmin1(delta,pnorm)
c
c           evaluate the function at x + p and calculate its norm.
c
	    iflag = 1
	    call fcn(m,n,wa2,wa4,iflag)
	    nfev = nfev + 1
	    if (iflag .lt. 0) go to 300
	    fnorm1 = enorm(m,wa4)
c
c           compute the scaled actual reduction.
c
	    actred = -one
	    if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
c
c           compute the scaled predicted reduction and
c           the scaled directional derivative.
c
	    do 230 j = 1, n
	       wa3(j) = zero
	       l = ipvt(j)
	       temp = wa1(l)
	       do 220 i = 1, j
		  wa3(i) = wa3(i) + fjac(i,j)*temp
  220             continue
  230          continue
	    temp1 = enorm(n,wa3)/fnorm
	    temp2 = (sqrt(par)*pnorm)/fnorm
	    prered = temp1**2 + temp2**2/p5
	    dirder = -(temp1**2 + temp2**2)
c
c           compute the ratio of the actual to the predicted
c           reduction.
c
	    ratio = zero
	    if (prered .ne. zero) ratio = actred/prered
c
c           update the step bound.
c
	    if (ratio .gt. p25) go to 240
	       if (actred .ge. zero) temp = p5
	       if (actred .lt. zero)
     *            temp = p5*dirder/(dirder + p5*actred)
	       if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1
	       delta = temp*dmin1(delta,pnorm/p1)
	       par = par/temp
	       go to 260
  240       continue
	       if (par .ne. zero .and. ratio .lt. p75) go to 250
	       delta = pnorm/p5
	       par = p5*par
  250          continue
  260       continue
c
c           test for successful iteration.
c
	    if (ratio .lt. p0001) go to 290
c
c           successful iteration. update x, fvec, and their norms.
c
	    do 270 j = 1, n
	       x(j) = wa2(j)
	       wa2(j) = diag(j)*x(j)
  270          continue
	    do 280 i = 1, m
	       fvec(i) = wa4(i)
  280          continue
	    xnorm = enorm(n,wa2)
	    fnorm = fnorm1
	    iter = iter + 1
  290       continue
c
c           tests for convergence.
c
	    if (abs(actred) .le. ftol .and. prered .le. ftol
     *          .and. p5*ratio .le. one) info = 1
	    if (delta .le. xtol*xnorm) info = 2
	    if (abs(actred) .le. ftol .and. prered .le. ftol
     *          .and. p5*ratio .le. one .and. info .eq. 2) info = 3
	    if (info .ne. 0) go to 300
c
c           tests for termination and stringent tolerances.
c
	    if (nfev .ge. maxfev) info = 5
	    if (abs(actred) .le. epsmch .and. prered .le. epsmch
     *          .and. p5*ratio .le. one) info = 6
	    if (delta .le. epsmch*xnorm) info = 7
	    if (gnorm .le. epsmch) info = 8
	    if (info .ne. 0) go to 300
c
c           end of the inner loop. repeat if iteration unsuccessful.
c
	    if (ratio .lt. p0001) go to 200
c
c        end of the outer loop.
c
	 go to 30
  300 continue
c
c     termination, either normal or user imposed.
c
      if (iflag .lt. 0) info = iflag
      iflag = 0
      if (nprint .gt. 0) call fcn(m,n,x,fvec,iflag)
      if (info .lt. 0) call xerror(63hsnls   -- execution terminated bec
     1ause user set iflag negative.,63,1,1)
      if (info .eq. 0) call xerror(34hsnls   -- invalid input parameter.
     1,34,2,1)
      if (info .eq. 4) call xerror(70hsnls   -- third convergence condit
     1ion, check results before accepting.,70,1,1)
      if (info .eq. 5) call xerror(40hsnls   -- too many function evalua
     1tions.,40,9,1)
      if (info .ge. 6) call xerror(64hsnls   -- tolerances too small, no
     1 further improvement possible.,64,3,1)
      return
c
c     last card of subroutine snls.
c
      end
      subroutine snlse(fcn,jac,iopt,m,n,x,fvec,tol,nprint,info,iw,
     *                  wa,lwa)
c***begin prologue  snlse
c***revision date  811015   (yymmdd)
c***category no.  e2g1b1,e2g1b2
c***keywords easy-to-use,nonlinear least squares,nonlinear data
c            fitting,levenberg-marquardt
c***date written  march 1980
c***author  hiebert k.l. (sla)
c***purpose
c       snlse is the easy-to-use version of snls which minimizes the sum
c       of the squares of m nonlinear functions in n variables by a
c       modification of the levenberg-marquardt algorithm.  this code is
c       the combination of the minpack codes (argonne) lmder1, lmdif1,
c       and lmstr.
c***description
c
c
c 1. purpose.
c
c       the purpose of snlse is to minimize the sum of the squares of m
c       nonlinear functions in n variables by a modification of the
c       levenberg-marquardt algorithm.  this is done by using the more
c       general least-squares solver snls.  the user must provide a sub-
c       routine which calculates the functions.  the user has the option
c       of how the jacobian will be supplied.  the user can supply the
c       full jacobian, or the rows of the jacobian (to avoid storing
c       the full jacobian), or let the code approximate the jacobian by
c       forward-differencing.   this code is the combination of the
c       minpack codes (argonne) lmder1, lmdif1, and lmstr1.
c
c
c 2. subroutine and type statements.
c
c       subroutine snlse(fcn,jac,iopt,m,n,x,fvec,tol,nprint,
c      *                  info,iw,wa,lwa)
c       integer iopt,m,n,nprint,info,lwa
c       integer iw(n)
c       double precision tol
c       double precision x(n),fvec(m),wa(lwa)
c       external fcn,jac
c
c
c 3. parameters.
c
c       parameters designated as input parameters must be specified on
c       entry to snlse and are not changed on exit, while parameters
c       designated as output parameters need not be specified on entry
c       and are set to appropriate values on exit from snlse.
c
c       fcn is the name of the user-supplied subroutine which calculates
c         the functions.  fcn must be declared in an external statement
c         in the user calling program, and should be written as follows.
c
c         subroutine fcn(m,n,x,fvec,iflag)
c         integer m,n,iflag
c         double precision x(n),fvec(m)
c         ----------
c         calculate the functions at x and
c         return this vector in fvec.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by fcn unless the
c         user wants to terminate execution of snlse.  in this case set
c         iflag to a negative integer.
c
c       jac is the name of the user-supplied subroutine which calculates
c         the jacobian. if iopt=1 or 3, then jac must be declared in an
c         external statement in the user calling program, and should be
c         written as follows.
c
c         for iopt = 1
c         subroutine jac(m,n,x,fvec,fjac,ldfjac,iflag)
c         integer m,n,ldfjac,iflag
c         double precision x(n),fvec(m),fjac(ldfjac,n)
c         ----------
c         calculate the jacobian at x and return this
c         matrix in fjac.  fvec contains the function
c         values at x and should not be altered.
c         ----------
c         return
c         end
c
c         for iopt = 3
c         subroutine jac(m,n,x,fvec,fjrow,iflag)
c         integer m,n,iflag
c         double precision x(n),fvec(m),fjrow(n)
c         ----------
c         if iflag = i, then calculate the i-th row of
c         jacobian at x and return this vector in fjrow.
c         fvec contains the function values at x and
c         should not be altered.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by jac unless the
c         user wants to terminate execution of snlse.  in this case set
c         iflag to a negative integer.
c
c         if iopt=2, jac can be ignored (treat it as a dummy argument).
c
c       iopt is an input variable which specifies how the jacobian will
c         be calculated.  if iopt=1 or 3, then the user must supply the
c         jacobian through the subroutine jac.  if iopt=1, then the user
c         supplies the full jacobian with each call to jac.  if iopt=3,
c         then the user supplies one row of the jacobian with each call.
c         (in this manner, storage can be saved because the full
c         jacobian is not stored.)   if iopt=2, then the code will
c         approximate the jacobian by forward-differencing.
c
c       m is a positive integer input variable set to the number of
c         functions.
c
c       n is a positive integer input variable set to the number of
c         variables.  n must not exceed m.
c
c       x is an array of length n.  on input x must contain an initial
c         estimate of the solution vector.  on output x contains the
c         final estimate of the solution vector.
c
c       fvec is an output array of length m which contains the functions
c         evaluated at the output x.
c
c       tol is a nonnegative input variable.  termination occurs when
c         the algorithm estimates either that the relative error in the
c         sum of squares is at most tol or that the relative error
c         between x and the solution is at most tol.  section 4 contains
c         more details about tol.
c
c       nprint is an integer input variable that enables controlled
c         printing of iterates if it is positive.  in this case, fcn is
c         called with iflag = 0 at the beginning of the first iteration
c         and every nprint iterations thereafter and immediately prior
c         to return, with x and fvec available for printing. appropriate
c         print statements must be added to fcn (see example) and
c         fvec should not be altered.  if nprint is not positive, no
c         special calls of fcn with iflag = 0 are made.
c
c       info is an integer output variable.  if the user has terminated
c         execution, info is set to the (negative) value of iflag.  see
c         description of fcn and jac. otherwise, info is set as follows.
c
c         info = 0  improper input parameters.
c
c         info = 1  algorithm estimates that the relative error in the
c                   sum of squares is at most tol.
c
c         info = 2  algorithm estimates that the relative error between
c                   x and the solution is at most tol.
c
c         info = 3  conditions for info = 1 and info = 2 both hold.
c
c         info = 4  fvec is orthogonal to the columns of the jacobian to
c                   machine precision.
c
c         info = 5  number of calls to fcn has reached 100*(n+1)
c                   for iopt=1 or 3 or 200*(n+1) for iopt=2.
c
c         info = 6  tol is too small.  no further reduction in the sum
c                   of squares is possible.
c
c         info = 7  tol is too small.  no further improvement in the
c                   approximate solution x is possible.
c
c         sections 4 and 5 contain more details about info.
c
c       iw is an integer work array of length n.
c
c       wa is a work array of length lwa.
c
c       lwa is a positive integer input variable not less than
c         n*(m+5)+m for iopt=1 and 2 or n*(n+5)+m for iopt=3.
c
c
c 4. successful completion.
c
c       the accuracy of snlse is controlled by the convergence parame-
c       ter tol.  this parameter is used in tests which make three types
c       of comparisons between the approximation x and a solution xsol.
c       snlse terminates when any of the tests is satisfied.  if tol is
c       less than the machine precision (as defined by the function
c       r1mach(4)), then snlse only attempts to satisfy the test
c       defined by the machine precision.  further progress is not usu-
c       ally possible.  unless high precision solutions are required,
c       the recommended value for tol is the square root of the machine
c       precision.
c
c       the tests assume that the functions are reasonably well behaved,
c       and, if the jacobian is supplied by the user, that the functions
c       and the jacobian are coded consistently.  if these conditions
c       are not satisfied, then snlse may incorrectly indicate conver-
c       gence.  the coding of the jacobian can be checked by the
c       subroutine chkder. if the jacobian is coded correctly or iopt=2,
c       then the validity of the answer can be checked, for example, by
c       rerunning snlse with tighter tolerances.
c
c       first convergence test.  if enorm(z) denotes the euclidean norm
c         of a vector z, then this test attempts to guarantee that
c
c               enorm(fvec) .le. (1+tol)*enorm(fvecs),
c
c         where fvecs denotes the functions evaluated at xsol.  if this
c         condition is satisfied with tol = 10**(-k), then the final
c         residual norm enorm(fvec) has k significant decimal digits and
c         info is set to 1 (or to 3 if the second test is also satis-
c         fied).
c
c       second convergence test.  if d is a diagonal matrix (implicitly
c         generated by snlse) whose entries contain scale factors for
c         the variables, then this test attempts to guarantee that
c
c               enorm(d*(x-xsol)) .le. tol*enorm(d*xsol).
c
c         if this condition is satisfied with tol = 10**(-k), then the
c         larger components of d*x have k significant decimal digits and
c         info is set to 2 (or to 3 if the first test is also satis-
c         fied).  there is a danger that the smaller components of d*x
c         may have large relative errors, but the choice of d is such
c         that the accuracy of the components of x is usually related to
c         their sensitivity.
c
c       third convergence test.  this test is satisfied when fvec is
c         orthogonal to the columns of the jacobian to machine preci-
c         sion.  there is no clear relationship between this test and
c         the accuracy of snlse, and furthermore, the test is equally
c         well satisfied at other critical points, namely maximizers and
c         saddle points.  therefore, termination caused by this test
c         (info = 4) should be examined carefully.
c
c
c 5. unsuccessful completion.
c
c       unsuccessful termination of snlse can be due to improper input
c       parameters, arithmetic interrupts, or an excessive number of
c       function evaluations.
c
c       improper input parameters.  info is set to 0 if iopt .lt. 1
c         or iopt .gt. 3, or n .le. 0, or m .lt. n, or tol .lt. 0.e0,
c         or for iopt=1 or 2 lwa .lt. n*(m+5)+m, or for iopt=3
c         lwa .lt. n*(n+5)+m.
c
c       arithmetic interrupts.  if these interrupts occur in the fcn
c         subroutine during an early stage of the computation, they may
c         be caused by an unacceptable choice of x by snlse.  in this
c         case, it may be possible to remedy the situation by not evalu-
c         ating the functions here, but instead setting the components
c         of fvec to numbers that exceed those in the initial fvec.
c
c       excessive number of function evaluations.  if the number of
c         calls to fcn reaches 100*(n+1) for iopt=1 or 3 or 200*(n+1)
c         for iopt=2, then this indicates that the routine is converging
c         very slowly as measured by the progress of fvec, and info is
c         set to 5.  in this case, it may be helpful to restart snlse
c         thereby forcing it to disregaad old (and possibly harmful)
c         information.
c
c
c 6. characteristics of the algorithm.
c
c       snlse is a modification of the levenberg-marquardt algorithm.
c       two of its main characteristics involve the proper use of
c       implicitly scaled variables and an optimal choice for the cor-
c       rection.  the use of implicitly scaled variables achieves scale
c       invariance of snlse and limits the size of the correction in
c       any direction where the functions are changing rapidly.  the
c       optimal choice of the correction guarantees (under reasonable
c       conditions) global convergence from starting points far from the
c       solution and a fast rate of convergence for problems with small
c       residuals.
c
c       timing.  the time required by snlse to solve a given problem
c         depends on m and n, the behavior of the functions, the accu-
c         racy requested, and the starting point.  the number of arith-
c         metic operations needed by snlse is about n**3 to process each
c         evaluation of the functions (call to fcn) and to process each
c         evaluation of the jacobian snlse takes m*n**2 for iopt=1 (one
c         call to jac), m*n**2 for iopt=2 (n calls to fcn) and
c         1.5*m*n**2 for iopt=3 (m calls to jac).  unless fcn and jac
c         can be evaluated quickly, the timing of snlse will be
c         strongly influenced by the time spent in fcn and jac.
c
c       storage.  snlse requires (m*n + 2*m + 6*n) for iopt=1 or 2 and
c         (n**2 + 2*m + 6*n) for iopt=3 single precision storage
c         locations and n integer storage locations, in addition to
c         the storage required by the program.  there are no internally
c         declared storage arrays.
c
c
c
c 7. example.
c
c       the problem is to determine the values of x(1), x(2), and x(3)
c       which provide the best fit (in the least squares sense) of
c
c             x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)),  i = 1, 15
c
c       to the data
c
c             y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,
c                  0.37,0.58,0.73,0.96,1.34,2.10,4.39),
c
c       where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)).  the
c       i-th component of fvec is thus defined by
c
c             y(i) - (x(1) + u(i)/(v(i)*x(2) + w(i)*x(3))).
c
c       **********
c
c       program test(input,output,tape6=output)
c c
c c     driver for snlse example.
c c
c       integer j,iopt,m,n,nprint,info,lwa,nwrite
c       integer iw(3)
c       double precision tol,fnorm
c       double precision x(3),fvec(15),wa(75)
c       double precision enorm,r1mach
c       external fcn,jac
c       data nwrite /6/
c c
c       iopt = 1
c       m = 15
c       n = 3
c c
c c     the following starting values provide a rough fit.
c c
c       x(1) = 1.e0
c       x(2) = 1.e0
c       x(3) = 1.e0
c c
c       lwa = 75
c       nprint = 0
c c
c c     set tol to the square root of the machine precision.
c c     unless high precision solutions are required,
c c     this is the recommended setting.
c c
c       tol = sqrt(r1mach(4))
c c
c       call snlse(fcn,jac,iopt,m,n,x,fvec,tol,nprint,
c      *            info,iw,wa,lwa)
c       fnorm = enorm(m,fvec)
c       write (nwrite,1000) fnorm,info,(x(j),j=1,n)
c       stop
c  1000 format (5x,31h final l2 norm of the residuals,e15.7 //
c      *        5x,15h exit parameter,16x,i10 //
c      *        5x,27h final approximate solution // 5x,3e15.7)
c       end
c       subroutine fcn(m,n,x,fvec,iflag)
c       integer m,n,iflag
c       double precision x(n),fvec(m)
c       integer i
c       double precision tmp1,tmp2,tmp3,tmp4
c       double precision y(15)
c       data y(1),y(2),y(3),y(4),y(5),y(6),y(7),y(8),
c      *     y(9),y(10),y(11),y(12),y(13),y(14),y(15)
c      *     /1.4e-1,1.8e-1,2.2e-1,2.5e-1,2.9e-1,3.2e-1,3.5e-1,3.9e-1,
c      *      3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0/
c c
c       do 10 i = 1, 15
c          tmp1 = i
c          tmp2 = 16 - i
c          tmp3 = tmp1
c          if (i .gt. 8) tmp3 = tmp2
c          fvec(i) = y(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3))
c    10    continue
c       return
c       end
c       subroutine jac(m,n,x,fvec,fjac,ldfjac,iflag)
c       integer m,n,ldfjac,iflag
c       double precision x(n),fvec(m),fjac(ldfjac,n)
c       integer i
c       double precision tmp1,tmp2,tmp3,tmp4
c c
c       do 10 i = 1, 15
c          tmp1 = i
c          tmp2 = 16 - i
c          tmp3 = tmp1
c          if (i .gt. 8) tmp3 = tmp2
c          tmp4 = (x(2)*tmp2 + x(3)*tmp3)**2
c          fjac(i,1) = -1.e0
c          fjac(i,2) = tmp1*tmp2/tmp4
c          fjac(i,3) = tmp1*tmp3/tmp4
c    10    continue
c       return
c       end
c
c       results obtained with different compilers or machines
c       may be slightly different.
c
c       final l2 norm of the residuals  0.9063596e-01
c
c       exit parameter                         1
c
c       final approximate solution
c
c        0.8241058e-01  0.1133037e+01  0.2343695e+01
c
c
c       for iopt = 2, fcn would be the same and jac would not be
c         supplied.
c
c       for iopt = 3, fcn would be the same, fjac wolud be dimensioned
c         as fjac(3,3), ldfjac would be set to 3, and jac would be
c         written as follows.
c
c       subroutine jac(m,n,x,fvec,fjrow,iflag)
c       integer m,n,iflag
c       double precision x(n),fvec(m),fjrow(n)
c       integer i
c       double precision tmp1,tmp2,tmp3,tmp4
c
c       i = iflag
c          tmp1 = i
c          tmp2 = 16 - i
c          tmp3 = tmp1
c          if (i .gt. 8) tmp3 = tmp2
c          tmp4 = (x(2)*tmp2 + x(3)*tmp3)**2
c          fjrow(1) = -1.e0
c          fjrow(2) = tmp1*tmp2/tmp4
c          fjrow(3) = tmp1*tmp3/tmp4
c       return
c       end
c
c
c***references
c       more, jorge j.
c       the levenberg-marquardt algorithm, implementation and theory.
c       numerical analysis, g. a. watson, editor.
c       lecture notes in mathematics 630, springer-verlag, 1977.
c***routines called  snls,xerror
c***end prologue  snlse
      integer m,n,nprint,info,lwa,iopt
      integer iw(n)
      double precision tol
      double precision x(n),fvec(m),wa(lwa)
      external fcn,jac
      integer maxfev,mode,nfev,njev
      double precision factor,ftol,gtol,xtol,zero,epsfcn
      data factor,zero /1.0d2,0.0d0/
c***first executable statement  snlse
      info = 0
c
c     check the input parameters for errors.
c
      if (iopt .lt. 1 .or. iopt .gt. 3 .or.
     *    n .le. 0 .or. m .lt. n .or. tol .lt. zero
     *    .or. lwa .lt. n*(n+5) + m) go to 10
      if (iopt .lt. 3 .and. lwa .lt. n*(m+5) + m) go to 10
c
c     call snls.
c
      maxfev = 100*(n + 1)
      if (iopt .eq. 2) maxfev = 2*maxfev
      ftol = tol
      xtol = tol
      gtol = zero
      epsfcn = zero
      mode = 1
      index = 5*n+m
      call snls(fcn,jac,iopt,m,n,x,fvec,wa(index+1),m,ftol,xtol,gtol,
     *           maxfev,epsfcn,wa(1),mode,factor,nprint,info,nfev,njev,
     *           iw,wa(n+1),wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
      if (info .eq. 8) info = 4
   10 continue
      if (info .eq. 0) call xerror(34hsnlse  -- invalid input parameter.
     1,34,2,1)
      return
c
c     last card of subroutine snlse.
c
      end
      subroutine snsq(fcn,jac,iopt,n,x,fvec,fjac,ldfjac,xtol,maxfev,
     *                 ml,mu,epsfcn,diag,mode,factor,nprint,info,nfev,
     *                 njev,r,lr,qtf,wa1,wa2,wa3,wa4)
c***begin prologue  snsq
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  c5a,c5b
c***keywords nonlinear square system, zero, powell hybrid method
c***date written  march 1980
c***author  hiebert k.l. (sla)
c***purpose
c       snsq finds to find a zero of a system of n nonlinear functions
c       in n variables by a modification of the powell hybrid method.
c       this code is the combination of the minpack codes (argonne)
c       hybrd and hybrdj.
c***description
c
c
c 1. purpose.
c
c       the purpose of snsq is to find a zero of a system of n non-
c       linear functions in n variables by a modification of the powell
c       hybrid method.  the user must provide a subroutine which calcu-
c       lates the functions.  the user has the option of either to
c       provide a subroutine which calculates the jacobian or to let the
c       code calculate it by a forward-difference approximation.
c       this code is the combination of the minpack codes (argonne)
c       hybrd and hybrdj.
c
c
c 2. subroutine and type statements.
c
c       subroutine snsq(fcn,jac,iopt,n,x,fvec,fjac,ldfjac,xtol,maxfev,
c      *                 ml,mu,epsfcn,diag,mode,factor,nprint,info,nfev,
c      *                 njev,r,lr,qtf,wa1,wa2,wa3,wa4)
c       integer iopt,n,maxfev,ml,mu,mode,nprint,info,nfev,ldfjac,njev,lr
c       real xtol,epsfcn,factor
c       real x(n),fvec(n),diag(n),fjac(ldfjac,n),r(lr),qtf(n),
c      *     wa1(n),wa2(n),wa3(n),wa4(n)
c       external fcn,jac
c
c
c 3. parameters.
c
c       parameters designated as input parameters must be specified on
c       entry to snsq and are not changed on exit, while parameters
c       designated as output parameters need not be specified on entry
c       and are set to appropriate values on exit from snsq.
c
c       fcn is the name of the user-supplied subroutine which calculates
c         the functions.  fcn must be declared in an external statement
c         in the user calling program, and should be written as follows.
c
c         subroutine fcn(n,x,fvec,iflag)
c         integer n,iflag
c         real x(n),fvec(n)
c         ----------
c         calculate the functions at x and
c         return this vector in fvec.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by fcn unless the
c         user wants to terminate execution of snsq.  in this case set
c         iflag to a negative integer.
c
c       jac is the name of the user-supplied subroutine which calculates
c         the jacobian.  if iopt=1, then jac must be declared in an
c         external statement in the user calling program, and should be
c         written as follows.
c
c         subroutine jac(n,x,fvec,fjac,ldfjac,iflag)
c         integer n,ldfjac,iflag
c         real x(n),fvec(n),fjac(ldfjac,n)
c         ----------
c         calculate the jacobian at x and return this
c         matrix in fjac.  fvec contains the function
c         values at x and should not be altered.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by jac unless the
c         user wants to terminate execution of snsq.  in this case set
c         iflag to a negative integer.
c
c         if iopt=2, jac can be ignored (treat it as a dummy argument).
c
c       iopt is an input variable which specifies how the jacobian will
c         be calculated.  if iopt=1, then the user must supply the
c         jacobian through the subroutine jac.  if iopt=2, then the
c         code will approximate the jacobian by forward-differencing.
c
c       n is a positive integer input variable set to the number of
c         functions and variables.
c
c       x is an array of length n.  on input x must contain an initial
c         estimate of the solution vector.  on output x contains the
c         final estimate of the solution vector.
c
c       fvec is an output array of length n which contains the functions
c         evaluated at the output x.
c
c       fjac is an output n by n array which contains the orthogonal
c         matrix q produced by the qr factorization of the final approx-
c         imate jacobian.
c
c       ldfjac is a positive integer input variable not less than n
c         which specifies the leading dimension of the array fjac.
c
c       xtol is a nonnegative input variable.  termination occurs when
c         the relative error between two consecutive iterates is at most
c         xtol.  therefore, xtol measures the relative error desired in
c         the approximate solution.  section 4 contains more details
c         about xtol.
c
c       maxfev is a positive integer input variable.  termination occurs
c         when the number of calls to fcn is at least maxfev by the end
c         of an iteration.
c
c       ml is a nonnegative integer input variable which specifies the
c         number of subdiagonals within the band of the jacobian matrix.
c         if the jacobian is not banded or iopt=1, set ml to at
c         least n - 1.
c
c       mu is a nonnegative integer input variable which specifies the
c         number of superdiagonals within the band of the jacobian
c         matrix.  if the jacobian is not banded or iopt=1, set mu to at
c         least n - 1.
c
c       epsfcn is an input variable used in determining a suitable step
c         for the forward-difference approximation.  this approximation
c         assumes that the relative errors in the functions are of the
c         order of epsfcn.  if epsfcn is less than the machine preci-
c         sion, it is assumed that the relative errors in the functions
c         are of the order of the machine precision.  if iopt=1, then
c         epsfcn can be ignored (treat it as a dummy argument).
c
c       diag is an array of length n.  if mode = 1 (see below), diag is
c         internally set.  if mode = 2, diag must contain positive
c         entries that serve as implicit (multiplicative) scale factors
c         for the variables.
c
c       mode is an integer input variable.  if mode = 1, the variables
c         will be scaled internally.  if mode = 2, the scaling is speci-
c         fied by the input diag.  other values of mode are equivalent
c         to mode = 1.
c
c       factor is a positive input variable used in determining the ini-
c         tial step bound.  this bound is set to the product of factor
c         and the euclidean norm of diag*x if nonzero, or else to factor
c         itself.  in most cases factor should lie in the interval
c         (.1,100.).  100. is a generally recommended value.
c
c       nprint is an integer input variable that enables controlled
c         printing of iterates if it is positive.  in this case, fcn is
c         called with iflag = 0 at the beginning of the first iteration
c         and every nprint iterations thereafter and immediately prior
c         to return, with x and fvec available for printing. appropriate
c         print statements must be added to fcn(see example).  if nprint
c         is not positive, no special calls of fcn with iflag = 0 are
c         made.
c
c       info is an integer output variable.  if the user has terminated
c         execution, info is set to the (negative) value of iflag.  see
c         description of fcn and jac. otherwise, info is set as follows.
c
c         info = 0  improper input parameters.
c
c         info = 1  relative error between two consecutive iterates is
c                   at most xtol.
c
c         info = 2  number of calls to fcn has reached or exceeded
c                   maxfev.
c
c         info = 3  xtol is too small.  no further improvement in the
c                   approximate solution x is possible.
c
c         info = 4  iteration is not making good progress, as measured
c                   by the improvement from the last five jacobian eval-
c                   uations.
c
c         info = 5  iteration is not making good progress, as measured
c                   by the improvement from the last ten iterations.
c
c         sections 4 and 5 contain more details about info.
c
c       nfev is an integer output variable set to the number of calls to
c         fcn.
c
c       njev is an integer output variable set to the number of calls to
c         jac. (if iopt=2, then njev is set to zero.)
c
c       r is an output array of length lr which contains the upper
c         triangular matrix produced by the qr factorization of the
c         final approximate jacobian, stored rowwise.
c
c       lr is a positive integer input variable not less than
c         (n*(n+1))/2.
c
c       qtf is an output array of length n which contains the vector
c         (q transpose)*fvec.
c
c       wa1, wa2, wa3, and wa4 are work arrays of length n.
c
c
c 4. successful completion.
c
c       the accuracy of snsq is controlled by the convergence parameter
c       xtol.  this parameter is used in a test which makes a comparison
c       between the approximation x and a solution xsol.  snsq termi-
c       nates when the test is satisfied.  if the convergence parameter
c       is less than the machine precision (as defined by the function
c       r1mach(4)), then snsq only attempts to satisfy the test
c       defined by the machine precision.  further progress is not
c       usually possible.
c
c       the test assumes that the functions are reasonably well behaved,
c       and, if the jacobian is supplied by the user, that the functions
c       and the jacobian are coded consistently.  if these conditions
c       are not satisfied, then snsq may incorrectly indicate conver-
c       gence.  the coding of the jacobian can be checked by the
c       subroutine chkder. if the jacobian is coded correctly or iopt=2,
c       then the validity of the answer can be checked, for example, by
c       rerunning snsq with a tighter tolerance.
c
c       convergence test.  if enorm(z) denotes the euclidean norm of a
c         vector z and d is the diagonal matrix whose entries are
c         defined by the array diag, then this test attempts to guaran-
c         tee that
c
c               enorm(d*(x-xsol)) .le. xtol*enorm(d*xsol).
c
c         if this condition is satisfied with xtol = 10**(-k), then the
c         larger components of d*x have k significant decimal digits and
c         info is set to 1.  there is a danger that the smaller compo-
c         nents of d*x may have large relative errors, but the fast rate
c         of convergence of snsq usually avoids this possibility.
c         unless high precision solutions are required, the recommended
c         value for xtol is the square root of the machine precision.
c
c
c 5. unsuccessful completion.
c
c       unsuccessful termination of snsq can be due to improper input
c       parameters, arithmetic interrupts, an excessive number of func-
c       tion evaluations, or lack of good progress.
c
c       improper input parameters.  info is set to 0 if iopt .lt .1,
c         or iopt .gt. 2, or n .le. 0, or ldfjac .lt. n, or
c         xtol .lt. 0.e0, or maxfev .le. 0, or ml .lt. 0, or mu .lt. 0,
c         or factor .le. 0.e0, or lr .lt. (n*(n+1))/2.
c
c       arithmetic interrupts.  if these interrupts occur in the fcn
c         subroutine during an early stage of the computation, they may
c         be caused by an unacceptable choice of x by snsq.  in this
c         case, it may be possible to remedy the situation by rerunning
c         snsq with a smaller value of factor.
c
c       excessive number of function evaluations.  a reasonable value
c         for maxfev is 100*(n+1) for iopt=1 and 200*(n+1) for iopt=2.
c         if the number of calls to fcn reaches maxfev, then this
c         indicates that the routine is converging very slowly as
c         measured by the progress of fvec, and info is set to 2. this
c         situation should be unusual because, as indicated below, lack
c         of good progress is usually diagnosed earlier by snsq,
c         causing termination with info = 4 or info = 5.
c
c       lack of good progress.  snsq searches for a zero of the system
c         by minimizing the sum of the squares of the functions.  in so
c         doing, it can become trapped in a region where the minimum
c         does not correspond to a zero of the system and, in this situ-
c         ation, the iteration eventually fails to make good progress.
c         in particular, this will happen if the system does not have a
c         zero.  if the system has a zero, rerunning snsq from a dif-
c         ferent starting point may be helpful.
c
c
c 6. characteristics of the algorithm.
c
c       snsq is a modification of the powell hybrid method.  two of its
c       main characteristics involve the choice of the correction as a
c       convex combination of the newton and scaled gradient directions,
c       and the updating of the jacobian by the rank-1 method of broy-
c       den.  the choice of the correction guarantees (under reasonable
c       conditions) global convergence for starting points far from the
c       solution and a fast rate of convergence.  the jacobian is
c       calculated at the starting point by either the user-supplied
c       subroutine or a forward-difference approximation, but it is not
c       recalculated until the rank-1 method fails to produce satis-
c       factory progress.
c
c       timing.  the time required by snsq to solve a given problem
c         depends on n, the behavior of the functions, the accuracy
c         requested, and the starting point.  the number of arithmetic
c         operations needed by snsq is about 11.5*(n**2) to process
c         each evaluation of the functions (call to fcn) and 1.3*(n**3)
c         to process each evaluation of the jacobian (call to jac,
c         if iopt = 1).  unless fcn and jac can be evaluated quickly,
c         the timing of snsq will be strongly influenced by the time
c         spent in fcn and jac.
c
c       storage.  snsq requires (3*n**2 + 17*n)/2 single precision
c         storage locations, in addition to the storage required by the
c         program.  there are no internally declared storage arrays.
c
c
c 7. example.
c
c       the problem is to determine the values of x(1), x(2), ..., x(9),
c       which solve the system of tridiagonal equations
c
c       (3-2*x(1))*x(1)           -2*x(2)                   = -1
c               -x(i-1) + (3-2*x(i))*x(i)         -2*x(i+1) = -1, i=2-8
c                                   -x(8) + (3-2*x(9))*x(9) = -1
c c     **********
c
c       program test(input,output,tape6=output)
c c
c c     driver for snsq example.
c c
c       integer j,iopt,n,maxfev,ml,mu,mode,nprint,info,nfev,ldfjac,lr,
c      *        nwrite
c       real xtol,epsfcn,factor,fnorm
c       real x(9),fvec(9),diag(9),fjac(9,9),r(45),qtf(9),
c      *     wa1(9),wa2(9),wa3(9),wa4(9)
c       real enorm,r1mach
c       external fcn
c       data nwrite /6/
c c
c       iopt = 2
c       n = 9
c c
c c     the following starting values provide a rough solution.
c c
c       do 10 j = 1, 9
c          x(j) = -1.e0
c    10    continue
c c
c       ldfjac = 9
c       lr = 45
c c
c c     set xtol to the square root of the machine precision.
c c     unless high precision solutions are required,
c c     this is the recommended setting.
c c
c       xtol = sqrt(r1mach(4))
c c
c       maxfev = 2000
c       ml = 1
c       mu = 1
c       epsfcn = 0.e0
c       mode = 2
c       do 20 j = 1, 9
c          diag(j) = 1.e0
c    20    continue
c       factor = 1.e2
c       nprint = 0
c c
c       call snsq(fcn,jac,iopt,n,x,fvec,fjac,ldfjac,xtol,maxfev,ml,mu,
c      *           epsfcn,diag,mode,factor,nprint,info,nfev,njev,
c      *           r,lr,qtf,wa1,wa2,wa3,wa4)
c       fnorm = enorm(n,fvec)
c       write (nwrite,1000) fnorm,nfev,info,(x(j),j=1,n)
c       stop
c  1000 format (5x,31h final l2 norm of the residuals,e15.7 //
c      *        5x,31h number of function evaluations,i10 //
c      *        5x,15h exit parameter,16x,i10 //
c      *        5x,27h final approximate solution // (5x,3e15.7))
c       end
c       subroutine fcn(n,x,fvec,iflag)
c       integer n,iflag
c       real x(n),fvec(n)
c       integer k
c       real one,temp,temp1,temp2,three,two,zero
c       data zero,one,two,three /0.e0,1.e0,2.e0,3.e0/
c c
c       if (iflag .ne. 0) go to 5
c c
c c     insert print statements here when nprint is positive.
c c
c       return
c     5 continue
c       do 10 k = 1, n
c          temp = (three - two*x(k))*x(k)
c          temp1 = zero
c          if (k .ne. 1) temp1 = x(k-1)
c          temp2 = zero
c          if (k .ne. n) temp2 = x(k+1)
c          fvec(k) = temp - temp1 - two*temp2 + one
c    10    continue
c       return
c       end
c
c       results obtained with different compilers or machines
c       may be slightly different.
c
c       final l2 norm of the residuals  0.1192636e-07
c
c       number of function evaluations        14
c
c       exit parameter                         1
c
c       final approximate solution
c
c       -0.5706545e+00 -0.6816283e+00 -0.7017325e+00
c       -0.7042129e+00 -0.7013690e+00 -0.6918656e+00
c       -0.6657920e+00 -0.5960342e+00 -0.4164121e+00
c
c
c***references
c       powell, m. j. d.
c       a hybrid method for nonlinear equations.
c       numerical methods for nonlinear algebraic equations,
c       p. rabinowitz, editor.  gordon and breach, 1970.
c***routines called  dogleg,enorm,fdjac1,qform,qrfac,r1mach,
c                    r1mpyq,r1updt,xerror
c***end prologue  snsq
      implicit real*8 (a-h,o-z)
      integer iopt,n,maxfev,ml,mu,mode,nprint,info,nfev,ldfjac,lr,njev
      real*8 xtol
      real*8 epsfcn,factor
      real*8 x(n),fvec(n),diag(n),fjac(ldfjac,n),r(lr),qtf(n),wa1(n),
     *     wa2(n),wa3(n),wa4(n)
      external fcn
      integer i,iflag,iter,j,jm1,l,ncfail,ncsuc,nslow1,nslow2
      integer iwa(1)
      logical jeval,sing
      real*8 actred,delta,epsmch,fnorm,fnorm1,one,pnorm,prered,p1,p5,
     *     p001,p0001,ratio,sum,temp,xnorm,zero
      real*8 r1mach,enorm
      data one,p1,p5,p001,p0001,zero
     *     /1.0d0,1.0d-1,5.0d-1,1.0d-3,1.0d-4,0.0d0/
c
c***first executable statement  snsq
      epsmch = r1mach(4)
c
      info = 0
      iflag = 0
      nfev = 0
      njev = 0
c
c     check the input parameters for errors.
c
      if (iopt .lt. 1 .or. iopt .gt. 2 .or.
     *    n .le. 0 .or. xtol .lt. zero .or. maxfev .le. 0
     *    .or. ml .lt. 0 .or. mu .lt. 0 .or. factor .le. zero
     *    .or. ldfjac .lt. n .or. lr .lt. (n*(n + 1))/2) go to 300
      if (mode .ne. 2) go to 20
      do 10 j = 1, n
	 if (diag(j) .le. zero) go to 300
   10    continue
   20 continue
c
c     evaluate the function at the starting point
c     and calculate its norm.
c
      iflag = 1
      call fcn(n,x,fvec,iflag)
      nfev = 1
      if (iflag .lt. 0) go to 300
      fnorm = enorm(n,fvec)
c
c     initialize iteration counter and monitors.
c
      iter = 1
      ncsuc = 0
      ncfail = 0
      nslow1 = 0
      nslow2 = 0
c
c     beginning of the outer loop.
c
   30 continue
	 jeval = .true.
c
c        calculate the jacobian matrix.
c
	 if (iopt .eq. 2) go to 31
c
c        user supplies jacobian
c
	    call jac(n,x,fvec,fjac,ldfjac,iflag)
	    njev = njev+1
	    go to 32
c
c        code approximates the jacobian
c
   31       iflag = 2
	    call fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,wa1,
     *               wa2)
	    nfev = nfev + min0(ml+mu+1,n)
c
   32    if (iflag .lt. 0) go to 300
c
c        compute the qr factorization of the jacobian.
c
	 call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
c
c        on the first iteration and if mode is 1, scale according
c        to the norms of the columns of the initial jacobian.
c
	 if (iter .ne. 1) go to 70
	 if (mode .eq. 2) go to 50
	 do 40 j = 1, n
	    diag(j) = wa2(j)
	    if (wa2(j) .eq. zero) diag(j) = one
   40       continue
   50    continue
c
c        on the first iteration, calculate the norm of the scaled x
c        and initialize the step bound delta.
c
	 do 60 j = 1, n
	    wa3(j) = diag(j)*x(j)
   60       continue
	 xnorm = enorm(n,wa3)
	 delta = factor*xnorm
	 if (delta .eq. zero) delta = factor
   70    continue
c
c        form (q transpose)*fvec and store in qtf.
c
	 do 80 i = 1, n
	    qtf(i) = fvec(i)
   80       continue
	 do 120 j = 1, n
	    if (fjac(j,j) .eq. zero) go to 110
	    sum = zero
	    do 90 i = j, n
	       sum = sum + fjac(i,j)*qtf(i)
   90          continue
	    temp = -sum/fjac(j,j)
	    do 100 i = j, n
	       qtf(i) = qtf(i) + fjac(i,j)*temp
  100          continue
  110       continue
  120       continue
c
c        copy the triangular factor of the qr factorization into r.
c
	 sing = .false.
	 do 150 j = 1, n
	    l = j
	    jm1 = j - 1
	    if (jm1 .lt. 1) go to 140
	    do 130 i = 1, jm1
	       r(l) = fjac(i,j)
	       l = l + n - i
  130          continue
  140       continue
	    r(l) = wa1(j)
	    if (wa1(j) .eq. zero) sing = .true.
  150       continue
c
c        accumulate the orthogonal factor in fjac.
c
	 call qform(n,n,fjac,ldfjac,wa1)
c
c        rescale if necessary.
c
	 if (mode .eq. 2) go to 170
	 do 160 j = 1, n
	    diag(j) = dmax1(diag(j),wa2(j))
  160       continue
  170    continue
c
c        beginning of the inner loop.
c
  180    continue
c
c           if requested, call fcn to enable printing of iterates.
c
	    if (nprint .le. 0) go to 190
	    iflag = 0
	    if (mod(iter-1,nprint) .eq. 0) call fcn(n,x,fvec,iflag)
	    if (iflag .lt. 0) go to 300
  190       continue
c
c           determine the direction p.
c
	    call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
c
c           store the direction p and x + p. calculate the norm of p.
c
	    do 200 j = 1, n
	       wa1(j) = -wa1(j)
	       wa2(j) = x(j) + wa1(j)
	       wa3(j) = diag(j)*wa1(j)
  200          continue
	    pnorm = enorm(n,wa3)
c
c           on the first iteration, adjust the initial step bound.
c
	    if (iter .eq. 1) delta = dmin1(delta,pnorm)
c
c           evaluate the function at x + p and calculate its norm.
c
	    iflag = 1
	    call fcn(n,wa2,wa4,iflag)
	    nfev = nfev + 1
	    if (iflag .lt. 0) go to 300
	    fnorm1 = enorm(n,wa4)
c
c           compute the scaled actual reduction.
c
	    actred = -one
	    if (fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
c
c           compute the scaled predicted reduction.
c
	    l = 1
	    do 220 i = 1, n
	       sum = zero
	       do 210 j = i, n
		  sum = sum + r(l)*wa1(j)
		  l = l + 1
  210             continue
	       wa3(i) = qtf(i) + sum
  220          continue
	    temp = enorm(n,wa3)
	    prered = zero
	    if (temp .lt. fnorm) prered = one - (temp/fnorm)**2
c
c           compute the ratio of the actual to the predicted
c           reduction.
c
	    ratio = zero
	    if (prered .gt. zero) ratio = actred/prered
c
c           update the step bound.
c
	    if (ratio .ge. p1) go to 230
	       ncsuc = 0
	       ncfail = ncfail + 1
	       delta = p5*delta
	       go to 240
  230       continue
	       ncfail = 0
	       ncsuc = ncsuc + 1
	       if (ratio .ge. p5 .or. ncsuc .gt. 1)
     *            delta = dmax1(delta,pnorm/p5)
	       if (abs(ratio-one) .le. p1) delta = pnorm/p5
  240       continue
c
c           test for successful iteration.
c
	    if (ratio .lt. p0001) go to 260
c
c           successful iteration. update x, fvec, and their norms.
c
	    do 250 j = 1, n
	       x(j) = wa2(j)
	       wa2(j) = diag(j)*x(j)
	       fvec(j) = wa4(j)
  250          continue
	    xnorm = enorm(n,wa2)
	    fnorm = fnorm1
	    iter = iter + 1
  260       continue
c
c           determine the progress of the iteration.
c
	    nslow1 = nslow1 + 1
	    if (actred .ge. p001) nslow1 = 0
	    if (jeval) nslow2 = nslow2 + 1
	    if (actred .ge. p1) nslow2 = 0
c
c           test for convergence.
c
	    if (delta .le. xtol*xnorm .or. fnorm .eq. zero) info = 1
	    if (info .ne. 0) go to 300
c
c           tests for termination and stringent tolerances.
c
	    if (nfev .ge. maxfev) info = 2
	    if (p1*dmax1(p1*delta,pnorm) .le. epsmch*xnorm) info = 3
	    if (nslow2 .eq. 5) info = 4
	    if (nslow1 .eq. 10) info = 5
	    if (info .ne. 0) go to 300
c
c           criterion for recalculating jacobian
c
	    if (ncfail .eq. 2) go to 290
c
c           calculate the rank one modification to the jacobian
c           and update qtf if necessary.
c
	    do 280 j = 1, n
	       sum = zero
	       do 270 i = 1, n
		  sum = sum + fjac(i,j)*wa4(i)
  270             continue
	       wa2(j) = (sum - wa3(j))/pnorm
	       wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
	       if (ratio .ge. p0001) qtf(j) = sum
  280          continue
c
c           compute the qr factorization of the updated jacobian.
c
	    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)
c
c           end of the inner loop.
c
	    jeval = .false.
	    go to 180
  290    continue
c
c        end of the outer loop.
c
	 go to 30
  300 continue
c
c     termination, either normal or user imposed.
c
      if (iflag .lt. 0) info = iflag
      iflag = 0
      if (nprint .gt. 0) call fcn(n,x,fvec,iflag)
      if (info .lt. 0) call xerror(63hsnsq   -- execution terminated bec
     1ause user set iflag negative.,63,1,1)
      if (info .eq. 0) call xerror(34hsnsq   -- invalid input parameter.
     1 ,34,2,1)
      if (info .eq. 2) call xerror(40hsnsq   -- too many function evalua
     1tions.,40,9,1)
      if (info .eq. 3) call xerror(58hsnsq   -- xtol too small. no furth
     1er improvement possible.,58,3,1)
      if (info .gt. 4) call xerror(45hsnsq   -- iteration not making goo
     1d progress.,45,1,1)
      return
c
c     last card of subroutine snsq.
c
      end
C=================================================================
      subroutine xerabt(messg,nmessg)
c***begin prologue  xerabt
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  q
c***keywords  xerror,error,abort
c***date written  august 1979
c***author  jones r.e. (sla)
c***purpose
c aborts program execution and prints error message
c***description
c     abstract
c        ***note*** machine dependent routine
c        xerabt aborts the execution of the program.
c        the error message causing the abort is given in the calling
c        sequence in case one needs it for printing on a dayfile,
c        for example.
c
c     description of parameters
c        messg and nmessg are as in xerror, except that nmessg may
c        be zero, in which case no message is being supplied.
c
c     written by ron jones, with slatec common math library subcommittee
c     latest revision ---  7 june 1978
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  (none)
c***end prologue  xerabt
c
      dimension messg(nmessg)
c***first executable statement  xerabt
      stop
c     return
      end
      subroutine xerctl(messg1,nmessg,nerr,level,kontrl)
c***begin prologue  xerctl
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  q
c***keywords  xerror,error,control
c***date written  august 1979
c***author  jones r.e. (sla)
c***purpose
c allows user control over handling of individual errors
c***description
c     abstract
c        allows user control over handling of individual errors.
c        just after each message is recorded, but before it is
c        processed any further (i.e., before it is printed or
c        a decision to abort is made) a call is made to xerctl.
c        if the user has provided his own version of xerctl, he
c        can then override the value of kontrol used in processing
c        this message by redefining its value.
c        kontrl may be set to any value from -2 to 2.
c        the meanings for kontrl are the same as in xsetf, except
c        that the value of kontrl changes only for this message.
c        if kontrl is set to a value outside the range from -2 to 2,
c        it will be moved back into that range.
c
c     description of parameters
c
c      --input--
c        messg1 - the first word (only) of the error message.
c        nmessg - same as in the call to xerror or xerrwv.
c        nerr   - same as in the call to xerror or xerrwv.
c        level  - same as in the call to xerror or xerrwv.
c        kontrl - the current value of the control flag as set
c                 by a call to xsetf.
c
c      --output--
c        kontrl - the new value of kontrl.  if kontrl is not
c                 defined, it will remain at its original value.
c                 this changed value of control affects only
c                 the current occurrence of the current message.
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  (none)
c***end prologue  xerctl
c
c***first executable statement  xerctl
      return
      end
      subroutine xerprt(messg,nmessg)
c***begin prologue  xerprt
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  q
c***keywords  xerror,error,print
c***date written  august 1979
c***author  jones r.e. (sla)
c***purpose
c prints error messages
c***description
c     abstract
c        print the hollerith message in messg, of length nmessg,
c        on each file indicated by xgetua.
c
c     prepare format
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  xgetua,s88fmt,i1mach
c***end prologue  xerprt
c
      implicit real*8 (a-h,o-z)
      integer f(10),g(14),lun(5)
      dimension messg(nmessg)
      data f(1),f(2),f(3),f(4),f(5),f(6),f(7),f(8),f(9),f(10)
     1   / 1h( ,1h1 ,1hx ,1h, ,1h  ,1h  ,1ha ,1h  ,1h  ,1h) /
      data g(1),g(2),g(3),g(4),g(5),g(6),g(7),g(8),g(9),g(10)
     1   / 1h( ,1h1 ,1hx ,1h  ,1h  ,1h  ,1h  ,1h  ,1h  ,1h  /
      data g(11),g(12),g(13),g(14)
     1   / 1h   ,1h   ,1h   ,1h)  /
      data la/1ha/,lcom/1h,/,lblank/1h /
c    prepare format for whole lines
c***first executable statement  xerprt
      nchar = i1mach(6)
      nfield = 72/nchar
      call s88fmt(2,nfield,f(5))
      call s88fmt(2,nchar,f(8))
c     prepare format for last, partial line, if needed
      ncharl = nfield*nchar
      nlines = nmessg/ncharl
      nword  = nlines*nfield
      nchrem = nmessg - nlines*ncharl
      if (nchrem.le.0) go to 40
	 do 10 i=4,13
10          g(i) = lblank
	 nfield = nchrem/nchar
	 if (nfield.le.0) go to 20
c        prepare whole word fields
	    g(4) = lcom
	    call s88fmt(2,nfield,g(5))
	    g(7) = la
	    call s88fmt(2,nchar,g(8))
20       continue
	 nchlst = mod(nchrem,nchar)
	 if (nchlst.le.0) go to 30
c        prepare partial word field
	    g(10) = lcom
	    g(11) = la
	    call s88fmt(2,nchlst,g(12))
30       continue
40    continue
c     print the message
      nword1 = nword+1
      nword2 = (nmessg+nchar-1)/nchar
      call xgetua(lun,nunit)
      do 50 kunit = 1,nunit
	 iunit = lun(kunit)
	 if (iunit.eq.0) iunit = i1mach(4)
	 if (nword.gt.0) write (iunit,f) (messg(i),i=1,nword)
	 if (nchrem.gt.0) write (iunit,g) (messg(i),i=nword1,nword2)
50    continue
      return
      end
      subroutine xerror(messg,nmessg,nerr,level)
c***begin prologue  xerror
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  q
c***keywords  xerror,error
c***date written  august 1979
c***author  jones r.e. (sla)
c***purpose
c processes an error (diagnostic) message
c***description
c     abstract
c        xerror processes a diagnostic message, in a manner
c        determined by the value of level and the current value
c        of the library error control flag, kontrl.
c        (see subroutine xsetf for details.)
c
c     description of parameters
c      --input--
c        messg - the hollerith message to be processed, containing
c                no more than 72 characters.
c        nmessg- the actual number of characters in messg.
c        nerr  - the error number associated with this message.
c                nerr must not be zero.
c        level - error category.
c                =2 means this is an unconditionally fatal error.
c                =1 means this is a recoverable error.  (i.e., it is
c                   non-fatal if xsetf has been appropriately called.)
c                =0 means this is a warning message only.
c                =-1 means this is a warning message which is to be
c                   printed at most once, regardless of how many
c                   times this call is executed.
c
c     examples
c        call xerror(23hsmooth -- num was zero.,23,1,2)
c        call xerror(43hinteg  -- less than full accuracy achieved.,
c                    43,2,1)
c        call xerror(65hrooter -- actual zero of f found before interval
c    1 fully collapsed.,65,3,0)
c        call xerror(39hexp    -- underflows being set to zero.,39,1,-1)
c
c     written by ron jones, with slatec common math library subcommittee
c     latest revision ---  7 feb 1979
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  xerrwv
c***end prologue  xerror
c
      implicit real*8 (a-h,o-z)
      dimension messg(nmessg)
c***first executable statement  xerror
      call xerrwv(messg,nmessg,nerr,level,0,0,0,0,0.d00,0.d00)
      return
      end
      subroutine xerrwv(messg,nmessg,nerr,level,ni,i1,i2,nr,r1,r2)
c***begin prologue  xerrwv
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  q
c***keywords  xerror,error,print with values
c***date written  march 19, 1980
c***author  jones r.e. (sla)
c***purpose
c processes error message allowing 2 integer and two real values
c  to be included in the message.
c***description
c     abstract
c        xerrwv processes a diagnostic message, in a manner
c        determined by the value of level and the current value
c        of the library error control flag, kontrl.
c        (see subroutine xsetf for details.)
c        in addition, up to two integer values and two real
c        values may be printed along with the message.
c
c     description of parameters
c      --input--
c        messg - the hollerith message to be processed.
c        nmessg- the actual number of characters in messg.
c        nerr  - the error number associated with this message.
c                nerr must not be zero.
c        level - error category.
c                =2 means this is an unconditionally fatal error.
c                =1 means this is a recoverable error.  (i.e., it is
c                   non-fatal if xsetf has been appropriately called.)
c                =0 means this is a warning message only.
c                =-1 means this is a warning message which is to be
c                   printed at most once, regardless of how many
c                   times this call is executed.
c        ni    - number of integer values to be printed. (o to 2)
c        i1    - first integer value.
c        i2    - second integer value.
c        nr    - number of real values to be printed. (0 to 2)
c        r1    - first real value.
c        r2    - second real value.
c
c     examples
c        call xerror(29hsmooth -- num (=i1) was zero.,29,1,2,
c    1   1,num,0,0,0.,0.)
c        call xerrwv(54hquadxy -- requested error (r1) less than minimum
c    1 (r2).,54,77,1,0,0,0,2,errreq,errmin)
c
c     written by ron jones, with slatec common math library subcommittee
c     latest revision ---  19 mar 1980
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  fdump,i1mach,j4save,xerabt,xerctl,xerprt,
c                      xersav
c***end prologue  xerrwv
c
      implicit real*8 (a-h,o-z)
      dimension messg(nmessg),lun(5)
c***first executable statement  xerrwv
      lkntrl = j4save(2,0,.false.)
      maxmes = j4save(4,0,.false.)
c     check for valid input
      if ((nmessg.gt.0).and.(nerr.ne.0).and.
     1    (level.ge.(-1)).and.(level.le.2)) go to 10
	 if (lkntrl.gt.0) call xerprt(17hfatal error in...,17)
	 call xerprt(23hxerror -- invalid input,23)
	 if (lkntrl.gt.0) call fdump
	 if (lkntrl.gt.0) call xerprt(29hjob abort due to fatal error.,
     1   29)
	 if (lkntrl.gt.0) call xersav(1h ,0,0,0,kdummy)
	 call xerabt(23hxerror -- invalid input,23)
	 return
   10 continue
c     record message
      junk = j4save(1,nerr,.true.)
      call xersav(messg,nmessg,nerr,level,kount)
c     let user override
      lfirst = messg(1)
      lmessg = nmessg
      lerr = nerr
      llevel = level
      call xerctl(lfirst,lmessg,lerr,llevel,lkntrl)
c     reset to original values
      lmessg = nmessg
      lerr = nerr
      llevel = level
      lkntrl = max0(-2,min0(2,lkntrl))
      mkntrl = iabs(lkntrl)
c     decide whether to print message
      if ((llevel.lt.2).and.(lkntrl.eq.0)) go to 100
      if (((llevel.eq.(-1)).and.(kount.gt.min0(1,maxmes)))
     1.or.((llevel.eq.0)   .and.(kount.gt.maxmes))
     2.or.((llevel.eq.1)   .and.(kount.gt.maxmes).and.(mkntrl.eq.1))
     3.or.((llevel.eq.2)   .and.(kount.gt.max0(1,maxmes)))) go to 100
	 if (lkntrl.le.0) go to 20
	    call xerprt(1h ,1)
c           introduction
	    if (llevel.eq.(-1)) call xerprt
     1(57hwarning message...this message will only be printed once.,57)
	    if (llevel.eq.0) call xerprt(13hwarning in...,13)
	    if (llevel.eq.1) call xerprt
     1      (23hrecoverable error in...,23)
	    if (llevel.eq.2) call xerprt(17hfatal error in...,17)
   20    continue
c        message
	 call xerprt(messg,lmessg)
	 call xgetua(lun,nunit)
	 do 50 kunit=1,nunit
	    iunit = lun(kunit)
	    if (iunit.eq.0) iunit = i1mach(4)
	    if (ni.ge.1) write (iunit,22) i1
	    if (ni.ge.2) write (iunit,23) i2
	    if (nr.ge.1) write (iunit,24) r1
	    if (nr.ge.2) write (iunit,25) r2
   22       format (11x,21hin above message, i1=,i10)
   23       format (11x,21hin above message, i2=,i10)
   24       format (11x,21hin above message, r1=,e20.10)
   25       format (11x,21hin above message, r2=,e20.10)
	    if (lkntrl.le.0) go to 40
c              error number
	       write (iunit,30) lerr
   30          format (15h error number =,i10)
   40       continue
   50    continue
c        trace-back
	 call fdump
  100 continue
      ifatal = 0
      if ((llevel.eq.2).or.((llevel.eq.1).and.(mkntrl.eq.2)))
     1ifatal = 1
c     quit here if message is not fatal
      if (ifatal.le.0) return
      if ((lkntrl.le.0).or.(kount.gt.max0(1,maxmes))) go to 120
c        print reason for abort
	 if (llevel.eq.1) call xerprt
     1   (35hjob abort due to unrecovered error.,35)
	 if (llevel.eq.2) call xerprt
     1   (29hjob abort due to fatal error.,29)
c        print error summary
	 call xersav(1h ,-1,0,0,kdummy)
  120 continue
c     abort
      if ((llevel.eq.2).and.(kount.gt.max0(1,maxmes))) lmessg = 0
      call xerabt(messg,lmessg)
      return
      end
      subroutine xersav(messg,nmessg,nerr,level,icount)
c***begin prologue  xersav
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  q
c***keywords  xerror,error,record
c***date written  march 19, 1980
c***author  jones r.e. (sla)
c***purpose
c records that an error occurred.
c***description
c     abstract
c        record that this error occurred.
c
c     description of parameters
c     --input--
c       messg, nmessg, nerr, level are as in xerror,
c       except that when nmessg=0 the tables will be
c       dumped and cleared, and when nmessg is less than zero the
c       tables will be dumped and not cleared.
c     --output--
c       icount will be the number of times this message has
c       been seen, or zero if the table has overflowed and
c       does not contain this message specifically.
c       when nmessg=0, icount will not be altered.
c
c     written by ron jones, with slatec common math library subcommittee
c     latest revision ---  19 mar 1980
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  xgetua,s88fmt,i1mach
c***end prologue  xersav
c
      implicit real*8 (a-h,o-z)
      integer f(17),lun(5)
      dimension messg(1)
      dimension mestab(10),nertab(10),levtab(10),kount(10)
      data mestab(1),mestab(2),mestab(3),mestab(4),mestab(5),
     1     mestab(6),mestab(7),mestab(8),mestab(9),mestab(10)
     2     /0,0,0,0,0,0,0,0,0,0/
      data nertab(1),nertab(2),nertab(3),nertab(4),nertab(5),
     1     nertab(6),nertab(7),nertab(8),nertab(9),nertab(10)
     2     /0,0,0,0,0,0,0,0,0,0/
      data levtab(1),levtab(2),levtab(3),levtab(4),levtab(5),
     1     levtab(6),levtab(7),levtab(8),levtab(9),levtab(10)
     2     /0,0,0,0,0,0,0,0,0,0/
      data kount(1),kount(2),kount(3),kount(4),kount(5),
     1     kount(6),kount(7),kount(8),kount(9),kount(10)
     2     /0,0,0,0,0,0,0,0,0,0/
      data kountx/0/
      data f(1),f(2),f(3),f(4),f(5),f(6),f(7),f(8),f(9),f(10),
     1     f(11),f(12),f(13),f(14),f(15),f(16),f(17)
     2     /1h( ,1h1 ,1hx ,1h, ,1ha ,1h  ,1h  ,1h, ,1hi ,1h  ,
     3      1h  ,1h, ,1h2 ,1hi ,1h1 ,1h0 ,1h) /
c***first executable statement  xersav
      if (nmessg.gt.0) go to 80
c     dump the table
	 if (kount(1).eq.0) return
c        prepare format
	 nchar = i1mach(6)
	 call s88fmt(2,nchar,f(6))
	 ncol = 20 - nchar
	 call s88fmt(2,ncol,f(10))
c        print to each unit
	 call xgetua(lun,nunit)
	 do 60 kunit=1,nunit
	    iunit = lun(kunit)
	    if (iunit.eq.0) iunit = i1mach(4)
c           print table header
	    write (iunit,10)
   10       format (32h0          error message summary/
     1              41h first word      nerr     level     count)
c           print body of table
	    do 20 i=1,10
	       if (kount(i).eq.0) go to 30
	       write (iunit,f) mestab(i),nertab(i),levtab(i),kount(i)
   20       continue
   30       continue
c           print number of other errors
	    if (kountx.ne.0) write (iunit,40) kountx
   40       format (41h0other errors not individually tabulated=,i10)
	    write (iunit,50)
   50       format (1x)
   60    continue
	 if (nmessg.lt.0) return
c        clear the error tables
	 do 70 i=1,10
   70       kount(i) = 0
	 kountx = 0
	 return
   80 continue
c     process a message...
c     search for this messg, or else an empty slot for this messg,
c     or else determine that the error table is full.
      do 90 i=1,10
	 ii = i
	 if (kount(i).eq.0) go to 110
	 if (messg(1).ne.mestab(i)) go to 90
	 if (nerr.ne.nertab(i)) go to 90
	 if (level.ne.levtab(i)) go to 90
	 go to 100
   90 continue
c     three possible cases...
c     table is full
	 kountx = kountx+1
	 icount = 1
	 return
c     message found in table
  100    kount(ii) = kount(ii) + 1
	 icount = kount(ii)
	 return
c     empty slot found for new message
  110    mestab(ii) = messg(1)
	 nertab(ii) = nerr
	 levtab(ii) = level
	 kount(ii)  = 1
	 icount = 1
	 return
      end
      subroutine xgetua(iunit,n)
c***begin prologue  xgetua
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  q
c***keywords  xerror,error,unit number
c***date written  august 1979
c***author  jones r.e. (sla)
c***purpose
c returns unit number(s) to which error messages are being sent
c***description
c     abstract
c        xgetua may be called to determine the unit number or numbers
c        to which error messages are being sent.
c        these unit numbers may have been set by a call to xsetun,
c        or a call to xsetua, or may be a default value.
c
c     description of parameters
c      --output--
c        iunit - an array of one to five unit numbers, depending
c                on the value of n.  a value of zero refers to the
c                default unit, as defined by the i1mach machine
c                constant routine.  only iunit(1),...,iunit(n) are
c                defined by xgetua.  the values of iunit(n+1),...,
c                iunit(5) are not defined (for n.lt.5) or altered
c                in any way by xgetua.
c        n     - the number of units to which copies of the
c                error messages are being sent.  n will be in the
c                range from 1 to 5.
c
c     written by ron jones, with slatec common math library subcommittee
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called   j4save
c***end prologue  xgetua
c
      implicit real*8 (a-h,o-z)
      dimension iunit(5)
c***first executable statement  xgetua
      n = j4save(5,0,.false.)
      do 30 i=1,n
	 index = i+4
	 if (i.eq.1) index = 3
	 iunit(i) = j4save(index,0,.false.)
   30 continue
      return
      end
C==============================================================