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
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)
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.
30 C 2. Subroutine and Type Statements.
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
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)
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.
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.
60 C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
61 C INTEGER IFLAG,LDFJAC,M,N
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.
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.
74 C If IFLAG=1, calculate the functions at X and return
75 C this vector in FVEC.
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).
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).
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.
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
110 C M is a positive integer input variable set to the number of
113 C N is a positive integer input variable set to the number of
114 C variables. N must not exceed M.
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.
120 C FVEC is an output array of length M which contains the functions
121 C evaluated at the output X.
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
129 C P *(JAC *JAC)*P = R *R,
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.
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
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.
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
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.
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.
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).
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
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
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.
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.
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.
201 C INFO = 0 improper input parameters.
203 C INFO = 1 both actual and predicted relative reductions in the
204 C sum of squares are at most FTOL.
206 C INFO = 2 relative error between two consecutive iterates is
209 C INFO = 3 conditions for INFO = 1 and INFO = 2 both hold.
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.
214 C INFO = 5 number of calls to FCN for function evaluation
215 C has reached MAXFEV.
217 C INFO = 6 FTOL is too small. No further reduction in the sum
218 C of squares is possible.
220 C INFO = 7 XTOL is too small. No further improvement in the
221 C approximate solution X is possible.
223 C INFO = 8 GTOL is too small. FVEC is orthogonal to the
224 C columns of the Jacobian to machine precision.
226 C Sections 4 and 5 contain more details about INFO.
228 C NFEV is an integer output variable set to the number of calls to
229 C FCN for function evaluation.
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.
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
244 C QTF is an output array of length N which contains the first N
245 C elements of the vector (Q transpose)*FVEC.
247 C WA1, WA2, and WA3 are work arrays of length N.
249 C WA4 is a work array of length M.
252 C 4. Successful Completion.
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.
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.
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
274 C ENORM(FVEC) .LE. (1+FTOL)*ENORM(FVECS),
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
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
288 C ENORM(D*(X-XSOL)) .LE. XTOL*ENORM(D*XSOL).
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
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.
310 C 5. Unsuccessful Completion.
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.
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
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.
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
337 C 6. Characteristics of the Algorithm.
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.
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.
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.
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
374 C X(1) + U(I)/(V(I)*X(2) + W(I)*X(3)), I = 1, 15
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),
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
384 C Y(I) - (X(1) + U(I)/(V(I)*X(2) + W(I)*X(3))).
390 C C Driver for SNLS1 example.
392 C INTEGER J,IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,
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)
406 C C The following starting values provide a rough fit.
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.
418 C FTOL = SQRT(R1MACH(4))
419 C XTOL = SQRT(R1MACH(4))
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)
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)
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.
446 C REAL TMP1,TMP2,TMP3,TMP4
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/
453 C IF (IFLAG .NE. 0) GO TO 5
455 C C Insert print statements here when NPRINT is positive.
463 C IF (I .GT. 8) TMP3 = TMP2
464 C FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
470 C Results obtained with different compilers or machines
471 C may be slightly different.
473 C FINAL L2 NORM OF THE RESIDUALS 0.9063596E-01
475 C NUMBER OF FUNCTION EVALUATIONS 25
477 C NUMBER OF JACOBIAN EVALUATIONS 0
481 C FINAL APPROXIMATE SOLUTION
483 C 0.8241058E-01 0.1133037E+01 0.2343695E+01
486 C For IOPT=2, FCN would be modified as follows to also
487 C calculate the full Jacobian when IFLAG=2.
489 C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
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.
494 C INTEGER LDFJAC,M,N,IFLAG
496 C REAL FJAC(LDFJAC,N)
498 C REAL TMP1,TMP2,TMP3,TMP4
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/
505 C IF (IFLAG .NE. 0) GO TO 5
507 C C Insert print statements here when NPRINT is positive.
511 C IF(IFLAG.NE.1) GO TO 20
516 C IF (I .GT. 8) TMP3 = TMP2
517 C FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
521 C C Below, calculate the full Jacobian.
529 C IF (I .GT. 8) TMP3 = TMP2
530 C TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2
532 C FJAC(I,2) = TMP1*TMP2/TMP4
533 C FJAC(I,3) = TMP1*TMP3/TMP4
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.
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.
550 C REAL TMP1,TMP2,TMP3,TMP4
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/
557 C IF (IFLAG .NE. 0) GO TO 5
559 C C Insert print statements here when NPRINT is positive.
563 C IF( IFLAG.NE.1) GO TO 20
568 C IF (I .GT. 8) TMP3 = TMP2
569 C FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
573 C C Below, calculate the LDFJAC-th row of the Jacobian.
581 C IF (I .GT. 8) TMP3 = TMP2
582 C TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2
584 C FJAC(2) = TMP1*TMP2/TMP4
585 C FJAC(3) = TMP1*TMP3/TMP4
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,
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(*),
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,
616 REAL R1MACH,ENORM,ERR,CHKLIM
620 SAVE CHKLIM, ONE, P1, P5, P25, P75, P0001, ZERO
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/
625 C***FIRST EXECUTABLE STATEMENT SNLS1
633 C CHECK THE INPUT PARAMETERS FOR ERRORS.
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
642 IF (DIAG(J) .LE. ZERO) GO TO 300
646 C EVALUATE THE FUNCTION AT THE STARTING POINT
647 C AND CALCULATE ITS NORM.
651 CALL FCN(IFLAG,M,N,X,FVEC,FJAC,IJUNK)
653 IF (IFLAG .LT. 0) GO TO 300
654 FNORM = ENORM(M,FVEC)
656 C INITIALIZE LEVENBERG-MARQUARDT PARAMETER AND ITERATION COUNTER.
661 C BEGINNING OF THE OUTER LOOP.
665 C IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.
667 IF (NPRINT .LE. 0) GO TO 40
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
674 C CALCULATE THE JACOBIAN MATRIX.
676 IF (IOPT .EQ. 3) GO TO 475
678 C STORE THE FULL JACOBIAN USING M*N STORAGE
680 IF (IOPT .EQ. 1) GO TO 410
682 C THE USER SUPPLIES THE JACOBIAN
685 CALL FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
688 C ON THE FIRST ITERATION, CHECK THE USER SUPPLIED JACOBIAN
690 IF (ITER .LE. 1) THEN
691 IF (IFLAG .LT. 0) GO TO 300
693 C GET THE INCREMENTED X-VALUES INTO WA1(*).
696 CALL CHKDER(M,N,X,FVEC,FJAC,LDFJAC,WA1,WA4,MODECH,ERR)
698 C EVALUATE FUNCTION AT INCREMENTED VALUE AND PUT IN WA4(*).
701 CALL FCN(IFLAG,M,N,WA1,WA4,FJAC,LDFJAC)
703 IF(IFLAG .LT. 0) GO TO 300
706 CALL CHKDER(1,N,X,FVEC(I),FJAC(I,1),LDFJAC,WA1,
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)
720 C THE CODE APPROXIMATES THE JACOBIAN
723 CALL FDJAC3(FCN,M,N,X,FVEC,FJAC,LDFJAC,IFLAG,EPSFCN,WA4)
725 420 IF (IFLAG .LT. 0) GO TO 300
727 C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.
729 CALL QRFAC(M,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3)
731 C FORM (Q TRANSPOSE)*FVEC AND STORE THE FIRST N COMPONENTS IN
738 IF (FJAC(J,J) .EQ. ZERO) GO TO 460
741 SUM = SUM + FJAC(I,J)*WA4(I)
743 TEMP = -SUM/FJAC(J,J)
745 WA4(I) = WA4(I) + FJAC(I,J)*TEMP
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.
768 CALL FCN(IFLAG,M,N,X,FVEC,WA3,NROW)
769 IF (IFLAG .LT. 0) GO TO 300
771 C ON THE FIRST ITERATION, CHECK THE USER SUPPLIED JACOBIAN.
773 IF(ITER .GT. 1) GO TO 498
775 C GET THE INCREMENTED X-VALUES INTO WA1(*).
778 CALL CHKDER(M,N,X,FVEC,FJAC,LDFJAC,WA1,WA4,MODECH,ERR)
780 C EVALUATE AT INCREMENTED VALUES, IF NOT ALREADY EVALUATED.
782 IF(I .NE. 1) GO TO 495
784 C EVALUATE FUNCTION AT INCREMENTED VALUE AND PUT INTO WA4(*).
787 CALL FCN(IFLAG,M,N,WA1,WA4,FJAC,NROW)
789 IF(IFLAG .LT. 0) GO TO 300
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)
803 CALL RWUPDT(N,FJAC,LDFJAC,WA3,QTF,TEMP,WA1,WA2)
807 C IF THE JACOBIAN IS RANK DEFICIENT, CALL QRFAC TO
808 C REORDER ITS COLUMNS AND UPDATE THE COMPONENTS OF QTF.
812 IF (FJAC(J,J) .EQ. ZERO) SING = .TRUE.
814 WA2(J) = ENORM(J,FJAC(1,J))
816 IF (.NOT.SING) GO TO 560
817 CALL QRFAC(N,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3)
819 IF (FJAC(J,J) .EQ. ZERO) GO TO 540
822 SUM = SUM + FJAC(I,J)*QTF(I)
824 TEMP = -SUM/FJAC(J,J)
826 QTF(I) = QTF(I) + FJAC(I,J)*TEMP
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.
836 IF (ITER .NE. 1) GO TO 80
837 IF (MODE .EQ. 2) GO TO 60
840 IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE
844 C ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X
845 C AND INITIALIZE THE STEP BOUND DELTA.
848 WA3(J) = DIAG(J)*X(J)
852 IF (DELTA .EQ. ZERO) DELTA = FACTOR
855 C COMPUTE THE NORM OF THE SCALED GRADIENT.
858 IF (FNORM .EQ. ZERO) GO TO 170
861 IF (WA2(L) .EQ. ZERO) GO TO 150
864 SUM = SUM + FJAC(I,J)*(QTF(I)/FNORM)
866 GNORM = MAX(GNORM,ABS(SUM/WA2(L)))
871 C TEST FOR CONVERGENCE OF THE GRADIENT NORM.
873 IF (GNORM .LE. GTOL) INFO = 4
874 IF (INFO .NE. 0) GO TO 300
876 C RESCALE IF NECESSARY.
878 IF (MODE .EQ. 2) GO TO 190
880 DIAG(J) = MAX(DIAG(J),WA2(J))
884 C BEGINNING OF THE INNER LOOP.
888 C DETERMINE THE LEVENBERG-MARQUARDT PARAMETER.
890 CALL LMPAR(N,FJAC,LDFJAC,IPVT,DIAG,QTF,DELTA,PAR,WA1,WA2,
893 C STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.
897 WA2(J) = X(J) + WA1(J)
898 WA3(J) = DIAG(J)*WA1(J)
902 C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.
904 IF (ITER .EQ. 1) DELTA = MIN(DELTA,PNORM)
906 C EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.
909 CALL FCN(IFLAG,M,N,WA2,WA4,FJAC,IJUNK)
911 IF (IFLAG .LT. 0) GO TO 300
912 FNORM1 = ENORM(M,WA4)
914 C COMPUTE THE SCALED ACTUAL REDUCTION.
917 IF (P1*FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2
919 C COMPUTE THE SCALED PREDICTED REDUCTION AND
920 C THE SCALED DIRECTIONAL DERIVATIVE.
927 WA3(I) = WA3(I) + FJAC(I,J)*TEMP
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)
935 C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED
939 IF (PRERED .NE. ZERO) RATIO = ACTRED/PRERED
941 C UPDATE THE STEP BOUND.
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)
952 IF (PAR .NE. ZERO .AND. RATIO .LT. P75) GO TO 250
958 C TEST FOR SUCCESSFUL ITERATION.
960 IF (RATIO .LT. P0001) GO TO 290
962 C SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.
966 WA2(J) = DIAG(J)*X(J)
976 C TESTS FOR CONVERGENCE.
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
985 C TESTS FOR TERMINATION AND STRINGENT TOLERANCES.
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
994 C END OF THE INNER LOOP. REPEAT IF ITERATION UNSUCCESSFUL.
996 IF (RATIO .LT. P0001) GO TO 200
998 C END OF THE OUTER LOOP.
1003 C TERMINATION, EITHER NORMAL OR USER IMPOSED.
1005 IF (IFLAG .LT. 0) INFO = IFLAG
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.',
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)
1021 C LAST CARD OF SUBROUTINE SNLS1.