binaries/src/disembl/Tisean_3.0.1/source_f/slatec/snls1.f

   1 *DECK SNLS1
   2       SUBROUTINE SNLS1 (FCN, IOPT, M, N, X, FVEC, FJAC, LDFJAC, FTOL,
   3      +   XTOL, GTOL, MAXFEV, EPSFCN, DIAG, MODE, FACTOR, NPRINT, INFO,
   4      +   NFEV, NJEV, IPVT, QTF, WA1, WA2, WA3, WA4)
   5 C***BEGIN PROLOGUE  SNLS1
   6 C***PURPOSE  Minimize the sum of the squares of M nonlinear functions
   7 C            in N variables by a modification of the Levenberg-Marquardt
   8 C            algorithm.
   9 C***LIBRARY   SLATEC
  10 C***CATEGORY  K1B1A1, K1B1A2
  11 C***TYPE      SINGLE PRECISION (SNLS1-S, DNLS1-D)
  12 C***KEYWORDS  LEVENBERG-MARQUARDT, NONLINEAR DATA FITTING,
  13 C             NONLINEAR LEAST SQUARES
  14 C***AUTHOR  Hiebert, K. L., (SNLA)
  15 C***DESCRIPTION
  16 C
  17 C 1. Purpose.
  18 C
  19 C       The purpose of SNLS1 is to minimize the sum of the squares of M
  20 C       nonlinear functions in N variables by a modification of the
  21 C       Levenberg-Marquardt algorithm.  The user must provide a subrou-
  22 C       tine which calculates the functions.  The user has the option
  23 C       of how the Jacobian will be supplied.  The user can supply the
  24 C       full Jacobian, or the rows of the Jacobian (to avoid storing
  25 C       the full Jacobian), or let the code approximate the Jacobian by
  26 C       forward-differencing.   This code is the combination of the
  27 C       MINPACK codes (Argonne) LMDER, LMDIF, and LMSTR.
  28 C
  29 C
  30 C 2. Subroutine and Type Statements.
  31 C
  32 C       SUBROUTINE SNLS1(FCN,IOPT,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,
  33 C      *                 GTOL,MAXFEV,EPSFCN,DIAG,MODE,FACTOR,NPRINT,INFO
  34 C      *                 ,NFEV,NJEV,IPVT,QTF,WA1,WA2,WA3,WA4)
  35 C       INTEGER IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV
  36 C       INTEGER IPVT(N)
  37 C       REAL FTOL,XTOL,GTOL,EPSFCN,FACTOR
  38 C       REAL X(N),FVEC(M),FJAC(LDFJAC,N),DIAG(N),QTF(N),
  39 C      *     WA1(N),WA2(N),WA3(N),WA4(M)
  40 C
  41 C
  42 C 3. Parameters.
  43 C
  44 C       Parameters designated as input parameters must be specified on
  45 C       entry to SNLS1 and are not changed on exit, while parameters
  46 C       designated as output parameters need not be specified on entry
  47 C       and are set to appropriate values on exit from SNLS1.
  48 C
  49 C       FCN is the name of the user-supplied subroutine which calculates
  50 C         the functions.  If the user wants to supply the Jacobian
  51 C         (IOPT=2 or 3), then FCN must be written to calculate the
  52 C         Jacobian, as well as the functions.  See the explanation
  53 C         of the IOPT argument below.
  54 C         If the user wants the iterates printed (NPRINT positive), then
  55 C         FCN must do the printing.  See the explanation of NPRINT
  56 C         below.  FCN must be declared in an EXTERNAL statement in the
  57 C         calling program and should be written as follows.
  58 C
  59 C
  60 C         SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
  61 C         INTEGER IFLAG,LDFJAC,M,N
  62 C         REAL X(N),FVEC(M)
  63 C         ----------
  64 C         FJAC and LDFJAC may be ignored     , if IOPT=1.
  65 C         REAL FJAC(LDFJAC,N)                , if IOPT=2.
  66 C         REAL FJAC(N)                       , if IOPT=3.
  67 C         ----------
  68 C           If IFLAG=0, the values in X and FVEC are available
  69 C           for printing.  See the explanation of NPRINT below.
  70 C           IFLAG will never be zero unless NPRINT is positive.
  71 C           The values of X and FVEC must not be changed.
  72 C         RETURN
  73 C         ----------
  74 C           If IFLAG=1, calculate the functions at X and return
  75 C           this vector in FVEC.
  76 C         RETURN
  77 C         ----------
  78 C           If IFLAG=2, calculate the full Jacobian at X and return
  79 C           this matrix in FJAC.  Note that IFLAG will never be 2 unless
  80 C           IOPT=2.  FVEC contains the function values at X and must
  81 C           not be altered.  FJAC(I,J) must be set to the derivative
  82 C           of FVEC(I) with respect to X(J).
  83 C         RETURN
  84 C         ----------
  85 C           If IFLAG=3, calculate the LDFJAC-th row of the Jacobian
  86 C           and return this vector in FJAC.  Note that IFLAG will
  87 C           never be 3 unless IOPT=3.  FVEC contains the function
  88 C           values at X and must not be altered.  FJAC(J) must be
  89 C           set to the derivative of FVEC(LDFJAC) with respect to X(J).
  90 C         RETURN
  91 C         ----------
  92 C         END
  93 C
  94 C
  95 C         The value of IFLAG should not be changed by FCN unless the
  96 C         user wants to terminate execution of SNLS1.  In this case, set
  97 C         IFLAG to a negative integer.
  98 C
  99 C
 100 C       IOPT is an input variable which specifies how the Jacobian will
 101 C         be calculated.  If IOPT=2 or 3, then the user must supply the
 102 C         Jacobian, as well as the function values, through the
 103 C         subroutine FCN.  If IOPT=2, the user supplies the full
 104 C         Jacobian with one call to FCN.  If IOPT=3, the user supplies
 105 C         one row of the Jacobian with each call.  (In this manner,
 106 C         storage can be saved because the full Jacobian is not stored.)
 107 C         If IOPT=1, the code will approximate the Jacobian by forward
 108 C         differencing.
 109 C
 110 C       M is a positive integer input variable set to the number of
 111 C         functions.
 112 C
 113 C       N is a positive integer input variable set to the number of
 114 C         variables.  N must not exceed M.
 115 C
 116 C       X is an array of length N.  On input, X must contain an initial
 117 C         estimate of the solution vector.  On output, X contains the
 118 C         final estimate of the solution vector.
 119 C
 120 C       FVEC is an output array of length M which contains the functions
 121 C         evaluated at the output X.
 122 C
 123 C       FJAC is an output array.  For IOPT=1 and 2, FJAC is an M by N
 124 C         array.  For IOPT=3, FJAC is an N by N array.  The upper N by N
 125 C         submatrix of FJAC contains an upper triangular matrix R with
 126 C         diagonal elements of nonincreasing magnitude such that
 127 C
 128 C                T     T           T
 129 C               P *(JAC *JAC)*P = R *R,
 130 C
 131 C         where P is a permutation matrix and JAC is the final calcu-
 132 C         lated Jacobian.  Column J of P is column IPVT(J) (see below)
 133 C         of the identity matrix.  The lower part of FJAC contains
 134 C         information generated during the computation of R.
 135 C
 136 C       LDFJAC is a positive integer input variable which specifies
 137 C         the leading dimension of the array FJAC.  For IOPT=1 and 2,
 138 C         LDFJAC must not be less than M.  For IOPT=3, LDFJAC must not
 139 C         be less than N.
 140 C
 141 C       FTOL is a non-negative input variable.  Termination occurs when
 142 C         both the actual and predicted relative reductions in the sum
 143 C         of squares are at most FTOL.  Therefore, FTOL measures the
 144 C         relative error desired in the sum of squares.  Section 4 con-
 145 C         tains more details about FTOL.
 146 C
 147 C       XTOL is a non-negative input variable.  Termination occurs when
 148 C         the relative error between two consecutive iterates is at most
 149 C         XTOL.  Therefore, XTOL measures the relative error desired in
 150 C         the approximate solution.  Section 4 contains more details
 151 C         about XTOL.
 152 C
 153 C       GTOL is a non-negative input variable.  Termination occurs when
 154 C         the cosine of the angle between FVEC and any column of the
 155 C         Jacobian is at most GTOL in absolute value.  Therefore, GTOL
 156 C         measures the orthogonality desired between the function vector
 157 C         and the columns of the Jacobian.  Section 4 contains more
 158 C         details about GTOL.
 159 C
 160 C       MAXFEV is a positive integer input variable.  Termination occurs
 161 C         when the number of calls to FCN to evaluate the functions
 162 C         has reached MAXFEV.
 163 C
 164 C       EPSFCN is an input variable used in determining a suitable step
 165 C         for the forward-difference approximation.  This approximation
 166 C         assumes that the relative errors in the functions are of the
 167 C         order of EPSFCN.  If EPSFCN is less than the machine preci-
 168 C         sion, it is assumed that the relative errors in the functions
 169 C         are of the order of the machine precision.  If IOPT=2 or 3,
 170 C         then EPSFCN can be ignored (treat it as a dummy argument).
 171 C
 172 C       DIAG is an array of length N.  If MODE = 1 (see below), DIAG is
 173 C         internally set.  If MODE = 2, DIAG must contain positive
 174 C         entries that serve as implicit (multiplicative) scale factors
 175 C         for the variables.
 176 C
 177 C       MODE is an integer input variable.  If MODE = 1, the variables
 178 C         will be scaled internally.  If MODE = 2, the scaling is speci-
 179 C         fied by the input DIAG.  Other values of MODE are equivalent
 180 C         to MODE = 1.
 181 C
 182 C       FACTOR is a positive input variable used in determining the ini-
 183 C         tial step bound.  This bound is set to the product of FACTOR
 184 C         and the Euclidean norm of DIAG*X if nonzero, or else to FACTOR
 185 C         itself.  In most cases FACTOR should lie in the interval
 186 C         (.1,100.).  100. is a generally recommended value.
 187 C
 188 C       NPRINT is an integer input variable that enables controlled
 189 C         printing of iterates if it is positive.  In this case, FCN is
 190 C         called with IFLAG = 0 at the beginning of the first iteration
 191 C         and every NPRINT iterations thereafter and immediately prior
 192 C         to return, with X and FVEC available for printing. Appropriate
 193 C         print statements must be added to FCN (see example) and
 194 C         FVEC should not be altered.  If NPRINT is not positive, no
 195 C         special calls to FCN with IFLAG = 0 are made.
 196 C
 197 C       INFO is an integer output variable.  If the user has terminated
 198 C         execution, INFO is set to the (negative) value of IFLAG.  See
 199 C         description of FCN and JAC. Otherwise, INFO is set as follows.
 200 C
 201 C         INFO = 0  improper input parameters.
 202 C
 203 C         INFO = 1  both actual and predicted relative reductions in the
 204 C                   sum of squares are at most FTOL.
 205 C
 206 C         INFO = 2  relative error between two consecutive iterates is
 207 C                   at most XTOL.
 208 C
 209 C         INFO = 3  conditions for INFO = 1 and INFO = 2 both hold.
 210 C
 211 C         INFO = 4  the cosine of the angle between FVEC and any column
 212 C                   of the Jacobian is at most GTOL in absolute value.
 213 C
 214 C         INFO = 5  number of calls to FCN for function evaluation
 215 C                   has reached MAXFEV.
 216 C
 217 C         INFO = 6  FTOL is too small.  No further reduction in the sum
 218 C                   of squares is possible.
 219 C
 220 C         INFO = 7  XTOL is too small.  No further improvement in the
 221 C                   approximate solution X is possible.
 222 C
 223 C         INFO = 8  GTOL is too small.  FVEC is orthogonal to the
 224 C                   columns of the Jacobian to machine precision.
 225 C
 226 C         Sections 4 and 5 contain more details about INFO.
 227 C
 228 C       NFEV is an integer output variable set to the number of calls to
 229 C         FCN for function evaluation.
 230 C
 231 C       NJEV is an integer output variable set to the number of
 232 C         evaluations of the full Jacobian.  If IOPT=2, only one call to
 233 C         FCN is required for each evaluation of the full Jacobian.
 234 C         If IOPT=3, the M calls to FCN are required.
 235 C         If IOPT=1, then NJEV is set to zero.
 236 C
 237 C       IPVT is an integer output array of length N.  IPVT defines a
 238 C         permutation matrix P such that JAC*P = Q*R, where JAC is the
 239 C         final calculated Jacobian, Q is orthogonal (not stored), and R
 240 C         is upper triangular with diagonal elements of nonincreasing
 241 C         magnitude.  Column J of P is column IPVT(J) of the identity
 242 C         matrix.
 243 C
 244 C       QTF is an output array of length N which contains the first N
 245 C         elements of the vector (Q transpose)*FVEC.
 246 C
 247 C       WA1, WA2, and WA3 are work arrays of length N.
 248 C
 249 C       WA4 is a work array of length M.
 250 C
 251 C
 252 C 4. Successful Completion.
 253 C
 254 C       The accuracy of SNLS1 is controlled by the convergence parame-
 255 C       ters FTOL, XTOL, and GTOL.  These parameters are used in tests
 256 C       which make three types of comparisons between the approximation
 257 C       X and a solution XSOL.  SNLS1 terminates when any of the tests
 258 C       is satisfied.  If any of the convergence parameters is less than
 259 C       the machine precision (as defined by the function R1MACH(4)),
 260 C       then SNLS1 only attempts to satisfy the test defined by the
 261 C       machine precision.  Further progress is not usually possible.
 262 C
 263 C       The tests assume that the functions are reasonably well behaved,
 264 C       and, if the Jacobian is supplied by the user, that the functions
 265 C       and the Jacobian are coded consistently.  If these conditions
 266 C       are not satisfied, then SNLS1 may incorrectly indicate conver-
 267 C       gence.  If the Jacobian is coded correctly or IOPT=1,
 268 C       then the validity of the answer can be checked, for example, by
 269 C       rerunning SNLS1 with tighter tolerances.
 270 C
 271 C       First Convergence Test.  If ENORM(Z) denotes the Euclidean norm
 272 C         of a vector Z, then this test attempts to guarantee that
 273 C
 274 C               ENORM(FVEC) .LE. (1+FTOL)*ENORM(FVECS),
 275 C
 276 C         where FVECS denotes the functions evaluated at XSOL.  If this
 277 C         condition is satisfied with FTOL = 10**(-K), then the final
 278 C         residual norm ENORM(FVEC) has K significant decimal digits and
 279 C         INFO is set to 1 (or to 3 if the second test is also satis-
 280 C         fied).  Unless high precision solutions are required, the
 281 C         recommended value for FTOL is the square root of the machine
 282 C         precision.
 283 C
 284 C       Second Convergence Test.  If D is the diagonal matrix whose
 285 C         entries are defined by the array DIAG, then this test attempts
 286 C         to guarantee that
 287 C
 288 C               ENORM(D*(X-XSOL)) .LE. XTOL*ENORM(D*XSOL).
 289 C
 290 C         If this condition is satisfied with XTOL = 10**(-K), then the
 291 C         larger components of D*X have K significant decimal digits and
 292 C         INFO is set to 2 (or to 3 if the first test is also satis-
 293 C         fied).  There is a danger that the smaller components of D*X
 294 C         may have large relative errors, but if MODE = 1, then the
 295 C         accuracy of the components of X is usually related to their
 296 C         sensitivity.  Unless high precision solutions are required,
 297 C         the recommended value for XTOL is the square root of the
 298 C         machine precision.
 299 C
 300 C       Third Convergence Test.  This test is satisfied when the cosine
 301 C         of the angle between FVEC and any column of the Jacobian at X
 302 C         is at most GTOL in absolute value.  There is no clear rela-
 303 C         tionship between this test and the accuracy of SNLS1, and
 304 C         furthermore, the test is equally well satisfied at other crit-
 305 C         ical points, namely maximizers and saddle points.  Therefore,
 306 C         termination caused by this test (INFO = 4) should be examined
 307 C         carefully.  The recommended value for GTOL is zero.
 308 C
 309 C
 310 C 5. Unsuccessful Completion.
 311 C
 312 C       Unsuccessful termination of SNLS1 can be due to improper input
 313 C       parameters, arithmetic interrupts, or an excessive number of
 314 C       function evaluations.
 315 C
 316 C       Improper Input Parameters.  INFO is set to 0 if IOPT .LT. 1
 317 C         or IOPT .GT. 3, or N .LE. 0, or M .LT. N, or for IOPT=1 or 2
 318 C         LDFJAC .LT. M, or for IOPT=3 LDFJAC .LT. N, or FTOL .LT. 0.E0,
 319 C         or XTOL .LT. 0.E0, or GTOL .LT. 0.E0, or MAXFEV .LE. 0, or
 320 C         FACTOR .LE. 0.E0.
 321 C
 322 C       Arithmetic Interrupts.  If these interrupts occur in the FCN
 323 C         subroutine during an early stage of the computation, they may
 324 C         be caused by an unacceptable choice of X by SNLS1.  In this
 325 C         case, it may be possible to remedy the situation by rerunning
 326 C         SNLS1 with a smaller value of FACTOR.
 327 C
 328 C       Excessive Number of Function Evaluations.  A reasonable value
 329 C         for MAXFEV is 100*(N+1) for IOPT=2 or 3 and 200*(N+1) for
 330 C         IOPT=1.  If the number of calls to FCN reaches MAXFEV, then
 331 C         this indicates that the routine is converging very slowly
 332 C         as measured by the progress of FVEC, and INFO is set to 5.
 333 C         In this case, it may be helpful to restart SNLS1 with MODE
 334 C         set to 1.
 335 C
 336 C
 337 C 6. Characteristics of the Algorithm.
 338 C
 339 C       SNLS1 is a modification of the Levenberg-Marquardt algorithm.
 340 C       Two of its main characteristics involve the proper use of
 341 C       implicitly scaled variables (if MODE = 1) and an optimal choice
 342 C       for the correction.  The use of implicitly scaled variables
 343 C       achieves scale invariance of SNLS1 and limits the size of the
 344 C       correction in any direction where the functions are changing
 345 C       rapidly.  The optimal choice of the correction guarantees (under
 346 C       reasonable conditions) global convergence from starting points
 347 C       far from the solution and a fast rate of convergence for
 348 C       problems with small residuals.
 349 C
 350 C       Timing.  The time required by SNLS1 to solve a given problem
 351 C         depends on M and N, the behavior of the functions, the accu-
 352 C         racy requested, and the starting point.  The number of arith-
 353 C         metic operations needed by SNLS1 is about N**3 to process each
 354 C         evaluation of the functions (call to FCN) and to process each
 355 C         evaluation of the Jacobian it takes M*N**2 for IOPT=2 (one
 356 C         call to FCN), M*N**2 for IOPT=1 (N calls to FCN) and
 357 C         1.5*M*N**2 for IOPT=3 (M calls to FCN).  Unless FCN
 358 C         can be evaluated quickly, the timing of SNLS1 will be
 359 C         strongly influenced by the time spent in FCN.
 360 C
 361 C       Storage.  SNLS1 requires (M*N + 2*M + 6*N) for IOPT=1 or 2 and
 362 C         (N**2 + 2*M + 6*N) for IOPT=3 single precision storage
 363 C         locations and N integer storage locations, in addition to
 364 C         the storage required by the program.  There are no internally
 365 C         declared storage arrays.
 366 C
 367 C *Long Description:
 368 C
 369 C 7. Example.
 370 C
 371 C       The problem is to determine the values of X(1), X(2), and X(3)
 372 C       which provide the best fit (in the least squares sense) of
 373 C
 374 C             X(1) + U(I)/(V(I)*X(2) + W(I)*X(3)),  I = 1, 15
 375 C
 376 C       to the data
 377 C
 378 C             Y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,
 379 C                  0.37,0.58,0.73,0.96,1.34,2.10,4.39),
 380 C
 381 C       where U(I) = I, V(I) = 16 - I, and W(I) = MIN(U(I),V(I)).  The
 382 C       I-th component of FVEC is thus defined by
 383 C
 384 C             Y(I) - (X(1) + U(I)/(V(I)*X(2) + W(I)*X(3))).
 385 C
 386 C       **********
 387 C
 388 C       PROGRAM TEST
 389 C C
 390 C C     Driver for SNLS1 example.
 391 C C
 392 C       INTEGER J,IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,
 393 C      *        NWRITE
 394 C       INTEGER IPVT(3)
 395 C       REAL FTOL,XTOL,GTOL,FACTOR,FNORM,EPSFCN
 396 C       REAL X(3),FVEC(15),FJAC(15,3),DIAG(3),QTF(3),
 397 C      *     WA1(3),WA2(3),WA3(3),WA4(15)
 398 C       REAL ENORM,R1MACH
 399 C       EXTERNAL FCN
 400 C       DATA NWRITE /6/
 401 C C
 402 C       IOPT = 1
 403 C       M = 15
 404 C       N = 3
 405 C C
 406 C C     The following starting values provide a rough fit.
 407 C C
 408 C       X(1) = 1.E0
 409 C       X(2) = 1.E0
 410 C       X(3) = 1.E0
 411 C C
 412 C       LDFJAC = 15
 413 C C
 414 C C     Set FTOL and XTOL to the square root of the machine precision
 415 C C     and GTOL to zero.  Unless high precision solutions are
 416 C C     required, these are the recommended settings.
 417 C C
 418 C       FTOL = SQRT(R1MACH(4))
 419 C       XTOL = SQRT(R1MACH(4))
 420 C       GTOL = 0.E0
 421 C C
 422 C       MAXFEV = 400
 423 C       EPSFCN = 0.0
 424 C       MODE = 1
 425 C       FACTOR = 1.E2
 426 C       NPRINT = 0
 427 C C
 428 C       CALL SNLS1(FCN,IOPT,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,
 429 C      *           GTOL,MAXFEV,EPSFCN,DIAG,MODE,FACTOR,NPRINT,
 430 C      *           INFO,NFEV,NJEV,IPVT,QTF,WA1,WA2,WA3,WA4)
 431 C       FNORM = ENORM(M,FVEC)
 432 C       WRITE (NWRITE,1000) FNORM,NFEV,NJEV,INFO,(X(J),J=1,N)
 433 C       STOP
 434 C  1000 FORMAT (5X,' FINAL L2 NORM OF THE RESIDUALS',E15.7 //
 435 C      *        5X,' NUMBER OF FUNCTION EVALUATIONS',I10 //
 436 C      *        5X,' NUMBER OF JACOBIAN EVALUATIONS',I10 //
 437 C      *        5X,' EXIT PARAMETER',16X,I10 //
 438 C      *        5X,' FINAL APPROXIMATE SOLUTION' // 5X,3E15.7)
 439 C       END
 440 C       SUBROUTINE FCN(IFLAG,M,N,X,FVEC,DUM,IDUM)
 441 C C     This is the form of the FCN routine if IOPT=1,
 442 C C     that is, if the user does not calculate the Jacobian.
 443 C       INTEGER M,N,IFLAG
 444 C       REAL X(N),FVEC(M)
 445 C       INTEGER I
 446 C       REAL TMP1,TMP2,TMP3,TMP4
 447 C       REAL Y(15)
 448 C       DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
 449 C      *     Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
 450 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,
 451 C      *      3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
 452 C C
 453 C       IF (IFLAG .NE. 0) GO TO 5
 454 C C
 455 C C     Insert print statements here when NPRINT is positive.
 456 C C
 457 C       RETURN
 458 C     5 CONTINUE
 459 C       DO 10 I = 1, M
 460 C          TMP1 = I
 461 C          TMP2 = 16 - I
 462 C          TMP3 = TMP1
 463 C          IF (I .GT. 8) TMP3 = TMP2
 464 C          FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
 465 C    10    CONTINUE
 466 C       RETURN
 467 C       END
 468 C
 469 C
 470 C       Results obtained with different compilers or machines
 471 C       may be slightly different.
 472 C
 473 C       FINAL L2 NORM OF THE RESIDUALS  0.9063596E-01
 474 C
 475 C       NUMBER OF FUNCTION EVALUATIONS        25
 476 C
 477 C       NUMBER OF JACOBIAN EVALUATIONS         0
 478 C
 479 C       EXIT PARAMETER                         1
 480 C
 481 C       FINAL APPROXIMATE SOLUTION
 482 C
 483 C        0.8241058E-01  0.1133037E+01  0.2343695E+01
 484 C
 485 C
 486 C       For IOPT=2, FCN would be modified as follows to also
 487 C       calculate the full Jacobian when IFLAG=2.
 488 C
 489 C       SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
 490 C C
 491 C C     This is the form of the FCN routine if IOPT=2,
 492 C C     that is, if the user calculates the full Jacobian.
 493 C C
 494 C       INTEGER LDFJAC,M,N,IFLAG
 495 C       REAL X(N),FVEC(M)
 496 C       REAL FJAC(LDFJAC,N)
 497 C       INTEGER I
 498 C       REAL TMP1,TMP2,TMP3,TMP4
 499 C       REAL Y(15)
 500 C       DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
 501 C      *     Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
 502 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,
 503 C      *      3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
 504 C C
 505 C       IF (IFLAG .NE. 0) GO TO 5
 506 C C
 507 C C     Insert print statements here when NPRINT is positive.
 508 C C
 509 C       RETURN
 510 C     5 CONTINUE
 511 C       IF(IFLAG.NE.1) GO TO 20
 512 C       DO 10 I = 1, M
 513 C          TMP1 = I
 514 C          TMP2 = 16 - I
 515 C          TMP3 = TMP1
 516 C          IF (I .GT. 8) TMP3 = TMP2
 517 C          FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
 518 C    10    CONTINUE
 519 C       RETURN
 520 C C
 521 C C     Below, calculate the full Jacobian.
 522 C C
 523 C    20    CONTINUE
 524 C C
 525 C       DO 30 I = 1, M
 526 C          TMP1 = I
 527 C          TMP2 = 16 - I
 528 C          TMP3 = TMP1
 529 C          IF (I .GT. 8) TMP3 = TMP2
 530 C          TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2
 531 C          FJAC(I,1) = -1.E0
 532 C          FJAC(I,2) = TMP1*TMP2/TMP4
 533 C          FJAC(I,3) = TMP1*TMP3/TMP4
 534 C    30    CONTINUE
 535 C       RETURN
 536 C       END
 537 C
 538 C
 539 C       For IOPT = 3, FJAC would be dimensioned as FJAC(3,3),
 540 C         LDFJAC would be set to 3, and FCN would be written as
 541 C         follows to calculate a row of the Jacobian when IFLAG=3.
 542 C
 543 C       SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
 544 C C     This is the form of the FCN routine if IOPT=3,
 545 C C     that is, if the user calculates the Jacobian row by row.
 546 C       INTEGER M,N,IFLAG
 547 C       REAL X(N),FVEC(M)
 548 C       REAL FJAC(N)
 549 C       INTEGER I
 550 C       REAL TMP1,TMP2,TMP3,TMP4
 551 C       REAL Y(15)
 552 C       DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
 553 C      *     Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
 554 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,
 555 C      *      3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
 556 C C
 557 C       IF (IFLAG .NE. 0) GO TO 5
 558 C C
 559 C C     Insert print statements here when NPRINT is positive.
 560 C C
 561 C       RETURN
 562 C     5 CONTINUE
 563 C       IF( IFLAG.NE.1) GO TO 20
 564 C       DO 10 I = 1, M
 565 C          TMP1 = I
 566 C          TMP2 = 16 - I
 567 C          TMP3 = TMP1
 568 C          IF (I .GT. 8) TMP3 = TMP2
 569 C          FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
 570 C    10    CONTINUE
 571 C       RETURN
 572 C C
 573 C C     Below, calculate the LDFJAC-th row of the Jacobian.
 574 C C
 575 C    20 CONTINUE
 576 C
 577 C       I = LDFJAC
 578 C          TMP1 = I
 579 C          TMP2 = 16 - I
 580 C          TMP3 = TMP1
 581 C          IF (I .GT. 8) TMP3 = TMP2
 582 C          TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2
 583 C          FJAC(1) = -1.E0
 584 C          FJAC(2) = TMP1*TMP2/TMP4
 585 C          FJAC(3) = TMP1*TMP3/TMP4
 586 C       RETURN
 587 C       END
 588 C
 589 C***REFERENCES  Jorge J. More, The Levenberg-Marquardt algorithm:
 590 C                 implementation and theory.  In Numerical Analysis
 591 C                 Proceedings (Dundee, June 28 - July 1, 1977, G. A.
 592 C                 Watson, Editor), Lecture Notes in Mathematics 630,
 593 C                 Springer-Verlag, 1978.
 594 C***ROUTINES CALLED  CHKDER, ENORM, FDJAC3, LMPAR, QRFAC, R1MACH,
 595 C                    RWUPDT, XERMSG
 596 C***REVISION HISTORY  (YYMMDD)
 597 C   800301  DATE WRITTEN
 598 C   890531  Changed all specific intrinsics to generic.  (WRB)
 599 C   890531  REVISION DATE from Version 3.2
 600 C   891214  Prologue converted to Version 4.0 format.  (BAB)
 601 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 602 C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
 603 C   920501  Reformatted the REFERENCES section.  (WRB)
 604 C***END PROLOGUE  SNLS1
 605       INTEGER IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV
 606       INTEGER IJUNK,NROW,IPVT(*)
 607       REAL FTOL,XTOL,GTOL,FACTOR,EPSFCN
 608       REAL X(*),FVEC(*),FJAC(LDFJAC,*),DIAG(*),QTF(*),WA1(*),WA2(*),
 609      1     WA3(*),WA4(*)
 610       LOGICAL SING
 611       EXTERNAL FCN
 612       INTEGER I,IFLAG,ITER,J,L,MODECH
 613       REAL ACTRED,DELTA,DIRDER,EPSMCH,FNORM,FNORM1,GNORM,ONE,PAR,
 614      1     PNORM,PRERED,P1,P5,P25,P75,P0001,RATIO,SUM,TEMP,TEMP1,
 615      2     TEMP2,XNORM,ZERO
 616       REAL R1MACH,ENORM,ERR,CHKLIM
 617       CHARACTER*8 XERN1
 618       CHARACTER*16 XERN3
 619 C
 620       SAVE CHKLIM, ONE, P1, P5, P25, P75, P0001, ZERO
 621       DATA CHKLIM/.1E0/
 622       DATA ONE,P1,P5,P25,P75,P0001,ZERO
 623      1     /1.0E0,1.0E-1,5.0E-1,2.5E-1,7.5E-1,1.0E-4,0.0E0/
 624 C
 625 C***FIRST EXECUTABLE STATEMENT  SNLS1
 626       EPSMCH = R1MACH(4)
 627 C
 628       INFO = 0
 629       IFLAG = 0
 630       NFEV = 0
 631       NJEV = 0
 632 C
 633 C     CHECK THE INPUT PARAMETERS FOR ERRORS.
 634 C
 635       IF (IOPT .LT. 1 .OR. IOPT .GT. 3 .OR. N .LE. 0 .OR.
 636      1    M .LT. N .OR. LDFJAC .LT. N .OR. FTOL .LT. ZERO
 637      2    .OR. XTOL .LT. ZERO .OR. GTOL .LT. ZERO
 638      3    .OR. MAXFEV .LE. 0 .OR. FACTOR .LE. ZERO) GO TO 300
 639       IF (IOPT .LT. 3 .AND. LDFJAC .LT. M) GO TO 300
 640       IF (MODE .NE. 2) GO TO 20
 641       DO 10 J = 1, N
 642          IF (DIAG(J) .LE. ZERO) GO TO 300
 643    10    CONTINUE
 644    20 CONTINUE
 645 C
 646 C     EVALUATE THE FUNCTION AT THE STARTING POINT
 647 C     AND CALCULATE ITS NORM.
 648 C
 649       IFLAG = 1
 650       IJUNK = 1
 651       CALL FCN(IFLAG,M,N,X,FVEC,FJAC,IJUNK)
 652       NFEV = 1
 653       IF (IFLAG .LT. 0) GO TO 300
 654       FNORM = ENORM(M,FVEC)
 655 C
 656 C     INITIALIZE LEVENBERG-MARQUARDT PARAMETER AND ITERATION COUNTER.
 657 C
 658       PAR = ZERO
 659       ITER = 1
 660 C
 661 C     BEGINNING OF THE OUTER LOOP.
 662 C
 663    30 CONTINUE
 664 C
 665 C        IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.
 666 C
 667          IF (NPRINT .LE. 0) GO TO 40
 668          IFLAG = 0
 669          IF (MOD(ITER-1,NPRINT) .EQ. 0)
 670      1      CALL FCN(IFLAG,M,N,X,FVEC,FJAC,IJUNK)
 671          IF (IFLAG .LT. 0) GO TO 300
 672    40    CONTINUE
 673 C
 674 C        CALCULATE THE JACOBIAN MATRIX.
 675 C
 676       IF (IOPT .EQ. 3) GO TO 475
 677 C
 678 C     STORE THE FULL JACOBIAN USING M*N STORAGE
 679 C
 680       IF (IOPT .EQ. 1) GO TO 410
 681 C
 682 C     THE USER SUPPLIES THE JACOBIAN
 683 C
 684          IFLAG = 2
 685          CALL FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
 686          NJEV = NJEV + 1
 687 C
 688 C             ON THE FIRST ITERATION, CHECK THE USER SUPPLIED JACOBIAN
 689 C
 690          IF (ITER .LE. 1) THEN
 691             IF (IFLAG .LT. 0) GO TO 300
 692 C
 693 C           GET THE INCREMENTED X-VALUES INTO WA1(*).
 694 C
 695             MODECH = 1
 696             CALL CHKDER(M,N,X,FVEC,FJAC,LDFJAC,WA1,WA4,MODECH,ERR)
 697 C
 698 C           EVALUATE FUNCTION AT INCREMENTED VALUE AND PUT IN WA4(*).
 699 C
 700             IFLAG = 1
 701             CALL FCN(IFLAG,M,N,WA1,WA4,FJAC,LDFJAC)
 702             NFEV = NFEV + 1
 703             IF(IFLAG .LT. 0) GO TO 300
 704             DO 350 I = 1, M
 705                MODECH = 2
 706                CALL CHKDER(1,N,X,FVEC(I),FJAC(I,1),LDFJAC,WA1,
 707      1              WA4(I),MODECH,ERR)
 708                IF (ERR .LT. CHKLIM) THEN
 709                   WRITE (XERN1, '(I8)') I
 710                   WRITE (XERN3, '(1PE15.6)') ERR
 711                   CALL XERMSG ('SLATEC', 'SNLS1', 'DERIVATIVE OF ' //
 712      *               'FUNCTION ' // XERN1 // ' MAY BE WRONG, ERR = ' //
 713      *               XERN3 // ' TOO CLOSE TO 0.', 7, 0)
 714                ENDIF
 715   350       CONTINUE
 716          ENDIF
 717 C
 718          GO TO 420
 719 C
 720 C     THE CODE APPROXIMATES THE JACOBIAN
 721 C
 722 410      IFLAG = 1
 723          CALL FDJAC3(FCN,M,N,X,FVEC,FJAC,LDFJAC,IFLAG,EPSFCN,WA4)
 724          NFEV = NFEV + N
 725   420    IF (IFLAG .LT. 0) GO TO 300
 726 C
 727 C        COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.
 728 C
 729          CALL QRFAC(M,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3)
 730 C
 731 C        FORM (Q TRANSPOSE)*FVEC AND STORE THE FIRST N COMPONENTS IN
 732 C        QTF.
 733 C
 734          DO 430 I = 1, M
 735             WA4(I) = FVEC(I)
 736   430         CONTINUE
 737          DO 470 J = 1, N
 738             IF (FJAC(J,J) .EQ. ZERO) GO TO 460
 739             SUM = ZERO
 740             DO 440 I = J, M
 741                SUM = SUM + FJAC(I,J)*WA4(I)
 742   440          CONTINUE
 743             TEMP = -SUM/FJAC(J,J)
 744             DO 450 I = J, M
 745                WA4(I) = WA4(I) + FJAC(I,J)*TEMP
 746   450          CONTINUE
 747   460       CONTINUE
 748             FJAC(J,J) = WA1(J)
 749             QTF(J) = WA4(J)
 750   470       CONTINUE
 751          GO TO 560
 752 C
 753 C        ACCUMULATE THE JACOBIAN BY ROWS IN ORDER TO SAVE STORAGE.
 754 C        COMPUTE THE QR FACTORIZATION OF THE JACOBIAN MATRIX
 755 C        CALCULATED ONE ROW AT A TIME, WHILE SIMULTANEOUSLY
 756 C        FORMING (Q TRANSPOSE)*FVEC AND STORING THE FIRST
 757 C        N COMPONENTS IN QTF.
 758 C
 759   475    DO 490 J = 1, N
 760             QTF(J) = ZERO
 761             DO 480 I = 1, N
 762                FJAC(I,J) = ZERO
 763   480          CONTINUE
 764   490        CONTINUE
 765          DO 500 I = 1, M
 766             NROW = I
 767             IFLAG = 3
 768             CALL FCN(IFLAG,M,N,X,FVEC,WA3,NROW)
 769             IF (IFLAG .LT. 0) GO TO 300
 770 C
 771 C            ON THE FIRST ITERATION, CHECK THE USER SUPPLIED JACOBIAN.
 772 C
 773             IF(ITER .GT. 1) GO TO 498
 774 C
 775 C            GET THE INCREMENTED X-VALUES INTO WA1(*).
 776 C
 777             MODECH = 1
 778             CALL CHKDER(M,N,X,FVEC,FJAC,LDFJAC,WA1,WA4,MODECH,ERR)
 779 C
 780 C            EVALUATE AT INCREMENTED VALUES, IF NOT ALREADY EVALUATED.
 781 C
 782             IF(I .NE. 1) GO TO 495
 783 C
 784 C            EVALUATE FUNCTION AT INCREMENTED VALUE AND PUT INTO WA4(*).
 785 C
 786             IFLAG = 1
 787             CALL FCN(IFLAG,M,N,WA1,WA4,FJAC,NROW)
 788             NFEV = NFEV + 1
 789             IF(IFLAG .LT. 0) GO TO 300
 790 495         CONTINUE
 791             MODECH = 2
 792             CALL CHKDER(1,N,X,FVEC(I),WA3,1,WA1,WA4(I),MODECH,ERR)
 793             IF (ERR .LT. CHKLIM) THEN
 794                WRITE (XERN1, '(I8)') I
 795                WRITE (XERN3, '(1PE15.6)') ERR
 796                CALL XERMSG ('SLATEC', 'SNLS1', 'DERIVATIVE OF FUNCTION '
 797      *            // XERN1 // ' MAY BE WRONG, ERR = ' // XERN3 //
 798      *            ' TOO CLOSE TO 0.', 7, 0)
 799             ENDIF
 800 498         CONTINUE
 801 C
 802             TEMP = FVEC(I)
 803             CALL RWUPDT(N,FJAC,LDFJAC,WA3,QTF,TEMP,WA1,WA2)
 804   500       CONTINUE
 805          NJEV = NJEV + 1
 806 C
 807 C        IF THE JACOBIAN IS RANK DEFICIENT, CALL QRFAC TO
 808 C        REORDER ITS COLUMNS AND UPDATE THE COMPONENTS OF QTF.
 809 C
 810          SING = .FALSE.
 811          DO 510 J = 1, N
 812             IF (FJAC(J,J) .EQ. ZERO) SING = .TRUE.
 813             IPVT(J) = J
 814             WA2(J) = ENORM(J,FJAC(1,J))
 815   510       CONTINUE
 816          IF (.NOT.SING) GO TO 560
 817          CALL QRFAC(N,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3)
 818          DO 550 J = 1, N
 819             IF (FJAC(J,J) .EQ. ZERO) GO TO 540
 820             SUM = ZERO
 821             DO 520 I = J, N
 822                SUM = SUM + FJAC(I,J)*QTF(I)
 823   520         CONTINUE
 824             TEMP = -SUM/FJAC(J,J)
 825             DO 530 I = J, N
 826                QTF(I) = QTF(I) + FJAC(I,J)*TEMP
 827   530          CONTINUE
 828   540       CONTINUE
 829             FJAC(J,J) = WA1(J)
 830   550       CONTINUE
 831   560    CONTINUE
 832 C
 833 C        ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING
 834 C        TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.
 835 C
 836          IF (ITER .NE. 1) GO TO 80
 837          IF (MODE .EQ. 2) GO TO 60
 838          DO 50 J = 1, N
 839             DIAG(J) = WA2(J)
 840             IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE
 841    50       CONTINUE
 842    60    CONTINUE
 843 C
 844 C        ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X
 845 C        AND INITIALIZE THE STEP BOUND DELTA.
 846 C
 847          DO 70 J = 1, N
 848             WA3(J) = DIAG(J)*X(J)
 849    70       CONTINUE
 850          XNORM = ENORM(N,WA3)
 851          DELTA = FACTOR*XNORM
 852          IF (DELTA .EQ. ZERO) DELTA = FACTOR
 853    80    CONTINUE
 854 C
 855 C        COMPUTE THE NORM OF THE SCALED GRADIENT.
 856 C
 857          GNORM = ZERO
 858          IF (FNORM .EQ. ZERO) GO TO 170
 859          DO 160 J = 1, N
 860             L = IPVT(J)
 861             IF (WA2(L) .EQ. ZERO) GO TO 150
 862             SUM = ZERO
 863             DO 140 I = 1, J
 864                SUM = SUM + FJAC(I,J)*(QTF(I)/FNORM)
 865   140          CONTINUE
 866             GNORM = MAX(GNORM,ABS(SUM/WA2(L)))
 867   150       CONTINUE
 868   160       CONTINUE
 869   170    CONTINUE
 870 C
 871 C        TEST FOR CONVERGENCE OF THE GRADIENT NORM.
 872 C
 873          IF (GNORM .LE. GTOL) INFO = 4
 874          IF (INFO .NE. 0) GO TO 300
 875 C
 876 C        RESCALE IF NECESSARY.
 877 C
 878          IF (MODE .EQ. 2) GO TO 190
 879          DO 180 J = 1, N
 880             DIAG(J) = MAX(DIAG(J),WA2(J))
 881   180       CONTINUE
 882   190    CONTINUE
 883 C
 884 C        BEGINNING OF THE INNER LOOP.
 885 C
 886   200    CONTINUE
 887 C
 888 C           DETERMINE THE LEVENBERG-MARQUARDT PARAMETER.
 889 C
 890             CALL LMPAR(N,FJAC,LDFJAC,IPVT,DIAG,QTF,DELTA,PAR,WA1,WA2,
 891      1                 WA3,WA4)
 892 C
 893 C           STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.
 894 C
 895             DO 210 J = 1, N
 896                WA1(J) = -WA1(J)
 897                WA2(J) = X(J) + WA1(J)
 898                WA3(J) = DIAG(J)*WA1(J)
 899   210          CONTINUE
 900             PNORM = ENORM(N,WA3)
 901 C
 902 C           ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.
 903 C
 904             IF (ITER .EQ. 1) DELTA = MIN(DELTA,PNORM)
 905 C
 906 C           EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.
 907 C
 908             IFLAG = 1
 909             CALL FCN(IFLAG,M,N,WA2,WA4,FJAC,IJUNK)
 910             NFEV = NFEV + 1
 911             IF (IFLAG .LT. 0) GO TO 300
 912             FNORM1 = ENORM(M,WA4)
 913 C
 914 C           COMPUTE THE SCALED ACTUAL REDUCTION.
 915 C
 916             ACTRED = -ONE
 917             IF (P1*FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2
 918 C
 919 C           COMPUTE THE SCALED PREDICTED REDUCTION AND
 920 C           THE SCALED DIRECTIONAL DERIVATIVE.
 921 C
 922             DO 230 J = 1, N
 923                WA3(J) = ZERO
 924                L = IPVT(J)
 925                TEMP = WA1(L)
 926                DO 220 I = 1, J
 927                   WA3(I) = WA3(I) + FJAC(I,J)*TEMP
 928   220             CONTINUE
 929   230          CONTINUE
 930             TEMP1 = ENORM(N,WA3)/FNORM
 931             TEMP2 = (SQRT(PAR)*PNORM)/FNORM
 932             PRERED = TEMP1**2 + TEMP2**2/P5
 933             DIRDER = -(TEMP1**2 + TEMP2**2)
 934 C
 935 C           COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED
 936 C           REDUCTION.
 937 C
 938             RATIO = ZERO
 939             IF (PRERED .NE. ZERO) RATIO = ACTRED/PRERED
 940 C
 941 C           UPDATE THE STEP BOUND.
 942 C
 943             IF (RATIO .GT. P25) GO TO 240
 944                IF (ACTRED .GE. ZERO) TEMP = P5
 945                IF (ACTRED .LT. ZERO)
 946      1            TEMP = P5*DIRDER/(DIRDER + P5*ACTRED)
 947                IF (P1*FNORM1 .GE. FNORM .OR. TEMP .LT. P1) TEMP = P1
 948                DELTA = TEMP*MIN(DELTA,PNORM/P1)
 949                PAR = PAR/TEMP
 950                GO TO 260
 951   240       CONTINUE
 952                IF (PAR .NE. ZERO .AND. RATIO .LT. P75) GO TO 250
 953                DELTA = PNORM/P5
 954                PAR = P5*PAR
 955   250          CONTINUE
 956   260       CONTINUE
 957 C
 958 C           TEST FOR SUCCESSFUL ITERATION.
 959 C
 960             IF (RATIO .LT. P0001) GO TO 290
 961 C
 962 C           SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.
 963 C
 964             DO 270 J = 1, N
 965                X(J) = WA2(J)
 966                WA2(J) = DIAG(J)*X(J)
 967   270          CONTINUE
 968             DO 280 I = 1, M
 969                FVEC(I) = WA4(I)
 970   280          CONTINUE
 971             XNORM = ENORM(N,WA2)
 972             FNORM = FNORM1
 973             ITER = ITER + 1
 974   290       CONTINUE
 975 C
 976 C           TESTS FOR CONVERGENCE.
 977 C
 978             IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL
 979      1          .AND. P5*RATIO .LE. ONE) INFO = 1
 980             IF (DELTA .LE. XTOL*XNORM) INFO = 2
 981             IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL
 982      1          .AND. P5*RATIO .LE. ONE .AND. INFO .EQ. 2) INFO = 3
 983             IF (INFO .NE. 0) GO TO 300
 984 C
 985 C           TESTS FOR TERMINATION AND STRINGENT TOLERANCES.
 986 C
 987             IF (NFEV .GE. MAXFEV) INFO = 5
 988             IF (ABS(ACTRED) .LE. EPSMCH .AND. PRERED .LE. EPSMCH
 989      1          .AND. P5*RATIO .LE. ONE) INFO = 6
 990             IF (DELTA .LE. EPSMCH*XNORM) INFO = 7
 991             IF (GNORM .LE. EPSMCH) INFO = 8
 992             IF (INFO .NE. 0) GO TO 300
 993 C
 994 C           END OF THE INNER LOOP. REPEAT IF ITERATION UNSUCCESSFUL.
 995 C
 996             IF (RATIO .LT. P0001) GO TO 200
 997 C
 998 C        END OF THE OUTER LOOP.
 999 C
1000          GO TO 30
1001   300 CONTINUE
1002 C
1003 C     TERMINATION, EITHER NORMAL OR USER IMPOSED.
1004 C
1005       IF (IFLAG .LT. 0) INFO = IFLAG
1006       IFLAG = 0
1007       IF (NPRINT .GT. 0) CALL FCN(IFLAG,M,N,X,FVEC,FJAC,IJUNK)
1008       IF (INFO .LT. 0) CALL XERMSG ('SLATEC', 'SNLS1',
1009      +   'EXECUTION TERMINATED BECAUSE USER SET IFLAG NEGATIVE.', 1, 1)
1010       IF (INFO .EQ. 0) CALL XERMSG ('SLATEC', 'SNLS1',
1011      +   'INVALID INPUT PARAMETER.', 2, 1)
1012       IF (INFO .EQ. 4) CALL XERMSG ('SLATEC', 'SNLS1',
1013      +   'THIRD CONVERGENCE CONDITION, CHECK RESULTS BEFORE ACCEPTING.',
1014      +   1, 1)
1015       IF (INFO .EQ. 5) CALL XERMSG ('SLATEC', 'SNLS1',
1016      +   'TOO MANY FUNCTION EVALUATIONS.', 9, 1)
1017       IF (INFO .GE. 6) CALL XERMSG ('SLATEC', 'SNLS1',
1018      +   'TOLERANCES TOO SMALL, NO FURTHER IMPROVEMENT POSSIBLE.', 3, 1)
1019       RETURN
1020 C
1021 C     LAST CARD OF SUBROUTINE SNLS1.
1022 C
1023       END