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