relibc_openlibm/slatec/snls1.f

*DECK SNLS1
      SUBROUTINE SNLS1 (FCN, 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  SNLS1
C***PURPOSE  Minimize the sum of the squares of M nonlinear functions
C            in N variables by a modification of the Levenberg-Marquardt
C            algorithm.
C***LIBRARY   SLATEC
C***CATEGORY  K1B1A1, K1B1A2
C***TYPE      SINGLE PRECISION (SNLS1-S, DNLS1-D)
C***KEYWORDS  LEVENBERG-MARQUARDT, NONLINEAR DATA FITTING,
C             NONLINEAR LEAST SQUARES
C***AUTHOR  Hiebert, K. L., (SNLA)
C***DESCRIPTION
C
C 1. Purpose.
C
C       The purpose of SNLS1 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 SNLS1(FCN,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       REAL FTOL,XTOL,GTOL,EPSFCN,FACTOR
C       REAL 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 SNLS1 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 SNLS1.
C
C       FCN is the name of the user-supplied subroutine which calculates
C         the functions.  If the user wants to supply the Jacobian
C         (IOPT=2 or 3), then FCN must be written to calculate the
C         Jacobian, as well as the functions.  See the explanation
C         of the IOPT argument below.
C         If the user wants the iterates printed (NPRINT positive), then
C         FCN must do the printing.  See the explanation of NPRINT
C         below.  FCN must be declared in an EXTERNAL statement in the
C         calling program and should be written as follows.
C
C
C         SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
C         INTEGER IFLAG,LDFJAC,M,N
C         REAL X(N),FVEC(M)
C         ----------
C         FJAC and LDFJAC may be ignored     , if IOPT=1.
C         REAL FJAC(LDFJAC,N)                , if IOPT=2.
C         REAL FJAC(N)                       , if IOPT=3.
C         ----------
C           If IFLAG=0, the values in X and FVEC are available
C           for printing.  See the explanation of NPRINT below.
C           IFLAG will never be zero unless NPRINT is positive.
C           The values of X and FVEC must not be changed.
C         RETURN
C         ----------
C           If IFLAG=1, calculate the functions at X and return
C           this vector in FVEC.
C         RETURN
C         ----------
C           If IFLAG=2, calculate the full Jacobian at X and return
C           this matrix in FJAC.  Note that IFLAG will never be 2 unless
C           IOPT=2.  FVEC contains the function values at X and must
C           not be altered.  FJAC(I,J) must be set to the derivative
C           of FVEC(I) with respect to X(J).
C         RETURN
C         ----------
C           If IFLAG=3, calculate the LDFJAC-th row of the Jacobian
C           and return this vector in FJAC.  Note that IFLAG will
C           never be 3 unless IOPT=3.  FVEC contains the function
C           values at X and must not be altered.  FJAC(J) must be
C           set to the derivative of FVEC(LDFJAC) with respect to X(J).
C         RETURN
C         ----------
C         END
C
C
C         The value of IFLAG should not be changed by FCN unless the
C         user wants to terminate execution of SNLS1.  In this case, set
C         IFLAG to a negative integer.
C
C
C       IOPT is an input variable which specifies how the Jacobian will
C         be calculated.  If IOPT=2 or 3, then the user must supply the
C         Jacobian, as well as the function values, through the
C         subroutine FCN.  If IOPT=2, the user supplies the full
C         Jacobian with one call to FCN.  If IOPT=3, the user supplies
C         one row of the Jacobian with each call.  (In this manner,
C         storage can be saved because the full Jacobian is not stored.)
C         If IOPT=1, the code will approximate the Jacobian by forward
C         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 non-negative 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 non-negative 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 non-negative 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 to evaluate the functions
C         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=2 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 to 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 for function evaluation
C                   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 for function evaluation.
C
C       NJEV is an integer output variable set to the number of
C         evaluations of the full Jacobian.  If IOPT=2, only one call to
C         FCN is required for each evaluation of the full Jacobian.
C         If IOPT=3, the M calls to FCN are required.
C         If IOPT=1, 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 SNLS1 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.  SNLS1 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 SNLS1 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 SNLS1 may incorrectly indicate conver-
C       gence.  If the Jacobian is coded correctly or IOPT=1,
C       then the validity of the answer can be checked, for example, by
C       rerunning SNLS1 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 SNLS1, 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 SNLS1 can be due to improper input
C       parameters, arithmetic interrupts, or an excessive number of
C       function evaluations.
C
C       Improper Input Parameters.  INFO is set to 0 if IOPT .LT. 1
C         or IOPT .GT. 3, or N .LE. 0, or M .LT. N, or for IOPT=1 or 2
C         LDFJAC .LT. M, or for IOPT=3 LDFJAC .LT. N, or FTOL .LT. 0.E0,
C         or XTOL .LT. 0.E0, or GTOL .LT. 0.E0, or MAXFEV .LE. 0, or
C         FACTOR .LE. 0.E0.
C
C       Arithmetic Interrupts.  If these interrupts occur in the FCN
C         subroutine during an early stage of the computation, they may
C         be caused by an unacceptable choice of X by SNLS1.  In this
C         case, it may be possible to remedy the situation by rerunning
C         SNLS1 with a smaller value of FACTOR.
C
C       Excessive Number of Function Evaluations.  A reasonable value
C         for MAXFEV is 100*(N+1) for IOPT=2 or 3 and 200*(N+1) for
C         IOPT=1.  If the number of calls to FCN reaches MAXFEV, then
C         this indicates that the routine is converging very slowly
C         as measured by the progress of FVEC, and INFO is set to 5.
C         In this case, it may be helpful to restart SNLS1 with MODE
C         set to 1.
C
C
C 6. Characteristics of the Algorithm.
C
C       SNLS1 is a modification of the Levenberg-Marquardt algorithm.
C       Two of its main characteristics involve the proper use of
C       implicitly scaled variables (if MODE = 1) and an optimal choice
C       for the correction.  The use of implicitly scaled variables
C       achieves scale invariance of SNLS1 and limits the size of the
C       correction in any direction where the functions are changing
C       rapidly.  The optimal choice of the correction guarantees (under
C       reasonable conditions) global convergence from starting points
C       far from the solution and a fast rate of convergence for
C       problems with small residuals.
C
C       Timing.  The time required by SNLS1 to solve a given problem
C         depends on M and N, the behavior of the functions, the accu-
C         racy requested, and the starting point.  The number of arith-
C         metic operations needed by SNLS1 is about N**3 to process each
C         evaluation of the functions (call to FCN) and to process each
C         evaluation of the Jacobian it takes M*N**2 for IOPT=2 (one
C         call to FCN), M*N**2 for IOPT=1 (N calls to FCN) and
C         1.5*M*N**2 for IOPT=3 (M calls to FCN).  Unless FCN
C         can be evaluated quickly, the timing of SNLS1 will be
C         strongly influenced by the time spent in FCN.
C
C       Storage.  SNLS1 requires (M*N + 2*M + 6*N) for IOPT=1 or 2 and
C         (N**2 + 2*M + 6*N) for IOPT=3 single precision storage
C         locations and N integer storage locations, in addition to
C         the storage required by the program.  There are no internally
C         declared storage arrays.
C
C *Long Description:
C
C 7. Example.
C
C       The problem is to determine the values of X(1), X(2), and X(3)
C       which provide the best fit (in the least squares sense) of
C
C             X(1) + U(I)/(V(I)*X(2) + W(I)*X(3)),  I = 1, 15
C
C       to the data
C
C             Y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,
C                  0.37,0.58,0.73,0.96,1.34,2.10,4.39),
C
C       where U(I) = I, V(I) = 16 - I, and W(I) = MIN(U(I),V(I)).  The
C       I-th component of FVEC is thus defined by
C
C             Y(I) - (X(1) + U(I)/(V(I)*X(2) + W(I)*X(3))).
C
C       **********
C
C       PROGRAM TEST
C C
C C     Driver for SNLS1 example.
C C
C       INTEGER J,IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,
C      *        NWRITE
C       INTEGER IPVT(3)
C       REAL FTOL,XTOL,GTOL,FACTOR,FNORM,EPSFCN
C       REAL X(3),FVEC(15),FJAC(15,3),DIAG(3),QTF(3),
C      *     WA1(3),WA2(3),WA3(3),WA4(15)
C       REAL ENORM,R1MACH
C       EXTERNAL FCN
C       DATA NWRITE /6/
C C
C       IOPT = 1
C       M = 15
C       N = 3
C C
C C     The following starting values provide a rough fit.
C C
C       X(1) = 1.E0
C       X(2) = 1.E0
C       X(3) = 1.E0
C C
C       LDFJAC = 15
C C
C C     Set FTOL and XTOL to the square root of the machine precision
C C     and GTOL to zero.  Unless high precision solutions are
C C     required, these are the recommended settings.
C C
C       FTOL = SQRT(R1MACH(4))
C       XTOL = SQRT(R1MACH(4))
C       GTOL = 0.E0
C C
C       MAXFEV = 400
C       EPSFCN = 0.0
C       MODE = 1
C       FACTOR = 1.E2
C       NPRINT = 0
C C
C       CALL SNLS1(FCN,IOPT,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,
C      *           GTOL,MAXFEV,EPSFCN,DIAG,MODE,FACTOR,NPRINT,
C      *           INFO,NFEV,NJEV,IPVT,QTF,WA1,WA2,WA3,WA4)
C       FNORM = ENORM(M,FVEC)
C       WRITE (NWRITE,1000) FNORM,NFEV,NJEV,INFO,(X(J),J=1,N)
C       STOP
C  1000 FORMAT (5X,' FINAL L2 NORM OF THE RESIDUALS',E15.7 //
C      *        5X,' NUMBER OF FUNCTION EVALUATIONS',I10 //
C      *        5X,' NUMBER OF JACOBIAN EVALUATIONS',I10 //
C      *        5X,' EXIT PARAMETER',16X,I10 //
C      *        5X,' FINAL APPROXIMATE SOLUTION' // 5X,3E15.7)
C       END
C       SUBROUTINE FCN(IFLAG,M,N,X,FVEC,DUM,IDUM)
C C     This is the form of the FCN routine if IOPT=1,
C C     that is, if the user does not calculate the Jacobian.
C       INTEGER M,N,IFLAG
C       REAL X(N),FVEC(M)
C       INTEGER I
C       REAL TMP1,TMP2,TMP3,TMP4
C       REAL Y(15)
C       DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
C      *     Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
C      *     /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1,
C      *      3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
C C
C       IF (IFLAG .NE. 0) GO TO 5
C C
C C     Insert print statements here when NPRINT is positive.
C C
C       RETURN
C     5 CONTINUE
C       DO 10 I = 1, M
C          TMP1 = I
C          TMP2 = 16 - I
C          TMP3 = TMP1
C          IF (I .GT. 8) TMP3 = TMP2
C          FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
C    10    CONTINUE
C       RETURN
C       END
C
C
C       Results obtained with different compilers or machines
C       may be slightly different.
C
C       FINAL L2 NORM OF THE RESIDUALS  0.9063596E-01
C
C       NUMBER OF FUNCTION EVALUATIONS        25
C
C       NUMBER OF JACOBIAN EVALUATIONS         0
C
C       EXIT PARAMETER                         1
C
C       FINAL APPROXIMATE SOLUTION
C
C        0.8241058E-01  0.1133037E+01  0.2343695E+01
C
C
C       For IOPT=2, FCN would be modified as follows to also
C       calculate the full Jacobian when IFLAG=2.
C
C       SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
C C
C C     This is the form of the FCN routine if IOPT=2,
C C     that is, if the user calculates the full Jacobian.
C C
C       INTEGER LDFJAC,M,N,IFLAG
C       REAL X(N),FVEC(M)
C       REAL FJAC(LDFJAC,N)
C       INTEGER I
C       REAL TMP1,TMP2,TMP3,TMP4
C       REAL Y(15)
C       DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
C      *     Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
C      *     /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1,
C      *      3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
C C
C       IF (IFLAG .NE. 0) GO TO 5
C C
C C     Insert print statements here when NPRINT is positive.
C C
C       RETURN
C     5 CONTINUE
C       IF(IFLAG.NE.1) GO TO 20
C       DO 10 I = 1, M
C          TMP1 = I
C          TMP2 = 16 - I
C          TMP3 = TMP1
C          IF (I .GT. 8) TMP3 = TMP2
C          FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
C    10    CONTINUE
C       RETURN
C C
C C     Below, calculate the full Jacobian.
C C
C    20    CONTINUE
C C
C       DO 30 I = 1, M
C          TMP1 = I
C          TMP2 = 16 - I
C          TMP3 = TMP1
C          IF (I .GT. 8) TMP3 = TMP2
C          TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2
C          FJAC(I,1) = -1.E0
C          FJAC(I,2) = TMP1*TMP2/TMP4
C          FJAC(I,3) = TMP1*TMP3/TMP4
C    30    CONTINUE
C       RETURN
C       END
C
C
C       For IOPT = 3, FJAC would be dimensioned as FJAC(3,3),
C         LDFJAC would be set to 3, and FCN would be written as
C         follows to calculate a row of the Jacobian when IFLAG=3.
C
C       SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
C C     This is the form of the FCN routine if IOPT=3,
C C     that is, if the user calculates the Jacobian row by row.
C       INTEGER M,N,IFLAG
C       REAL X(N),FVEC(M)
C       REAL FJAC(N)
C       INTEGER I
C       REAL TMP1,TMP2,TMP3,TMP4
C       REAL Y(15)
C       DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
C      *     Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
C      *     /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1,
C      *      3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
C C
C       IF (IFLAG .NE. 0) GO TO 5
C C
C C     Insert print statements here when NPRINT is positive.
C C
C       RETURN
C     5 CONTINUE
C       IF( IFLAG.NE.1) GO TO 20
C       DO 10 I = 1, M
C          TMP1 = I
C          TMP2 = 16 - I
C          TMP3 = TMP1
C          IF (I .GT. 8) TMP3 = TMP2
C          FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
C    10    CONTINUE
C       RETURN
C C
C C     Below, calculate the LDFJAC-th row of the Jacobian.
C C
C    20 CONTINUE
C
C       I = LDFJAC
C          TMP1 = I
C          TMP2 = 16 - I
C          TMP3 = TMP1
C          IF (I .GT. 8) TMP3 = TMP2
C          TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2
C          FJAC(1) = -1.E0
C          FJAC(2) = TMP1*TMP2/TMP4
C          FJAC(3) = TMP1*TMP3/TMP4
C       RETURN
C       END
C
C***REFERENCES  Jorge J. More, The Levenberg-Marquardt algorithm:
C                 implementation and theory.  In Numerical Analysis
C                 Proceedings (Dundee, June 28 - July 1, 1977, G. A.
C                 Watson, Editor), Lecture Notes in Mathematics 630,
C                 Springer-Verlag, 1978.
C***ROUTINES CALLED  CHKDER, ENORM, FDJAC3, LMPAR, QRFAC, R1MACH,
C                    RWUPDT, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   800301  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  SNLS1
      INTEGER IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV
      INTEGER IJUNK,NROW,IPVT(*)
      REAL FTOL,XTOL,GTOL,FACTOR,EPSFCN
      REAL X(*),FVEC(*),FJAC(LDFJAC,*),DIAG(*),QTF(*),WA1(*),WA2(*),
     1     WA3(*),WA4(*)
      LOGICAL SING
      EXTERNAL FCN
      INTEGER I,IFLAG,ITER,J,L,MODECH
      REAL ACTRED,DELTA,DIRDER,EPSMCH,FNORM,FNORM1,GNORM,ONE,PAR,
     1     PNORM,PRERED,P1,P5,P25,P75,P0001,RATIO,SUM,TEMP,TEMP1,
     2     TEMP2,XNORM,ZERO
      REAL R1MACH,ENORM,ERR,CHKLIM
      CHARACTER*8 XERN1
      CHARACTER*16 XERN3
C
      SAVE CHKLIM, ONE, P1, P5, P25, P75, P0001, ZERO
      DATA CHKLIM/.1E0/
      DATA ONE,P1,P5,P25,P75,P0001,ZERO
     1     /1.0E0,1.0E-1,5.0E-1,2.5E-1,7.5E-1,1.0E-4,0.0E0/
C
C***FIRST EXECUTABLE STATEMENT  SNLS1
      EPSMCH = R1MACH(4)
C
      INFO = 0
      IFLAG = 0
      NFEV = 0
      NJEV = 0
C
C     CHECK THE INPUT PARAMETERS FOR ERRORS.
C
      IF (IOPT .LT. 1 .OR. IOPT .GT. 3 .OR. N .LE. 0 .OR.
     1    M .LT. N .OR. LDFJAC .LT. N .OR. FTOL .LT. ZERO
     2    .OR. XTOL .LT. ZERO .OR. GTOL .LT. ZERO
     3    .OR. MAXFEV .LE. 0 .OR. FACTOR .LE. ZERO) GO TO 300
      IF (IOPT .LT. 3 .AND. LDFJAC .LT. M) GO TO 300
      IF (MODE .NE. 2) GO TO 20
      DO 10 J = 1, N
         IF (DIAG(J) .LE. ZERO) GO TO 300
   10    CONTINUE
   20 CONTINUE
C
C     EVALUATE THE FUNCTION AT THE STARTING POINT
C     AND CALCULATE ITS NORM.
C
      IFLAG = 1
      IJUNK = 1
      CALL FCN(IFLAG,M,N,X,FVEC,FJAC,IJUNK)
      NFEV = 1
      IF (IFLAG .LT. 0) GO TO 300
      FNORM = ENORM(M,FVEC)
C
C     INITIALIZE LEVENBERG-MARQUARDT PARAMETER AND ITERATION COUNTER.
C
      PAR = ZERO
      ITER = 1
C
C     BEGINNING OF THE OUTER LOOP.
C
   30 CONTINUE
C
C        IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.
C
         IF (NPRINT .LE. 0) GO TO 40
         IFLAG = 0
         IF (MOD(ITER-1,NPRINT) .EQ. 0)
     1      CALL FCN(IFLAG,M,N,X,FVEC,FJAC,IJUNK)
         IF (IFLAG .LT. 0) GO TO 300
   40    CONTINUE
C
C        CALCULATE THE JACOBIAN MATRIX.
C
      IF (IOPT .EQ. 3) GO TO 475
C
C     STORE THE FULL JACOBIAN USING M*N STORAGE
C
      IF (IOPT .EQ. 1) GO TO 410
C
C     THE USER SUPPLIES THE JACOBIAN
C
         IFLAG = 2
         CALL FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
         NJEV = NJEV + 1
C
C             ON THE FIRST ITERATION, CHECK THE USER SUPPLIED JACOBIAN
C
         IF (ITER .LE. 1) THEN
            IF (IFLAG .LT. 0) GO TO 300
C
C           GET THE INCREMENTED X-VALUES INTO WA1(*).
C
            MODECH = 1
            CALL CHKDER(M,N,X,FVEC,FJAC,LDFJAC,WA1,WA4,MODECH,ERR)
C
C           EVALUATE FUNCTION AT INCREMENTED VALUE AND PUT IN WA4(*).
C
            IFLAG = 1
            CALL FCN(IFLAG,M,N,WA1,WA4,FJAC,LDFJAC)
            NFEV = NFEV + 1
            IF(IFLAG .LT. 0) GO TO 300
            DO 350 I = 1, M
               MODECH = 2
               CALL CHKDER(1,N,X,FVEC(I),FJAC(I,1),LDFJAC,WA1,
     1              WA4(I),MODECH,ERR)
               IF (ERR .LT. CHKLIM) THEN
                  WRITE (XERN1, '(I8)') I
                  WRITE (XERN3, '(1PE15.6)') ERR
                  CALL XERMSG ('SLATEC', 'SNLS1', 'DERIVATIVE OF ' //
     *               'FUNCTION ' // XERN1 // ' MAY BE WRONG, ERR = ' //
     *               XERN3 // ' TOO CLOSE TO 0.', 7, 0)
               ENDIF
  350       CONTINUE
         ENDIF
C
         GO TO 420
C
C     THE CODE APPROXIMATES THE JACOBIAN
C
410      IFLAG = 1
         CALL FDJAC3(FCN,M,N,X,FVEC,FJAC,LDFJAC,IFLAG,EPSFCN,WA4)
         NFEV = NFEV + N
  420    IF (IFLAG .LT. 0) GO TO 300
C
C        COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.
C
         CALL QRFAC(M,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3)
C
C        FORM (Q TRANSPOSE)*FVEC AND STORE THE FIRST N COMPONENTS IN
C        QTF.
C
         DO 430 I = 1, M
            WA4(I) = FVEC(I)
  430         CONTINUE
         DO 470 J = 1, N
            IF (FJAC(J,J) .EQ. ZERO) GO TO 460
            SUM = ZERO
            DO 440 I = J, M
               SUM = SUM + FJAC(I,J)*WA4(I)
  440          CONTINUE
            TEMP = -SUM/FJAC(J,J)
            DO 450 I = J, M
               WA4(I) = WA4(I) + FJAC(I,J)*TEMP
  450          CONTINUE
  460       CONTINUE
            FJAC(J,J) = WA1(J)
            QTF(J) = WA4(J)
  470       CONTINUE
         GO TO 560
C
C        ACCUMULATE THE JACOBIAN BY ROWS IN ORDER TO SAVE STORAGE.
C        COMPUTE THE QR FACTORIZATION OF THE JACOBIAN MATRIX
C        CALCULATED ONE ROW AT A TIME, WHILE SIMULTANEOUSLY
C        FORMING (Q TRANSPOSE)*FVEC AND STORING THE FIRST
C        N COMPONENTS IN QTF.
C
  475    DO 490 J = 1, N
            QTF(J) = ZERO
            DO 480 I = 1, N
               FJAC(I,J) = ZERO
  480          CONTINUE
  490        CONTINUE
         DO 500 I = 1, M
            NROW = I
            IFLAG = 3
            CALL FCN(IFLAG,M,N,X,FVEC,WA3,NROW)
            IF (IFLAG .LT. 0) GO TO 300
C
C            ON THE FIRST ITERATION, CHECK THE USER SUPPLIED JACOBIAN.
C
            IF(ITER .GT. 1) GO TO 498
C
C            GET THE INCREMENTED X-VALUES INTO WA1(*).
C
            MODECH = 1
            CALL CHKDER(M,N,X,FVEC,FJAC,LDFJAC,WA1,WA4,MODECH,ERR)
C
C            EVALUATE AT INCREMENTED VALUES, IF NOT ALREADY EVALUATED.
C
            IF(I .NE. 1) GO TO 495
C
C            EVALUATE FUNCTION AT INCREMENTED VALUE AND PUT INTO WA4(*).
C
            IFLAG = 1
            CALL FCN(IFLAG,M,N,WA1,WA4,FJAC,NROW)
            NFEV = NFEV + 1
            IF(IFLAG .LT. 0) GO TO 300
495         CONTINUE
            MODECH = 2
            CALL CHKDER(1,N,X,FVEC(I),WA3,1,WA1,WA4(I),MODECH,ERR)
            IF (ERR .LT. CHKLIM) THEN
               WRITE (XERN1, '(I8)') I
               WRITE (XERN3, '(1PE15.6)') ERR
               CALL XERMSG ('SLATEC', 'SNLS1', 'DERIVATIVE OF FUNCTION '
     *            // XERN1 // ' MAY BE WRONG, ERR = ' // XERN3 //
     *            ' TOO CLOSE TO 0.', 7, 0)
            ENDIF
498         CONTINUE
C
            TEMP = FVEC(I)
            CALL RWUPDT(N,FJAC,LDFJAC,WA3,QTF,TEMP,WA1,WA2)
  500       CONTINUE
         NJEV = NJEV + 1
C
C        IF THE JACOBIAN IS RANK DEFICIENT, CALL QRFAC TO
C        REORDER ITS COLUMNS AND UPDATE THE COMPONENTS OF QTF.
C
         SING = .FALSE.
         DO 510 J = 1, N
            IF (FJAC(J,J) .EQ. ZERO) SING = .TRUE.
            IPVT(J) = J
            WA2(J) = ENORM(J,FJAC(1,J))
  510       CONTINUE
         IF (.NOT.SING) GO TO 560
         CALL QRFAC(N,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3)
         DO 550 J = 1, N
            IF (FJAC(J,J) .EQ. ZERO) GO TO 540
            SUM = ZERO
            DO 520 I = J, N
               SUM = SUM + FJAC(I,J)*QTF(I)
  520         CONTINUE
            TEMP = -SUM/FJAC(J,J)
            DO 530 I = J, N
               QTF(I) = QTF(I) + FJAC(I,J)*TEMP
  530          CONTINUE
  540       CONTINUE
            FJAC(J,J) = WA1(J)
  550       CONTINUE
  560    CONTINUE
C
C        ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING
C        TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.
C
         IF (ITER .NE. 1) GO TO 80
         IF (MODE .EQ. 2) GO TO 60
         DO 50 J = 1, N
            DIAG(J) = WA2(J)
            IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE
   50       CONTINUE
   60    CONTINUE
C
C        ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X
C        AND INITIALIZE THE STEP BOUND DELTA.
C
         DO 70 J = 1, N
            WA3(J) = DIAG(J)*X(J)
   70       CONTINUE
         XNORM = ENORM(N,WA3)
         DELTA = FACTOR*XNORM
         IF (DELTA .EQ. ZERO) DELTA = FACTOR
   80    CONTINUE
C
C        COMPUTE THE NORM OF THE SCALED GRADIENT.
C
         GNORM = ZERO
         IF (FNORM .EQ. ZERO) GO TO 170
         DO 160 J = 1, N
            L = IPVT(J)
            IF (WA2(L) .EQ. ZERO) GO TO 150
            SUM = ZERO
            DO 140 I = 1, J
               SUM = SUM + FJAC(I,J)*(QTF(I)/FNORM)
  140          CONTINUE
            GNORM = MAX(GNORM,ABS(SUM/WA2(L)))
  150       CONTINUE
  160       CONTINUE
  170    CONTINUE
C
C        TEST FOR CONVERGENCE OF THE GRADIENT NORM.
C
         IF (GNORM .LE. GTOL) INFO = 4
         IF (INFO .NE. 0) GO TO 300
C
C        RESCALE IF NECESSARY.
C
         IF (MODE .EQ. 2) GO TO 190
         DO 180 J = 1, N
            DIAG(J) = MAX(DIAG(J),WA2(J))
  180       CONTINUE
  190    CONTINUE
C
C        BEGINNING OF THE INNER LOOP.
C
  200    CONTINUE
C
C           DETERMINE THE LEVENBERG-MARQUARDT PARAMETER.
C
            CALL LMPAR(N,FJAC,LDFJAC,IPVT,DIAG,QTF,DELTA,PAR,WA1,WA2,
     1                 WA3,WA4)
C
C           STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.
C
            DO 210 J = 1, N
               WA1(J) = -WA1(J)
               WA2(J) = X(J) + WA1(J)
               WA3(J) = DIAG(J)*WA1(J)
  210          CONTINUE
            PNORM = ENORM(N,WA3)
C
C           ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.
C
            IF (ITER .EQ. 1) DELTA = MIN(DELTA,PNORM)
C
C           EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.
C
            IFLAG = 1
            CALL FCN(IFLAG,M,N,WA2,WA4,FJAC,IJUNK)
            NFEV = NFEV + 1
            IF (IFLAG .LT. 0) GO TO 300
            FNORM1 = ENORM(M,WA4)
C
C           COMPUTE THE SCALED ACTUAL REDUCTION.
C
            ACTRED = -ONE
            IF (P1*FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2
C
C           COMPUTE THE SCALED PREDICTED REDUCTION AND
C           THE SCALED DIRECTIONAL DERIVATIVE.
C
            DO 230 J = 1, N
               WA3(J) = ZERO
               L = IPVT(J)
               TEMP = WA1(L)
               DO 220 I = 1, J
                  WA3(I) = WA3(I) + FJAC(I,J)*TEMP
  220             CONTINUE
  230          CONTINUE
            TEMP1 = ENORM(N,WA3)/FNORM
            TEMP2 = (SQRT(PAR)*PNORM)/FNORM
            PRERED = TEMP1**2 + TEMP2**2/P5
            DIRDER = -(TEMP1**2 + TEMP2**2)
C
C           COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED
C           REDUCTION.
C
            RATIO = ZERO
            IF (PRERED .NE. ZERO) RATIO = ACTRED/PRERED
C
C           UPDATE THE STEP BOUND.
C
            IF (RATIO .GT. P25) GO TO 240
               IF (ACTRED .GE. ZERO) TEMP = P5
               IF (ACTRED .LT. ZERO)
     1            TEMP = P5*DIRDER/(DIRDER + P5*ACTRED)
               IF (P1*FNORM1 .GE. FNORM .OR. TEMP .LT. P1) TEMP = P1
               DELTA = TEMP*MIN(DELTA,PNORM/P1)
               PAR = PAR/TEMP
               GO TO 260
  240       CONTINUE
               IF (PAR .NE. ZERO .AND. RATIO .LT. P75) GO TO 250
               DELTA = PNORM/P5
               PAR = P5*PAR
  250          CONTINUE
  260       CONTINUE
C
C           TEST FOR SUCCESSFUL ITERATION.
C
            IF (RATIO .LT. P0001) GO TO 290
C
C           SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.
C
            DO 270 J = 1, N
               X(J) = WA2(J)
               WA2(J) = DIAG(J)*X(J)
  270          CONTINUE
            DO 280 I = 1, M
               FVEC(I) = WA4(I)
  280          CONTINUE
            XNORM = ENORM(N,WA2)
            FNORM = FNORM1
            ITER = ITER + 1
  290       CONTINUE
C
C           TESTS FOR CONVERGENCE.
C
            IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL
     1          .AND. P5*RATIO .LE. ONE) INFO = 1
            IF (DELTA .LE. XTOL*XNORM) INFO = 2
            IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL
     1          .AND. P5*RATIO .LE. ONE .AND. INFO .EQ. 2) INFO = 3
            IF (INFO .NE. 0) GO TO 300
C
C           TESTS FOR TERMINATION AND STRINGENT TOLERANCES.
C
            IF (NFEV .GE. MAXFEV) INFO = 5
            IF (ABS(ACTRED) .LE. EPSMCH .AND. PRERED .LE. EPSMCH
     1          .AND. P5*RATIO .LE. ONE) INFO = 6
            IF (DELTA .LE. EPSMCH*XNORM) INFO = 7
            IF (GNORM .LE. EPSMCH) INFO = 8
            IF (INFO .NE. 0) GO TO 300
C
C           END OF THE INNER LOOP. REPEAT IF ITERATION UNSUCCESSFUL.
C
            IF (RATIO .LT. P0001) GO TO 200
C
C        END OF THE OUTER LOOP.
C
         GO TO 30
  300 CONTINUE
C
C     TERMINATION, EITHER NORMAL OR USER IMPOSED.
C
      IF (IFLAG .LT. 0) INFO = IFLAG
      IFLAG = 0
      IF (NPRINT .GT. 0) CALL FCN(IFLAG,M,N,X,FVEC,FJAC,IJUNK)
      IF (INFO .LT. 0) CALL XERMSG ('SLATEC', 'SNLS1',
     +   'EXECUTION TERMINATED BECAUSE USER SET IFLAG NEGATIVE.', 1, 1)
      IF (INFO .EQ. 0) CALL XERMSG ('SLATEC', 'SNLS1',
     +   'INVALID INPUT PARAMETER.', 2, 1)
      IF (INFO .EQ. 4) CALL XERMSG ('SLATEC', 'SNLS1',
     +   'THIRD CONVERGENCE CONDITION, CHECK RESULTS BEFORE ACCEPTING.',
     +   1, 1)
      IF (INFO .EQ. 5) CALL XERMSG ('SLATEC', 'SNLS1',
     +   'TOO MANY FUNCTION EVALUATIONS.', 9, 1)
      IF (INFO .GE. 6) CALL XERMSG ('SLATEC', 'SNLS1',
     +   'TOLERANCES TOO SMALL, NO FURTHER IMPROVEMENT POSSIBLE.', 3, 1)
      RETURN
C
C     LAST CARD OF SUBROUTINE SNLS1.
C
      END