1 2 SUBROUTINE SNLS1(FCN,IOPT,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL, 3 1 MAXFEV,EPSFCN,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,IPVT,QTF, 4 2 WA1,WA2,WA3,WA4) 5C***BEGIN PROLOGUE SNLS1 6C***DATE WRITTEN 800301 (YYMMDD) 7C***REVISION DATE 840405 (YYMMDD) 8C***CATEGORY NO. K1B1A1,K1B1A2 9C***KEYWORDS LEVENBERG-MARQUARDT,NONLINEAR DATA FITTING, 10C NONLINEAR LEAST SQUARES 11C***AUTHOR HIEBERT, K. L., (SNLA) 12C***PURPOSE SNLS1 minimizes the sum of the squares of M nonlinear 13C functions in N variables by a modification of the 14C Levenberg-Marquardt algorithm. This code is a combination 15C of the MINPACK codes (Argonne) LMDER, LMDIF, and LMSTR. 16C***DESCRIPTION 17C 18C 1. Purpose. 19C 20C The purpose of SNLS1 is to minimize the sum of the squares of M 21C nonlinear functions in N variables by a modification of the 22C Levenberg-Marquardt algorithm. The user must provide a subrou- 23C tine which calculates the functions. The user has the option 24C of how the Jacobian will be supplied. The user can supply the 25C full Jacobian, or the rows of the Jacobian (to avoid storing 26C the full Jacobian), or let the code approximate the Jacobian by 27C forward-differencing. This code is the combination of the 28C MINPACK codes (Argonne) LMDER, LMDIF, and LMSTR. 29C 30C 31C 2. Subroutine and Type Statements. 32C 33C SUBROUTINE SNLS1(FCN,IOPT,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL, 34C * GTOL,MAXFEV,EPSFCN,DIAG,MODE,FACTOR,NPRINT,INFO 35C * ,NFEV,NJEV,IPVT,QTF,WA1,WA2,WA3,WA4) 36C INTEGER IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV 37C INTEGER IPVT(N) 38C REAL FTOL,XTOL,GTOL,EPSFCN,FACTOR 39C REAL X(N),FVEC(M),FJAC(LDFJAC,N),DIAG(N),QTF(N), 40C * WA1(N),WA2(N),WA3(N),WA4(M) 41C 42C 43C 3. Parameters. 44C 45C Parameters designated as input parameters must be specified on 46C entry to SNLS1 and are not changed on exit, while parameters 47C designated as output parameters need not be specified on entry 48C and are set to appropriate values on exit from SNLS1. 49C 50C FCN is the name of the user-supplied subroutine which calculates 51C the functions. If the user wants to supply the Jacobian 52C (IOPT=2 or 3), then FCN must be written to calculate the 53C Jacobian, as well as the functions. See the explanation 54C of the IOPT argument below. 55C If the user wants the iterates printed (NPRINT positive), then 56C FCN must do the printing. See the explanation of NPRINT 57C below. FCN must be declared in an EXTERNAL statement in the 58C calling program and should be written as follows. 59C 60C 61C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC) 62C INTEGER IFLAG,LDFJAC,M,N 63C REAL X(N),FVEC(M) 64C ---------- 65C FJAC and LDFJAC may be ignored , if IOPT=1. 66C REAL FJAC(LDFJAC,N) , if IOPT=2. 67C REAL FJAC(N) , if IOPT=3. 68C ---------- 69C If IFLAG=0, the values in X and FVEC are available 70C for printing. See the explanation of NPRINT below. 71C IFLAG will never be zero unless NPRINT is positive. 72C The values of X and FVEC must not be changed. 73C RETURN 74C ---------- 75C If IFLAG=1, calculate the functions at X and return 76C this vector in FVEC. 77C RETURN 78C ---------- 79C If IFLAG=2, calculate the full Jacobian at X and return 80C this matrix in FJAC. Note that IFLAG will never be 2 unless 81C IOPT=2. FVEC contains the function values at X and must 82C not be altered. FJAC(I,J) must be set to the derivative 83C of FVEC(I) with respect to X(J). 84C RETURN 85C ---------- 86C If IFLAG=3, calculate the LDFJAC-th row of the Jacobian 87C and return this vector in FJAC. Note that IFLAG will 88C never be 3 unless IOPT=3. FVEC contains the function 89C values at X and must not be altered. FJAC(J) must be 90C set to the derivative of FVEC(LDFJAC) with respect to X(J). 91C RETURN 92C ---------- 93C END 94C 95C 96C The value of IFLAG should not be changed by FCN unless the 97C user wants to terminate execution of SNLS1. In this case, set 98C IFLAG to a negative integer. 99C 100C 101C IOPT is an input variable which specifies how the Jacobian will 102C be calculated. If IOPT=2 or 3, then the user must supply the 103C Jacobian, as well as the function values, through the 104C subroutine FCN. If IOPT=2, the user supplies the full 105C Jacobian with one call to FCN. If IOPT=3, the user supplies 106C one row of the Jacobian with each call. (In this manner, 107C storage can be saved because the full Jacobian is not stored.) 108C If IOPT=1, the code will approximate the Jacobian by forward 109C differencing. 110C 111C M is a positive integer input variable set to the number of 112C functions. 113C 114C N is a positive integer input variable set to the number of 115C variables. N must not exceed M. 116C 117C X is an array of length N. On input, X must contain an initial 118C estimate of the solution vector. On output, X contains the 119C final estimate of the solution vector. 120C 121C FVEC is an output array of length M which contains the functions 122C evaluated at the output X. 123C 124C FJAC is an output array. For IOPT=1 and 2, FJAC is an M by N 125C array. For IOPT=3, FJAC is an N by N array. The upper N by N 126C submatrix of FJAC contains an upper triangular matrix R with 127C diagonal elements of nonincreasing magnitude such that 128C 129C T T T 130C P *(JAC *JAC)*P = R *R, 131C 132C where P is a permutation matrix and JAC is the final calcu- 133C lated Jacobian. Column J of P is column IPVT(J) (see below) 134C of the identity matrix. The lower part of FJAC contains 135C information generated during the computation of R. 136C 137C LDFJAC is a positive integer input variable which specifies 138C the leading dimension of the array FJAC. For IOPT=1 and 2, 139C LDFJAC must not be less than M. For IOPT=3, LDFJAC must not 140C be less than N. 141C 142C FTOL is a non-negative input variable. Termination occurs when 143C both the actual and predicted relative reductions in the sum 144C of squares are at most FTOL. Therefore, FTOL measures the 145C relative error desired in the sum of squares. Section 4 con- 146C tains more details about FTOL. 147C 148C XTOL is a non-negative input variable. Termination occurs when 149C the relative error between two consecutive iterates is at most 150C XTOL. Therefore, XTOL measures the relative error desired in 151C the approximate solution. Section 4 contains more details 152C about XTOL. 153C 154C GTOL is a non-negative input variable. Termination occurs when 155C the cosine of the angle between FVEC and any column of the 156C Jacobian is at most GTOL in absolute value. Therefore, GTOL 157C measures the orthogonality desired between the function vector 158C and the columns of the Jacobian. Section 4 contains more 159C details about GTOL. 160C 161C MAXFEV is a positive integer input variable. Termination occurs 162C when the number of calls to FCN to evaluate the functions 163C has reached MAXFEV. 164C 165C EPSFCN is an input variable used in determining a suitable step 166C for the forward-difference approximation. This approximation 167C assumes that the relative errors in the functions are of the 168C order of EPSFCN. If EPSFCN is less than the machine preci- 169C sion, it is assumed that the relative errors in the functions 170C are of the order of the machine precision. If IOPT=2 or 3, 171C then EPSFCN can be ignored (treat it as a dummy argument). 172C 173C DIAG is an array of length N. If MODE = 1 (see below), DIAG is 174C internally set. If MODE = 2, DIAG must contain positive 175C entries that serve as implicit (multiplicative) scale factors 176C for the variables. 177C 178C MODE is an integer input variable. If MODE = 1, the variables 179C will be scaled internally. If MODE = 2, the scaling is speci- 180C fied by the input DIAG. Other values of MODE are equivalent 181C to MODE = 1. 182C 183C FACTOR is a positive input variable used in determining the ini- 184C tial step bound. This bound is set to the product of FACTOR 185C and the Euclidean norm of DIAG*X if nonzero, or else to FACTOR 186C itself. In most cases FACTOR should lie in the interval 187C (.1,100.). 100. is a generally recommended value. 188C 189C NPRINT is an integer input variable that enables controlled 190C printing of iterates if it is positive. In this case, FCN is 191C called with IFLAG = 0 at the beginning of the first iteration 192C and every NPRINT iterations thereafter and immediately prior 193C to return, with X and FVEC available for printing. Appropriate 194C print statements must be added to FCN (see example) and 195C FVEC should not be altered. If NPRINT is not positive, no 196C special calls to FCN with IFLAG = 0 are made. 197C 198C INFO is an integer output variable. If the user has terminated 199C execution, INFO is set to the (negative) value of IFLAG. See 200C description of FCN and JAC. Otherwise, INFO is set as follows. 201C 202C INFO = 0 improper input parameters. 203C 204C INFO = 1 both actual and predicted relative reductions in the 205C sum of squares are at most FTOL. 206C 207C INFO = 2 relative error between two consecutive iterates is 208C at most XTOL. 209C 210C INFO = 3 conditions for INFO = 1 and INFO = 2 both hold. 211C 212C INFO = 4 the cosine of the angle between FVEC and any column 213C of the Jacobian is at most GTOL in absolute value. 214C 215C INFO = 5 number of calls to FCN for function evaluation 216C has reached MAXFEV. 217C 218C INFO = 6 FTOL is too small. No further reduction in the sum 219C of squares is possible. 220C 221C INFO = 7 XTOL is too small. No further improvement in the 222C approximate solution X is possible. 223C 224C INFO = 8 GTOL is too small. FVEC is orthogonal to the 225C columns of the Jacobian to machine precision. 226C 227C Sections 4 and 5 contain more details about INFO. 228C 229C NFEV is an integer output variable set to the number of calls to 230C FCN for function evaluation. 231C 232C NJEV is an integer output variable set to the number of 233C evaluations of the full Jacobian. If IOPT=2, only one call to 234C FCN is required for each evaluation of the full Jacobian. 235C If IOPT=3, the M calls to FCN are required. 236C If IOPT=1, then NJEV is set to zero. 237C 238C IPVT is an integer output array of length N. IPVT defines a 239C permutation matrix P such that JAC*P = Q*R, where JAC is the 240C final calculated Jacobian, Q is orthogonal (not stored), and R 241C is upper triangular with diagonal elements of nonincreasing 242C magnitude. Column J of P is column IPVT(J) of the identity 243C matrix. 244C 245C QTF is an output array of length N which contains the first N 246C elements of the vector (Q transpose)*FVEC. 247C 248C WA1, WA2, and WA3 are work arrays of length N. 249C 250C WA4 is a work array of length M. 251C 252C 253C 4. Successful Completion. 254C 255C The accuracy of SNLS1 is controlled by the convergence parame- 256C ters FTOL, XTOL, and GTOL. These parameters are used in tests 257C which make three types of comparisons between the approximation 258C X and a solution XSOL. SNLS1 terminates when any of the tests 259C is satisfied. If any of the convergence parameters is less than 260C the machine precision (as defined by the function R1MACH(4)), 261C then SNLS1 only attempts to satisfy the test defined by the 262C machine precision. Further progress is not usually possible. 263C 264C The tests assume that the functions are reasonably well behaved, 265C and, if the Jacobian is supplied by the user, that the functions 266C and the Jacobian are coded consistently. If these conditions 267C are not satisfied, then SNLS1 may incorrectly indicate conver- 268C gence. If the Jacobian is coded correctly or IOPT=1, 269C then the validity of the answer can be checked, for example, by 270C rerunning SNLS1 with tighter tolerances. 271C 272C First Convergence Test. If ENORM(Z) denotes the Euclidean norm 273C of a vector Z, then this test attempts to guarantee that 274C 275C ENORM(FVEC) .LE. (1+FTOL)*ENORM(FVECS), 276C 277C where FVECS denotes the functions evaluated at XSOL. If this 278C condition is satisfied with FTOL = 10**(-K), then the final 279C residual norm ENORM(FVEC) has K significant decimal digits and 280C INFO is set to 1 (or to 3 if the second test is also satis- 281C fied). Unless high precision solutions are required, the 282C recommended value for FTOL is the square root of the machine 283C precision. 284C 285C Second Convergence Test. If D is the diagonal matrix whose 286C entries are defined by the array DIAG, then this test attempts 287C to guarantee that 288C 289C ENORM(D*(X-XSOL)) .LE. XTOL*ENORM(D*XSOL). 290C 291C If this condition is satisfied with XTOL = 10**(-K), then the 292C larger components of D*X have K significant decimal digits and 293C INFO is set to 2 (or to 3 if the first test is also satis- 294C fied). There is a danger that the smaller components of D*X 295C may have large relative errors, but if MODE = 1, then the 296C accuracy of the components of X is usually related to their 297C sensitivity. Unless high precision solutions are required, 298C the recommended value for XTOL is the square root of the 299C machine precision. 300C 301C Third Convergence Test. This test is satisfied when the cosine 302C of the angle between FVEC and any column of the Jacobian at X 303C is at most GTOL in absolute value. There is no clear rela- 304C tionship between this test and the accuracy of SNLS1, and 305C furthermore, the test is equally well satisfied at other crit- 306C ical points, namely maximizers and saddle points. Therefore, 307C termination caused by this test (INFO = 4) should be examined 308C carefully. The recommended value for GTOL is zero. 309C 310C 311C 5. Unsuccessful Completion. 312C 313C Unsuccessful termination of SNLS1 can be due to improper input 314C parameters, arithmetic interrupts, or an excessive number of 315C function evaluations. 316C 317C Improper Input Parameters. INFO is set to 0 if IOPT .LT. 1 318C or IOPT .GT. 3, or N .LE. 0, or M .LT. N, or for IOPT=1 or 2 319C LDFJAC .LT. M, or for IOPT=3 LDFJAC .LT. N, or FTOL .LT. 0.E0, 320C or XTOL .LT. 0.E0, or GTOL .LT. 0.E0, or MAXFEV .LE. 0, or 321C FACTOR .LE. 0.E0. 322C 323C Arithmetic Interrupts. If these interrupts occur in the FCN 324C subroutine during an early stage of the computation, they may 325C be caused by an unacceptable choice of X by SNLS1. In this 326C case, it may be possible to remedy the situation by rerunning 327C SNLS1 with a smaller value of FACTOR. 328C 329C Excessive Number of Function Evaluations. A reasonable value 330C for MAXFEV is 100*(N+1) for IOPT=2 or 3 and 200*(N+1) for 331C IOPT=1. If the number of calls to FCN reaches MAXFEV, then 332C this indicates that the routine is converging very slowly 333C as measured by the progress of FVEC, and INFO is set to 5. 334C In this case, it may be helpful to restart SNLS1 with MODE 335C set to 1. 336C 337C 338C 6. Characteristics of the Algorithm. 339C 340C SNLS1 is a modification of the Levenberg-Marquardt algorithm. 341C Two of its main characteristics involve the proper use of 342C implicitly scaled variables (if MODE = 1) and an optimal choice 343C for the correction. The use of implicitly scaled variables 344C achieves scale invariance of SNLS1 and limits the size of the 345C correction in any direction where the functions are changing 346C rapidly. The optimal choice of the correction guarantees (under 347C reasonable conditions) global convergence from starting points 348C far from the solution and a fast rate of convergence for 349C problems with small residuals. 350C 351C Timing. The time required by SNLS1 to solve a given problem 352C depends on M and N, the behavior of the functions, the accu- 353C racy requested, and the starting point. The number of arith- 354C metic operations needed by SNLS1 is about N**3 to process each 355C evaluation of the functions (call to FCN) and to process each 356C evaluation of the Jacobian it takes M*N**2 for IOPT=2 (one 357C call to FCN), M*N**2 for IOPT=1 (N calls to FCN) and 358C 1.5*M*N**2 for IOPT=3 (M calls to FCN). Unless FCN 359C can be evaluated quickly, the timing of SNLS1 will be 360C strongly influenced by the time spent in FCN. 361C 362C Storage. SNLS1 requires (M*N + 2*M + 6*N) for IOPT=1 or 2 and 363C (N**2 + 2*M + 6*N) for IOPT=3 single precision storage 364C locations and N integer storage locations, in addition to 365C the storage required by the program. There are no internally 366C declared storage arrays. 367C 368C ***See the LONG DESCRIPTION for an example problem*** 369C***LONG DESCRIPTION 370C 371C 7. Example. 372C 373C The problem is to determine the values of X(1), X(2), and X(3) 374C which provide the best fit (in the least squares sense) of 375C 376C X(1) + U(I)/(V(I)*X(2) + W(I)*X(3)), I = 1, 15 377C 378C to the data 379C 380C Y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39, 381C 0.37,0.58,0.73,0.96,1.34,2.10,4.39), 382C 383C where U(I) = I, V(I) = 16 - I, and W(I) = MIN(U(I),V(I)). The 384C I-th component of FVEC is thus defined by 385C 386C Y(I) - (X(1) + U(I)/(V(I)*X(2) + W(I)*X(3))). 387C 388C ********** 389C 390C PROGRAM TEST(INPUT,OUTPUT,TAPE6=OUTPUT) 391C C 392C C Driver for SNLS1 example. 393C C 394C INTEGER J,IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV, 395C * NWRITE 396C INTEGER IPVT(3) 397C REAL FTOL,XTOL,GTOL,FACTOR,FNORM,EPSFCN 398C REAL X(3),FVEC(15),FJAC(15,3),DIAG(3),QTF(3), 399C * WA1(3),WA2(3),WA3(3),WA4(15) 400C REAL ENORM,R1MACH 401C EXTERNAL FCN 402C DATA NWRITE /6/ 403C C 404C IOPT = 1 405C M = 15 406C N = 3 407C C 408C C The following starting values provide a rough fit. 409C C 410C X(1) = 1.E0 411C X(2) = 1.E0 412C X(3) = 1.E0 413C C 414C LDFJAC = 15 415C C 416C C Set FTOL and XTOL to the square root of the machine precision 417C C and GTOL to zero. Unless high precision solutions are 418C C required, these are the recommended settings. 419C C 420C FTOL = SQRT(R1MACH(4)) 421C XTOL = SQRT(R1MACH(4)) 422C GTOL = 0.E0 423C C 424C MAXFEV = 400 425C EPSFCN = 0.0 426C MODE = 1 427C FACTOR = 1.E2 428C NPRINT = 0 429C C 430C CALL SNLS1(FCN,IOPT,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL, 431C * GTOL,MAXFEV,EPSFCN,DIAG,MODE,FACTOR,NPRINT, 432C * INFO,NFEV,NJEV,IPVT,QTF,WA1,WA2,WA3,WA4) 433C FNORM = ENORM(M,FVEC) 434C WRITE (NWRITE,1000) FNORM,NFEV,NJEV,INFO,(X(J),J=1,N) 435C STOP 436C 1000 FORMAT (5X,' FINAL L2 NORM OF THE RESIDUALS',E15.7 // 437C * 5X,' NUMBER OF FUNCTION EVALUATIONS',I10 // 438C * 5X,' NUMBER OF JACOBIAN EVALUATIONS',I10 // 439C * 5X,' EXIT PARAMETER',16X,I10 // 440C * 5X,' FINAL APPROXIMATE SOLUTION' // 5X,3E15.7) 441C END 442C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,DUM,IDUM) 443C C This is the form of the FCN routine if IOPT=1, 444C C that is, if the user does not calculate the Jacobian. 445C INTEGER M,N,IFLAG 446C REAL X(N),FVEC(M) 447C INTEGER I 448C REAL TMP1,TMP2,TMP3,TMP4 449C REAL Y(15) 450C DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8), 451C * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15) 452C * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1, 453C * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/ 454C C 455C IF (IFLAG .NE. 0) GO TO 5 456C C 457C C Insert print statements here when NPRINT is positive. 458C C 459C RETURN 460C 5 CONTINUE 461C DO 10 I = 1, M 462C TMP1 = I 463C TMP2 = 16 - I 464C TMP3 = TMP1 465C IF (I .GT. 8) TMP3 = TMP2 466C FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3)) 467C 10 CONTINUE 468C RETURN 469C END 470C 471C 472C Results obtained with different compilers or machines 473C may be slightly different. 474C 475C FINAL L2 NORM OF THE RESIDUALS 0.9063596E-01 476C 477C NUMBER OF FUNCTION EVALUATIONS 25 478C 479C NUMBER OF JACOBIAN EVALUATIONS 0 480C 481C EXIT PARAMETER 1 482C 483C FINAL APPROXIMATE SOLUTION 484C 485C 0.8241058E-01 0.1133037E+01 0.2343695E+01 486C 487C 488C For IOPT=2, FCN would be modified as follows to also 489C calculate the full Jacobian when IFLAG=2. 490C 491C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC) 492C C 493C C This is the form of the FCN routine if IOPT=2, 494C C that is, if the user calculates the full Jacobian. 495C C 496C INTEGER LDFJAC,M,N,IFLAG 497C REAL X(N),FVEC(M) 498C REAL FJAC(LDFJAC,N) 499C INTEGER I 500C REAL TMP1,TMP2,TMP3,TMP4 501C REAL Y(15) 502C DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8), 503C * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15) 504C * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1, 505C * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/ 506C C 507C IF (IFLAG .NE. 0) GO TO 5 508C C 509C C Insert print statements here when NPRINT is positive. 510C C 511C RETURN 512C 5 CONTINUE 513C IF(IFLAG.NE.1) GO TO 20 514C DO 10 I = 1, M 515C TMP1 = I 516C TMP2 = 16 - I 517C TMP3 = TMP1 518C IF (I .GT. 8) TMP3 = TMP2 519C FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3)) 520C 10 CONTINUE 521C RETURN 522C C 523C C Below, calculate the full Jacobian. 524C C 525C 20 CONTINUE 526C C 527C DO 30 I = 1, M 528C TMP1 = I 529C TMP2 = 16 - I 530C TMP3 = TMP1 531C IF (I .GT. 8) TMP3 = TMP2 532C TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2 533C FJAC(I,1) = -1.E0 534C FJAC(I,2) = TMP1*TMP2/TMP4 535C FJAC(I,3) = TMP1*TMP3/TMP4 536C 30 CONTINUE 537C RETURN 538C END 539C 540C 541C For IOPT = 3, FJAC would be dimensioned as FJAC(3,3), 542C LDFJAC would be set to 3, and FCN would be written as 543C follows to calculate a row of the Jacobian when IFLAG=3. 544C 545C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC) 546C C This is the form of the FCN routine if IOPT=3, 547C C that is, if the user calculates the Jacobian row by row. 548C INTEGER M,N,IFLAG 549C REAL X(N),FVEC(M) 550C REAL FJAC(N) 551C INTEGER I 552C REAL TMP1,TMP2,TMP3,TMP4 553C REAL Y(15) 554C DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8), 555C * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15) 556C * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1, 557C * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/ 558C C 559C IF (IFLAG .NE. 0) GO TO 5 560C C 561C C Insert print statements here when NPRINT is positive. 562C C 563C RETURN 564C 5 CONTINUE 565C IF( IFLAG.NE.1) GO TO 20 566C DO 10 I = 1, M 567C TMP1 = I 568C TMP2 = 16 - I 569C TMP3 = TMP1 570C IF (I .GT. 8) TMP3 = TMP2 571C FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3)) 572C 10 CONTINUE 573C RETURN 574C C 575C C Below, calculate the LDFJAC-th row of the Jacobian. 576C C 577C 20 CONTINUE 578C 579C I = LDFJAC 580C TMP1 = I 581C TMP2 = 16 - I 582C TMP3 = TMP1 583C IF (I .GT. 8) TMP3 = TMP2 584C TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2 585C FJAC(1) = -1.E0 586C FJAC(2) = TMP1*TMP2/TMP4 587C FJAC(3) = TMP1*TMP3/TMP4 588C RETURN 589C END 590C***REFERENCES MORE, JORGE J. 591C THE LEVENBERG-MARQUARDT ALGORITHM, IMPLEMENTATION AND 592C THEORY. NUMERICAL ANALYSIS, G. A. WATSON, EDITOR. 593C LECTURE NOTES IN MATHEMATICS 630, SPRINGER-VERLAG, 594C 1977. 595C***ROUTINES CALLED CHKDER,ENORM,FDJAC3,LMPAR,QRFAC,R1MACH,RWUPDT, 596C XERROR,XERRWV 597C***END PROLOGUE SNLS1 598 599 600