1 0 Page 0 Documentation for MINPACK subroutine HYBRD1 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of HYBRD1 is to find a zero of a system of N non- linear functions in N variables by a modification of the Powell hybrid method. This is done by using the more general nonlinea equation solver HYBRD. The user must provide a subroutine whic calculates the functions. The Jacobian is then calculated by a forward-difference approximation. 0 2. Subroutine and type statements. 0 SUBROUTINE HYBRD1(FCN,N,X,FVEC,TOL,INFO,WA,LWA) INTEGER N,INFO,LWA REAL TOL REAL X(N),FVEC(N),WA(LWA) EXTERNAL FCN 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to HYBRD1 and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from HYBRD1. 0 FCN is the name of the user-supplied subroutine which calculate the functions. FCN must be declared in an EXTERNAL statement in the user calling program, and should be written as follows 0 SUBROUTINE FCN(N,X,FVEC,IFLAG) INTEGER N,IFLAG REAL X(N),FVEC(N) ---------- CALCULATE THE FUNCTIONS AT X AND RETURN THIS VECTOR IN FVEC. ---------- RETURN END 0 The value of IFLAG should not be changed by FCN unless the user wants to terminate execution of HYBRD1. In this case se IFLAG to a negative integer. 1 0 Page 0 N is a positive integer input variable set to the number of functions and variables. 0 X is an array of length N. On input X must contain an initial estimate of the solution vector. On output X contains the final estimate of the solution vector. 0 FVEC is an output array of length N which contains the function evaluated at the output X. 0 TOL is a nonnegative input variable. Termination occurs when the algorithm estimates that the relative error between X and the solution is at most TOL. Section 4 contains more details about TOL. 0 INFO is an integer output variable. If the user has terminated execution, INFO is set to the (negative) value of IFLAG. See description of FCN. Otherwise, INFO is set as follows. 0 INFO = 0 Improper input parameters. 0 INFO = 1 Algorithm estimates that the relative error between X and the solution is at most TOL. 0 INFO = 2 Number of calls to FCN has reached or exceeded 200*(N+1). 0 INFO = 3 TOL is too small. No further improvement in the approximate solution X is possible. 0 INFO = 4 Iteration is not making good progress. 0 Sections 4 and 5 contain more details about INFO. 0 WA is a work array of length LWA. 0 LWA is a positive integer input variable not less than (N*(3*N+13))/2. 0 4. Successful completion. 0 The accuracy of HYBRD1 is controlled by the convergence parame- ter TOL. This parameter is used in a test which makes a compar ison between the approximation X and a solution XSOL. HYBRD1 terminates when the test is satisfied. If TOL is less than the machine precision (as defined by the MINPACK function SPMPAR(1)), then HYBRD1 only attempts to satisfy the test defined by the machine precision. Further progress is not usu- ally possible. Unless high precision solutions are required, the recommended value for TOL is the square root of the machine precision. 0 The test assumes that the functions are reasonably well behaved 1 0 Page 0 If this condition is not satisfied, then HYBRD1 may incorrectly indicate convergence. The validity of the answer can be checked, for example, by rerunning HYBRD1 with a tighter toler- ance. 0 Convergence test. If ENORM(Z) denotes the Euclidean norm of a vector Z, then this test attempts to guarantee that 0 ENORM(X-XSOL) .LE. TOL*ENORM(XSOL). 0 If this condition is satisfied with TOL = 10**(-K), then the larger components of X have K significant decimal digits and INFO is set to 1. There is a danger that the smaller compo- nents of X may have large relative errors, but the fast rate of convergence of HYBRD1 usually avoids this possibility. 0 5. Unsuccessful completion. 0 Unsuccessful termination of HYBRD1 can be due to improper input parameters, arithmetic interrupts, an excessive number of func- tion evaluations, errors in the functions, or lack of good prog ress. 0 Improper input parameters. INFO is set to 0 if N .LE. 0, or TOL .LT. 0.E0, or LWA .LT. (N*(3*N+13))/2. 0 Arithmetic interrupts. If these interrupts occur in the FCN subroutine during an early stage of the computation, they may be caused by an unacceptable choice of X by HYBRD1. In this case, it may be possible to remedy the situation by not evalu ating the functions here, but instead setting the components of FVEC to numbers that exceed those in the initial FVEC, thereby indirectly reducing the step length. The step length can be more directly controlled by using instead HYBRD, which includes in its calling sequence the step-length- governing parameter FACTOR. 0 Excessive number of function evaluations. If the number of calls to FCN reaches 200*(N+1), then this indicates that the routine is converging very slowly as measured by the progress of FVEC, and INFO is set to 2. This situation should be unu- sual because, as indicated below, lack of good progress is usually diagnosed earlier by HYBRD1, causing termination with INFO = 4. 0 Errors in the functions. The choice of step length in the for- ward-difference approximation to the Jacobian assumes that th relative errors in the functions are of the order of the machine precision. If this is not the case, HYBRD1 may fail (usually with INFO = 4). The user should then use HYBRD instead, or one of the programs which require the analytic Jacobian (HYBRJ1 and HYBRJ). 1 0 Page 0 Lack of good progress. HYBRD1 searches for a zero of the syste by minimizing the sum of the squares of the functions. In so doing, it can become trapped in a region where the minimum does not correspond to a zero of the system and, in this situ ation, the iteration eventually fails to make good progress. In particular, this will happen if the system does not have a zero. If the system has a zero, rerunning HYBRD1 from a dif- ferent starting point may be helpful. 0 6. Characteristics of the algorithm. 0 HYBRD1 is a modification of the Powell hybrid method. Two of its main characteristics involve the choice of the correction a a convex combination of the Newton and scaled gradient direc- tions, and the updating of the Jacobian by the rank-1 method of Broyden. The choice of the correction guarantees (under reason able conditions) global convergence for starting points far fro the solution and a fast rate of convergence. The Jacobian is approximated by forward differences at the starting point, but forward differences are not used again until the rank-1 method fails to produce satisfactory progress. 0 Timing. The time required by HYBRD1 to solve a given problem depends on N, the behavior of the functions, the accuracy requested, and the starting point. The number of arithmetic operations needed by HYBRD1 is about 11.5*(N**2) to process each call to FCN. Unless FCN can be evaluated quickly, the timing of HYBRD1 will be strongly influenced by the time spen in FCN. 0 Storage. HYBRD1 requires (3*N**2 + 17*N)/2 single precision storage locations, in addition to the storage required by the program. There are no internally declared storage arrays. 0 7. Subprograms required. 0 USER-supplied ...... FCN 0 MINPACK-supplied ... DOGLEG,SPMPAR,ENORM,FDJAC1,HYBRD, QFORM,QRFAC,R1MPYQ,R1UPDT 0 FORTRAN-supplied ... ABS,AMAX1,AMIN1,SQRT,MIN0,MOD 0 8. References. 0 M. J. D. Powell, A Hybrid Method for Nonlinear Equations. Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, editor. Gordon and Breach, 1970. 0 9. Example. 1 0 Page 0 The problem is to determine the values of x(1), x(2), ..., x(9) which solve the system of tridiagonal equations 0 (3-2*x(1))*x(1) -2*x(2) = -1 -x(i-1) + (3-2*x(i))*x(i) -2*x(i+1) = -1, i=2-8 -x(8) + (3-2*x(9))*x(9) = -1 0 C ********** C C DRIVER FOR HYBRD1 EXAMPLE. C SINGLE PRECISION VERSION C C ********** INTEGER J,N,INFO,LWA,NWRITE REAL TOL,FNORM REAL X(9),FVEC(9),WA(180) REAL ENORM,SPMPAR EXTERNAL FCN C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C N = 9 C C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH SOLUTION. C DO 10 J = 1, 9 X(J) = -1.E0 10 CONTINUE C LWA = 180 C C SET TOL TO THE SQUARE ROOT OF THE MACHINE PRECISION. C UNLESS HIGH PRECISION SOLUTIONS ARE REQUIRED, C THIS IS THE RECOMMENDED SETTING. C TOL = SQRT(SPMPAR(1)) C CALL HYBRD1(FCN,N,X,FVEC,TOL,INFO,WA,LWA) FNORM = ENORM(N,FVEC) WRITE (NWRITE,1000) FNORM,INFO,(X(J),J=1,N) STOP 1000 FORMAT (5X,31H FINAL L2 NORM OF THE RESIDUALS,E15.7 // * 5X,15H EXIT PARAMETER,16X,I10 // * 5X,27H FINAL APPROXIMATE SOLUTION // (5X,3E15.7)) C C LAST CARD OF DRIVER FOR HYBRD1 EXAMPLE. C END SUBROUTINE FCN(N,X,FVEC,IFLAG) INTEGER N,IFLAG REAL X(N),FVEC(N) C 1 0 Page 0 C SUBROUTINE FCN FOR HYBRD1 EXAMPLE. C INTEGER K REAL ONE,TEMP,TEMP1,TEMP2,THREE,TWO,ZERO DATA ZERO,ONE,TWO,THREE /0.E0,1.E0,2.E0,3.E0/ C DO 10 K = 1, N TEMP = (THREE - TWO*X(K))*X(K) TEMP1 = ZERO IF (K .NE. 1) TEMP1 = X(K-1) TEMP2 = ZERO IF (K .NE. N) TEMP2 = X(K+1) FVEC(K) = TEMP - TEMP1 - TWO*TEMP2 + ONE 10 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 0 Results obtained with different compilers or machines may be slightly different. 0 FINAL L2 NORM OF THE RESIDUALS 0.1192636E-07 0 EXIT PARAMETER 1 0 FINAL APPROXIMATE SOLUTION 0 -0.5706545E+00 -0.6816283E+00 -0.7017325E+00 -0.7042129E+00 -0.7013690E+00 -0.6918656E+00 -0.6657920E+00 -0.5960342E+00 -0.4164121E+00 1 0 1 0 Page 0 Documentation for MINPACK subroutine HYBRD 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of HYBRD is to find a zero of a system of N non- linear functions in N variables by a modification of the Powell hybrid method. The user must provide a subroutine which calcu- lates the functions. The Jacobian is then calculated by a for- ward-difference approximation. 0 2. Subroutine and type statements. 0 SUBROUTINE HYBRD(FCN,N,X,FVEC,XTOL,MAXFEV,ML,MU,EPSFCN,DIAG, * MODE,FACTOR,NPRINT,INFO,NFEV,FJAC,LDFJAC, * R,LR,QTF,WA1,WA2,WA3,WA4) INTEGER N,MAXFEV,ML,MU,MODE,NPRINT,INFO,NFEV,LDFJAC,LR REAL XTOL,EPSFCN,FACTOR REAL X(N),FVEC(N),DIAG(N),FJAC(LDFJAC,N),R(LR),QTF(N), * WA1(N),WA2(N),WA3(N),WA4(N) EXTERNAL FCN 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to HYBRD and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from HYBRD. 0 FCN is the name of the user-supplied subroutine which calculate the functions. FCN must be declared in an EXTERNAL statement in the user calling program, and should be written as follows 0 SUBROUTINE FCN(N,X,FVEC,IFLAG) INTEGER N,IFLAG REAL X(N),FVEC(N) ---------- CALCULATE THE FUNCTIONS AT X AND RETURN THIS VECTOR IN FVEC. ---------- RETURN END 0 The value of IFLAG should not be changed by FCN unless the 1 0 Page 0 user wants to terminate execution of HYBRD. In this case set IFLAG to a negative integer. 0 N is a positive integer input variable set to the number of functions and variables. 0 X is an array of length N. On input X must contain an initial estimate of the solution vector. On output X contains the final estimate of the solution vector. 0 FVEC is an output array of length N which contains the function evaluated at the output X. 0 XTOL is a nonnegative input variable. Termination occurs when the relative error between two consecutive iterates is at mos XTOL. Therefore, XTOL measures the relative error desired in the approximate solution. Section 4 contains more details about XTOL. 0 MAXFEV is a positive integer input variable. Termination occur when the number of calls to FCN is at least MAXFEV by the end of an iteration. 0 ML is a nonnegative integer input variable which specifies the number of subdiagonals within the band of the Jacobian matrix If the Jacobian is not banded, set ML to at least N - 1. 0 MU is a nonnegative integer input variable which specifies the number of superdiagonals within the band of the Jacobian matrix. If the Jacobian is not banded, set MU to at least N - 1. 0 EPSFCN is an input variable used in determining a suitable step for the forward-difference approximation. This approximation assumes that the relative errors in the functions are of the order of EPSFCN. If EPSFCN is less than the machine preci- sion, it is assumed that the relative errors in the functions are of the order of the machine precision. 0 DIAG is an array of length N. If MODE = 1 (see below), DIAG is internally set. If MODE = 2, DIAG must contain positive entries that serve as multiplicative scale factors for the variables. 0 MODE is an integer input variable. If MODE = 1, the variables will be scaled internally. If MODE = 2, the scaling is speci fied by the input DIAG. Other values of MODE are equivalent to MODE = 1. 0 FACTOR is a positive input variable used in determining the ini tial step bound. This bound is set to the product of FACTOR and the Euclidean norm of DIAG*X if nonzero, or else to FACTO itself. In most cases FACTOR should lie in the interval (.1,100.). 100. is a generally recommended value. 1 0 Page 0 NPRINT is an integer input variable that enables controlled printing of iterates if it is positive. In this case, FCN is called with IFLAG = 0 at the beginning of the first iteration and every NPRINT iterations thereafter and immediately prior to return, with X and FVEC available for printing. If NPRINT is not positive, no special calls of FCN with IFLAG = 0 are made. 0 INFO is an integer output variable. If the user has terminated execution, INFO is set to the (negative) value of IFLAG. See description of FCN. Otherwise, INFO is set as follows. 0 INFO = 0 Improper input parameters. 0 INFO = 1 Relative error between two consecutive iterates is at most XTOL. 0 INFO = 2 Number of calls to FCN has reached or exceeded MAXFEV. 0 INFO = 3 XTOL is too small. No further improvement in the approximate solution X is possible. 0 INFO = 4 Iteration is not making good progress, as measured by the improvement from the last five Jacobian eval uations. 0 INFO = 5 Iteration is not making good progress, as measured by the improvement from the last ten iterations. 0 Sections 4 and 5 contain more details about INFO. 0 NFEV is an integer output variable set to the number of calls t FCN. 0 FJAC is an output N by N array which contains the orthogonal matrix Q produced by the QR factorization of the final approx imate Jacobian. 0 LDFJAC is a positive integer input variable not less than N which specifies the leading dimension of the array FJAC. 0 R is an output array of length LR which contains the upper triangular matrix produced by the QR factorization of the final approximate Jacobian, stored rowwise. 0 LR is a positive integer input variable not less than (N*(N+1))/2. 0 QTF is an output array of length N which contains the vector (Q transpose)*FVEC. 0 WA1, WA2, WA3, and WA4 are work arrays of length N. 1 0 Page 0 4. Successful completion. 0 The accuracy of HYBRD is controlled by the convergence paramete XTOL. This parameter is used in a test which makes a compariso between the approximation X and a solution XSOL. HYBRD termi- nates when the test is satisfied. If the convergence parameter is less than the machine precision (as defined by the MINPACK function SPMPAR(1)), then HYBRD only attempts to satisfy the test defined by the machine precision. Further progress is not usually possible. 0 The test assumes that the functions are reasonably well behaved If this condition is not satisfied, then HYBRD may incorrectly indicate convergence. The validity of the answer can be checked, for example, by rerunning HYBRD with a tighter toler- ance. 0 Convergence test. If ENORM(Z) denotes the Euclidean norm of a vector Z and D is the diagonal matrix whose entries are defined by the array DIAG, then this test attempts to guaran- tee that 0 ENORM(D*(X-XSOL)) .LE. XTOL*ENORM(D*XSOL). 0 If this condition is satisfied with XTOL = 10**(-K), then the larger components of D*X have K significant decimal digits an INFO is set to 1. There is a danger that the smaller compo- nents of D*X may have large relative errors, but the fast rat of convergence of HYBRD usually avoids this possibility. Unless high precision solutions are required, the recommended value for XTOL is the square root of the machine precision. 0 5. Unsuccessful completion. 0 Unsuccessful termination of HYBRD can be due to improper input parameters, arithmetic interrupts, an excessive number of func- tion evaluations, or lack of good progress. 0 Improper input parameters. INFO is set to 0 if N .LE. 0, or XTOL .LT. 0.E0, or MAXFEV .LE. 0, or ML .LT. 0, or MU .LT. 0, or FACTOR .LE. 0.E0, or LDFJAC .LT. N, or LR .LT. (N*(N+1))/2 0 Arithmetic interrupts. If these interrupts occur in the FCN subroutine during an early stage of the computation, they may be caused by an unacceptable choice of X by HYBRD. In this case, it may be possible to remedy the situation by rerunning HYBRD with a smaller value of FACTOR. 0 Excessive number of function evaluations. A reasonable value for MAXFEV is 200*(N+1). If the number of calls to FCN reaches MAXFEV, then this indicates that the routine is con- verging very slowly as measured by the progress of FVEC, and 1 0 Page 0 INFO is set to 2. This situation should be unusual because, as indicated below, lack of good progress is usually diagnose earlier by HYBRD, causing termination with INFO = 4 or INFO = 5. 0 Lack of good progress. HYBRD searches for a zero of the system by minimizing the sum of the squares of the functions. In so doing, it can become trapped in a region where the minimum does not correspond to a zero of the system and, in this situ ation, the iteration eventually fails to make good progress. In particular, this will happen if the system does not have a zero. If the system has a zero, rerunning HYBRD from a dif- ferent starting point may be helpful. 0 6. Characteristics of the algorithm. 0 HYBRD is a modification of the Powell hybrid method. Two of it main characteristics involve the choice of the correction as a convex combination of the Newton and scaled gradient directions and the updating of the Jacobian by the rank-1 method of Broy- den. The choice of the correction guarantees (under reasonable conditions) global convergence for starting points far from the solution and a fast rate of convergence. The Jacobian is approximated by forward differences at the starting point, but forward differences are not used again until the rank-1 method fails to produce satisfactory progress. 0 Timing. The time required by HYBRD to solve a given problem depends on N, the behavior of the functions, the accuracy requested, and the starting point. The number of arithmetic operations needed by HYBRD is about 11.5*(N**2) to process each call to FCN. Unless FCN can be evaluated quickly, the timing of HYBRD will be strongly influenced by the time spent in FCN. 0 Storage. HYBRD requires (3*N**2 + 17*N)/2 single precision storage locations, in addition to the storage required by the program. There are no internally declared storage arrays. 0 7. Subprograms required. 0 USER-supplied ...... FCN 0 MINPACK-supplied ... DOGLEG,SPMPAR,ENORM,FDJAC1, QFORM,QRFAC,R1MPYQ,R1UPDT 0 FORTRAN-supplied ... ABS,AMAX1,AMIN1,SQRT,MIN0,MOD 0 8. References. 0 M. J. D. Powell, A Hybrid Method for Nonlinear Equations. 1 0 Page 0 Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, editor. Gordon and Breach, 1970. 0 9. Example. 0 The problem is to determine the values of x(1), x(2), ..., x(9) which solve the system of tridiagonal equations 0 (3-2*x(1))*x(1) -2*x(2) = -1 -x(i-1) + (3-2*x(i))*x(i) -2*x(i+1) = -1, i=2-8 -x(8) + (3-2*x(9))*x(9) = -1 0 C ********** C C DRIVER FOR HYBRD EXAMPLE. C SINGLE PRECISION VERSION C C ********** INTEGER J,N,MAXFEV,ML,MU,MODE,NPRINT,INFO,NFEV,LDFJAC,LR,NWRITE REAL XTOL,EPSFCN,FACTOR,FNORM REAL X(9),FVEC(9),DIAG(9),FJAC(9,9),R(45),QTF(9), * WA1(9),WA2(9),WA3(9),WA4(9) REAL ENORM,SPMPAR EXTERNAL FCN C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C N = 9 C C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH SOLUTION. C DO 10 J = 1, 9 X(J) = -1.E0 10 CONTINUE C LDFJAC = 9 LR = 45 C C SET XTOL TO THE SQUARE ROOT OF THE MACHINE PRECISION. C UNLESS HIGH PRECISION SOLUTIONS ARE REQUIRED, C THIS IS THE RECOMMENDED SETTING. C XTOL = SQRT(SPMPAR(1)) C MAXFEV = 2000 ML = 1 MU = 1 EPSFCN = 0.E0 MODE = 2 DO 20 J = 1, 9 DIAG(J) = 1.E0 1 0 Page 0 20 CONTINUE FACTOR = 1.E2 NPRINT = 0 C CALL HYBRD(FCN,N,X,FVEC,XTOL,MAXFEV,ML,MU,EPSFCN,DIAG, * MODE,FACTOR,NPRINT,INFO,NFEV,FJAC,LDFJAC, * R,LR,QTF,WA1,WA2,WA3,WA4) FNORM = ENORM(N,FVEC) WRITE (NWRITE,1000) FNORM,NFEV,INFO,(X(J),J=1,N) STOP 1000 FORMAT (5X,31H FINAL L2 NORM OF THE RESIDUALS,E15.7 // * 5X,31H NUMBER OF FUNCTION EVALUATIONS,I10 // * 5X,15H EXIT PARAMETER,16X,I10 // * 5X,27H FINAL APPROXIMATE SOLUTION // (5X,3E15.7)) C C LAST CARD OF DRIVER FOR HYBRD EXAMPLE. C END SUBROUTINE FCN(N,X,FVEC,IFLAG) INTEGER N,IFLAG REAL X(N),FVEC(N) C C SUBROUTINE FCN FOR HYBRD EXAMPLE. C INTEGER K REAL ONE,TEMP,TEMP1,TEMP2,THREE,TWO,ZERO DATA ZERO,ONE,TWO,THREE /0.E0,1.E0,2.E0,3.E0/ C IF (IFLAG .NE. 0) GO TO 5 C C INSERT PRINT STATEMENTS HERE WHEN NPRINT IS POSITIVE. C RETURN 5 CONTINUE DO 10 K = 1, N TEMP = (THREE - TWO*X(K))*X(K) TEMP1 = ZERO IF (K .NE. 1) TEMP1 = X(K-1) TEMP2 = ZERO IF (K .NE. N) TEMP2 = X(K+1) FVEC(K) = TEMP - TEMP1 - TWO*TEMP2 + ONE 10 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 0 Results obtained with different compilers or machines may be slightly different. 0 FINAL L2 NORM OF THE RESIDUALS 0.1192636E-07 0 NUMBER OF FUNCTION EVALUATIONS 14 1 0 Page 0 EXIT PARAMETER 1 0 FINAL APPROXIMATE SOLUTION 0 -0.5706545E+00 -0.6816283E+00 -0.7017325E+00 -0.7042129E+00 -0.7013690E+00 -0.6918656E+00 -0.6657920E+00 -0.5960342E+00 -0.4164121E+00 1 0 1 0 Page 0 Documentation for MINPACK subroutine HYBRJ1 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of HYBRJ1 is to find a zero of a system of N non- linear functions in N variables by a modification of the Powell hybrid method. This is done by using the more general nonlinea equation solver HYBRJ. The user must provide a subroutine whic calculates the functions and the Jacobian. 0 2. Subroutine and type statements. 0 SUBROUTINE HYBRJ1(FCN,N,X,FVEC,FJAC,LDFJAC,TOL,INFO,WA,LWA) INTEGER N,LDFJAC,INFO,LWA REAL TOL REAL X(N),FVEC(N),FJAC(LDFJAC,N),WA(LWA) EXTERNAL FCN 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to HYBRJ1 and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from HYBRJ1. 0 FCN is the name of the user-supplied subroutine which calculate the functions and the Jacobian. FCN must be declared in an EXTERNAL statement in the user calling program, and should be written as follows. 0 SUBROUTINE FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG) INTEGER N,LDFJAC,IFLAG REAL X(N),FVEC(N),FJAC(LDFJAC,N) ---------- IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC. IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC. ---------- RETURN END 0 The value of IFLAG should not be changed by FCN unless the 1 0 Page 0 user wants to terminate execution of HYBRJ1. In this case se IFLAG to a negative integer. 0 N is a positive integer input variable set to the number of functions and variables. 0 X is an array of length N. On input X must contain an initial estimate of the solution vector. On output X contains the final estimate of the solution vector. 0 FVEC is an output array of length N which contains the function evaluated at the output X. 0 FJAC is an output N by N array which contains the orthogonal matrix Q produced by the QR factorization of the final approx imate Jacobian. Section 6 contains more details about the approximation to the Jacobian. 0 LDFJAC is a positive integer input variable not less than N which specifies the leading dimension of the array FJAC. 0 TOL is a nonnegative input variable. Termination occurs when the algorithm estimates that the relative error between X and the solution is at most TOL. Section 4 contains more details about TOL. 0 INFO is an integer output variable. If the user has terminated execution, INFO is set to the (negative) value of IFLAG. See description of FCN. Otherwise, INFO is set as follows. 0 INFO = 0 Improper input parameters. 0 INFO = 1 Algorithm estimates that the relative error between X and the solution is at most TOL. 0 INFO = 2 Number of calls to FCN with IFLAG = 1 has reached 100*(N+1). 0 INFO = 3 TOL is too small. No further improvement in the approximate solution X is possible. 0 INFO = 4 Iteration is not making good progress. 0 Sections 4 and 5 contain more details about INFO. 0 WA is a work array of length LWA. 0 LWA is a positive integer input variable not less than (N*(N+13))/2. 0 4. Successful completion. 0 The accuracy of HYBRJ1 is controlled by the convergence 1 0 Page 0 parameter TOL. This parameter is used in a test which makes a comparison between the approximation X and a solution XSOL. HYBRJ1 terminates when the test is satisfied. If TOL is less than the machine precision (as defined by the MINPACK function SPMPAR(1)), then HYBRJ1 only attempts to satisfy the test defined by the machine precision. Further progress is not usu- ally possible. Unless high precision solutions are required, the recommended value for TOL is the square root of the machine precision. 0 The test assumes that the functions and the Jacobian are coded consistently, and that the functions are reasonably well behaved. If these conditions are not satisfied, then HYBRJ1 ma incorrectly indicate convergence. The coding of the Jacobian can be checked by the MINPACK subroutine CHKDER. If the Jaco- bian is coded correctly, then the validity of the answer can be checked, for example, by rerunning HYBRJ1 with a tighter toler- ance. 0 Convergence test. If ENORM(Z) denotes the Euclidean norm of a vector Z, then this test attempts to guarantee that 0 ENORM(X-XSOL) .LE. TOL*ENORM(XSOL). 0 If this condition is satisfied with TOL = 10**(-K), then the larger components of X have K significant decimal digits and INFO is set to 1. There is a danger that the smaller compo- nents of X may have large relative errors, but the fast rate of convergence of HYBRJ1 usually avoids this possibility. 0 5. Unsuccessful completion. 0 Unsuccessful termination of HYBRJ1 can be due to improper input parameters, arithmetic interrupts, an excessive number of func- tion evaluations, or lack of good progress. 0 Improper input parameters. INFO is set to 0 if N .LE. 0, or LDFJAC .LT. N, or TOL .LT. 0.E0, or LWA .LT. (N*(N+13))/2. 0 Arithmetic interrupts. If these interrupts occur in the FCN subroutine during an early stage of the computation, they may be caused by an unacceptable choice of X by HYBRJ1. In this case, it may be possible to remedy the situation by not evalu ating the functions here, but instead setting the components of FVEC to numbers that exceed those in the initial FVEC, thereby indirectly reducing the step length. The step length can be more directly controlled by using instead HYBRJ, which includes in its calling sequence the step-length- governing parameter FACTOR. 0 Excessive number of function evaluations. If the number of calls to FCN with IFLAG = 1 reaches 100*(N+1), then this indi cates that the routine is converging very slowly as measured 1 0 Page 0 by the progress of FVEC, and INFO is set to 2. This situatio should be unusual because, as indicated below, lack of good progress is usually diagnosed earlier by HYBRJ1, causing ter- mination with INFO = 4. 0 Lack of good progress. HYBRJ1 searches for a zero of the syste by minimizing the sum of the squares of the functions. In so doing, it can become trapped in a region where the minimum does not correspond to a zero of the system and, in this situ ation, the iteration eventually fails to make good progress. In particular, this will happen if the system does not have a zero. If the system has a zero, rerunning HYBRJ1 from a dif- ferent starting point may be helpful. 0 6. Characteristics of the algorithm. 0 HYBRJ1 is a modification of the Powell hybrid method. Two of its main characteristics involve the choice of the correction a a convex combination of the Newton and scaled gradient direc- tions, and the updating of the Jacobian by the rank-1 method of Broyden. The choice of the correction guarantees (under reason able conditions) global convergence for starting points far fro the solution and a fast rate of convergence. The Jacobian is calculated at the starting point, but it is not recalculated until the rank-1 method fails to produce satisfactory progress. 0 Timing. The time required by HYBRJ1 to solve a given problem depends on N, the behavior of the functions, the accuracy requested, and the starting point. The number of arithmetic operations needed by HYBRJ1 is about 11.5*(N**2) to process each evaluation of the functions (call to FCN with IFLAG = 1) and 1.3*(N**3) to process each evaluation of the Jacobian (call to FCN with IFLAG = 2). Unless FCN can be evaluated quickly, the timing of HYBRJ1 will be strongly influenced by the time spent in FCN. 0 Storage. HYBRJ1 requires (3*N**2 + 17*N)/2 single precision storage locations, in addition to the storage required by the program. There are no internally declared storage arrays. 0 7. Subprograms required. 0 USER-supplied ...... FCN 0 MINPACK-supplied ... DOGLEG,SPMPAR,ENORM,HYBRJ, QFORM,QRFAC,R1MPYQ,R1UPDT 0 FORTRAN-supplied ... ABS,AMAX1,AMIN1,SQRT,MIN0,MOD 0 8. References. 1 0 Page 0 M. J. D. Powell, A Hybrid Method for Nonlinear Equations. Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, editor. Gordon and Breach, 1970. 0 9. Example. 0 The problem is to determine the values of x(1), x(2), ..., x(9) which solve the system of tridiagonal equations 0 (3-2*x(1))*x(1) -2*x(2) = -1 -x(i-1) + (3-2*x(i))*x(i) -2*x(i+1) = -1, i=2-8 -x(8) + (3-2*x(9))*x(9) = -1 0 C ********** C C DRIVER FOR HYBRJ1 EXAMPLE. C SINGLE PRECISION VERSION C C ********** INTEGER J,N,LDFJAC,INFO,LWA,NWRITE REAL TOL,FNORM REAL X(9),FVEC(9),FJAC(9,9),WA(99) REAL ENORM,SPMPAR EXTERNAL FCN C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C N = 9 C C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH SOLUTION. C DO 10 J = 1, 9 X(J) = -1.E0 10 CONTINUE C LDFJAC = 9 LWA = 99 C C SET TOL TO THE SQUARE ROOT OF THE MACHINE PRECISION. C UNLESS HIGH PRECISION SOLUTIONS ARE REQUIRED, C THIS IS THE RECOMMENDED SETTING. C TOL = SQRT(SPMPAR(1)) C CALL HYBRJ1(FCN,N,X,FVEC,FJAC,LDFJAC,TOL,INFO,WA,LWA) FNORM = ENORM(N,FVEC) WRITE (NWRITE,1000) FNORM,INFO,(X(J),J=1,N) STOP 1000 FORMAT (5X,31H FINAL L2 NORM OF THE RESIDUALS,E15.7 // * 5X,15H EXIT PARAMETER,16X,I10 // * 5X,27H FINAL APPROXIMATE SOLUTION // (5X,3E15.7)) 1 0 Page 0 C C LAST CARD OF DRIVER FOR HYBRJ1 EXAMPLE. C END SUBROUTINE FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG) INTEGER N,LDFJAC,IFLAG REAL X(N),FVEC(N),FJAC(LDFJAC,N) C C SUBROUTINE FCN FOR HYBRJ1 EXAMPLE. C INTEGER J,K REAL ONE,TEMP,TEMP1,TEMP2,THREE,TWO,ZERO DATA ZERO,ONE,TWO,THREE,FOUR /0.E0,1.E0,2.E0,3.E0,4.E0/ C IF (IFLAG .EQ. 2) GO TO 20 DO 10 K = 1, N TEMP = (THREE - TWO*X(K))*X(K) TEMP1 = ZERO IF (K .NE. 1) TEMP1 = X(K-1) TEMP2 = ZERO IF (K .NE. N) TEMP2 = X(K+1) FVEC(K) = TEMP - TEMP1 - TWO*TEMP2 + ONE 10 CONTINUE GO TO 50 20 CONTINUE DO 40 K = 1, N DO 30 J = 1, N FJAC(K,J) = ZERO 30 CONTINUE FJAC(K,K) = THREE - FOUR*X(K) IF (K .NE. 1) FJAC(K,K-1) = -ONE IF (K .NE. N) FJAC(K,K+1) = -TWO 40 CONTINUE 50 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 0 Results obtained with different compilers or machines may be slightly different. 0 FINAL L2 NORM OF THE RESIDUALS 0.1192636E-07 0 EXIT PARAMETER 1 0 FINAL APPROXIMATE SOLUTION 0 -0.5706545E+00 -0.6816283E+00 -0.7017325E+00 -0.7042129E+00 -0.7013690E+00 -0.6918656E+00 -0.6657920E+00 -0.5960342E+00 -0.4164121E+00 1 0 1 0 Page 0 Documentation for MINPACK subroutine HYBRJ 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of HYBRJ is to find a zero of a system of N non- linear functions in N variables by a modification of the Powell hybrid method. The user must provide a subroutine which calcu- lates the functions and the Jacobian. 0 2. Subroutine and type statements. 0 SUBROUTINE HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,DIAG, * MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF, * WA1,WA2,WA3,WA4) INTEGER N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,LR REAL XTOL,FACTOR REAL X(N),FVEC(N),FJAC(LDFJAC,N),DIAG(N),R(LR),QTF(N), * WA1(N),WA2(N),WA3(N),WA4(N) 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to HYBRJ and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from HYBRJ. 0 FCN is the name of the user-supplied subroutine which calculate the functions and the Jacobian. FCN must be declared in an EXTERNAL statement in the user calling program, and should be written as follows. 0 SUBROUTINE FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG) INTEGER N,LDFJAC,IFLAG REAL X(N),FVEC(N),FJAC(LDFJAC,N) ---------- IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC. IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC. ---------- RETURN END 1 0 Page 0 The value of IFLAG should not be changed by FCN unless the user wants to terminate execution of HYBRJ. In this case set IFLAG to a negative integer. 0 N is a positive integer input variable set to the number of functions and variables. 0 X is an array of length N. On input X must contain an initial estimate of the solution vector. On output X contains the final estimate of the solution vector. 0 FVEC is an output array of length N which contains the function evaluated at the output X. 0 FJAC is an output N by N array which contains the orthogonal matrix Q produced by the QR factorization of the final approx imate Jacobian. Section 6 contains more details about the approximation to the Jacobian. 0 LDFJAC is a positive integer input variable not less than N which specifies the leading dimension of the array FJAC. 0 XTOL is a nonnegative input variable. Termination occurs when the relative error between two consecutive iterates is at mos XTOL. Therefore, XTOL measures the relative error desired in the approximate solution. Section 4 contains more details about XTOL. 0 MAXFEV is a positive integer input variable. Termination occur when the number of calls to FCN with IFLAG = 1 has reached MAXFEV. 0 DIAG is an array of length N. If MODE = 1 (see below), DIAG is internally set. If MODE = 2, DIAG must contain positive entries that serve as multiplicative scale factors for the variables. 0 MODE is an integer input variable. If MODE = 1, the variables will be scaled internally. If MODE = 2, the scaling is speci fied by the input DIAG. Other values of MODE are equivalent to MODE = 1. 0 FACTOR is a positive input variable used in determining the ini tial step bound. This bound is set to the product of FACTOR and the Euclidean norm of DIAG*X if nonzero, or else to FACTO itself. In most cases FACTOR should lie in the interval (.1,100.). 100. is a generally recommended value. 0 NPRINT is an integer input variable that enables controlled printing of iterates if it is positive. In this case, FCN is called with IFLAG = 0 at the beginning of the first iteration and every NPRINT iterations thereafter and immediately prior to return, with X and FVEC available for printing. FVEC and FJAC should not be altered. If NPRINT is not positive, no 1 0 Page 0 special calls of FCN with IFLAG = 0 are made. 0 INFO is an integer output variable. If the user has terminated execution, INFO is set to the (negative) value of IFLAG. See description of FCN. Otherwise, INFO is set as follows. 0 INFO = 0 Improper input parameters. 0 INFO = 1 Relative error between two consecutive iterates is at most XTOL. 0 INFO = 2 Number of calls to FCN with IFLAG = 1 has reached MAXFEV. 0 INFO = 3 XTOL is too small. No further improvement in the approximate solution X is possible. 0 INFO = 4 Iteration is not making good progress, as measured by the improvement from the last five Jacobian eval uations. 0 INFO = 5 Iteration is not making good progress, as measured by the improvement from the last ten iterations. 0 Sections 4 and 5 contain more details about INFO. 0 NFEV is an integer output variable set to the number of calls t FCN with IFLAG = 1. 0 NJEV is an integer output variable set to the number of calls t FCN with IFLAG = 2. 0 R is an output array of length LR which contains the upper triangular matrix produced by the QR factorization of the final approximate Jacobian, stored rowwise. 0 LR is a positive integer input variable not less than (N*(N+1))/2. 0 QTF is an output array of length N which contains the vector (Q transpose)*FVEC. 0 WA1, WA2, WA3, and WA4 are work arrays of length N. 0 4. Successful completion. 0 The accuracy of HYBRJ is controlled by the convergence paramete XTOL. This parameter is used in a test which makes a compariso between the approximation X and a solution XSOL. HYBRJ termi- nates when the test is satisfied. If the convergence parameter is less than the machine precision (as defined by the MINPACK function SPMPAR(1)), then HYBRJ only attempts to satisfy the test defined by the machine precision. Further progress is not 1 0 Page 0 usually possible. 0 The test assumes that the functions and the Jacobian are coded consistently, and that the functions are reasonably well behaved. If these conditions are not satisfied, then HYBRJ may incorrectly indicate convergence. The coding of the Jacobian can be checked by the MINPACK subroutine CHKDER. If the Jaco- bian is coded correctly, then the validity of the answer can be checked, for example, by rerunning HYBRJ with a tighter toler- ance. 0 Convergence test. If ENORM(Z) denotes the Euclidean norm of a vector Z and D is the diagonal matrix whose entries are defined by the array DIAG, then this test attempts to guaran- tee that 0 ENORM(D*(X-XSOL)) .LE. XTOL*ENORM(D*XSOL). 0 If this condition is satisfied with XTOL = 10**(-K), then the larger components of D*X have K significant decimal digits an INFO is set to 1. There is a danger that the smaller compo- nents of D*X may have large relative errors, but the fast rat of convergence of HYBRJ usually avoids this possibility. Unless high precision solutions are required, the recommended value for XTOL is the square root of the machine precision. 0 5. Unsuccessful completion. 0 Unsuccessful termination of HYBRJ can be due to improper input parameters, arithmetic interrupts, an excessive number of func- tion evaluations, or lack of good progress. 0 Improper input parameters. INFO is set to 0 if N .LE. 0, or LDFJAC .LT. N, or XTOL .LT. 0.E0, or MAXFEV .LE. 0, or FACTOR .LE. 0.E0, or LR .LT. (N*(N+1))/2. 0 Arithmetic interrupts. If these interrupts occur in the FCN subroutine during an early stage of the computation, they may be caused by an unacceptable choice of X by HYBRJ. In this case, it may be possible to remedy the situation by rerunning HYBRJ with a smaller value of FACTOR. 0 Excessive number of function evaluations. A reasonable value for MAXFEV is 100*(N+1). If the number of calls to FCN with IFLAG = 1 reaches MAXFEV, then this indicates that the routin is converging very slowly as measured by the progress of FVEC and INFO is set to 2. This situation should be unusual because, as indicated below, lack of good progress is usually diagnosed earlier by HYBRJ, causing termination with INFO = 4 or INFO = 5. 0 Lack of good progress. HYBRJ searches for a zero of the system by minimizing the sum of the squares of the functions. In so 1 0 Page 0 doing, it can become trapped in a region where the minimum does not correspond to a zero of the system and, in this situ ation, the iteration eventually fails to make good progress. In particular, this will happen if the system does not have a zero. If the system has a zero, rerunning HYBRJ from a dif- ferent starting point may be helpful. 0 6. Characteristics of the algorithm. 0 HYBRJ is a modification of the Powell hybrid method. Two of it main characteristics involve the choice of the correction as a convex combination of the Newton and scaled gradient directions and the updating of the Jacobian by the rank-1 method of Broy- den. The choice of the correction guarantees (under reasonable conditions) global convergence for starting points far from the solution and a fast rate of convergence. The Jacobian is calcu lated at the starting point, but it is not recalculated until the rank-1 method fails to produce satisfactory progress. 0 Timing. The time required by HYBRJ to solve a given problem depends on N, the behavior of the functions, the accuracy requested, and the starting point. The number of arithmetic operations needed by HYBRJ is about 11.5*(N**2) to process each evaluation of the functions (call to FCN with IFLAG = 1) and 1.3*(N**3) to process each evaluation of the Jacobian (call to FCN with IFLAG = 2). Unless FCN can be evaluated quickly, the timing of HYBRJ will be strongly influenced by the time spent in FCN. 0 Storage. HYBRJ requires (3*N**2 + 17*N)/2 single precision storage locations, in addition to the storage required by the program. There are no internally declared storage arrays. 0 7. Subprograms required. 0 USER-supplied ...... FCN 0 MINPACK-supplied ... DOGLEG,SPMPAR,ENORM, QFORM,QRFAC,R1MPYQ,R1UPDT 0 FORTRAN-supplied ... ABS,AMAX1,AMIN1,SQRT,MIN0,MOD 0 8. References. 0 M. J. D. Powell, A Hybrid Method for Nonlinear Equations. Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, editor. Gordon and Breach, 1970. 0 9. Example. 1 0 Page 0 The problem is to determine the values of x(1), x(2), ..., x(9) which solve the system of tridiagonal equations 0 (3-2*x(1))*x(1) -2*x(2) = -1 -x(i-1) + (3-2*x(i))*x(i) -2*x(i+1) = -1, i=2-8 -x(8) + (3-2*x(9))*x(9) = -1 0 C ********** C C DRIVER FOR HYBRJ EXAMPLE. C SINGLE PRECISION VERSION C C ********** INTEGER J,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,LR,NWRITE REAL XTOL,FACTOR,FNORM REAL X(9),FVEC(9),FJAC(9,9),DIAG(9),R(45),QTF(9), * WA1(9),WA2(9),WA3(9),WA4(9) REAL ENORM,SPMPAR EXTERNAL FCN C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C N = 9 C C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH SOLUTION. C DO 10 J = 1, 9 X(J) = -1.E0 10 CONTINUE C LDFJAC = 9 LR = 45 C C SET XTOL TO THE SQUARE ROOT OF THE MACHINE PRECISION. C UNLESS HIGH PRECISION SOLUTIONS ARE REQUIRED, C THIS IS THE RECOMMENDED SETTING. C XTOL = SQRT(SPMPAR(1)) C MAXFEV = 1000 MODE = 2 DO 20 J = 1, 9 DIAG(J) = 1.E0 20 CONTINUE FACTOR = 1.E2 NPRINT = 0 C CALL HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,DIAG, * MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF, * WA1,WA2,WA3,WA4) FNORM = ENORM(N,FVEC) WRITE (NWRITE,1000) FNORM,NFEV,NJEV,INFO,(X(J),J=1,N) 1 0 Page 0 STOP 1000 FORMAT (5X,31H FINAL L2 NORM OF THE RESIDUALS,E15.7 // * 5X,31H NUMBER OF FUNCTION EVALUATIONS,I10 // * 5X,31H NUMBER OF JACOBIAN EVALUATIONS,I10 // * 5X,15H EXIT PARAMETER,16X,I10 // * 5X,27H FINAL APPROXIMATE SOLUTION // (5X,3E15.7)) C C LAST CARD OF DRIVER FOR HYBRJ EXAMPLE. C END SUBROUTINE FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG) INTEGER N,LDFJAC,IFLAG REAL X(N),FVEC(N),FJAC(LDFJAC,N) C C SUBROUTINE FCN FOR HYBRJ EXAMPLE. C INTEGER J,K REAL ONE,TEMP,TEMP1,TEMP2,THREE,TWO,ZERO DATA ZERO,ONE,TWO,THREE,FOUR /0.E0,1.E0,2.E0,3.E0,4.E0/ C IF (IFLAG .NE. 0) GO TO 5 C C INSERT PRINT STATEMENTS HERE WHEN NPRINT IS POSITIVE. C RETURN 5 CONTINUE IF (IFLAG .EQ. 2) GO TO 20 DO 10 K = 1, N TEMP = (THREE - TWO*X(K))*X(K) TEMP1 = ZERO IF (K .NE. 1) TEMP1 = X(K-1) TEMP2 = ZERO IF (K .NE. N) TEMP2 = X(K+1) FVEC(K) = TEMP - TEMP1 - TWO*TEMP2 + ONE 10 CONTINUE GO TO 50 20 CONTINUE DO 40 K = 1, N DO 30 J = 1, N FJAC(K,J) = ZERO 30 CONTINUE FJAC(K,K) = THREE - FOUR*X(K) IF (K .NE. 1) FJAC(K,K-1) = -ONE IF (K .NE. N) FJAC(K,K+1) = -TWO 40 CONTINUE 50 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 0 Results obtained with different compilers or machines may be slightly different. 1 0 Page 0 FINAL L2 NORM OF THE RESIDUALS 0.1192636E-07 0 NUMBER OF FUNCTION EVALUATIONS 11 0 NUMBER OF JACOBIAN EVALUATIONS 1 0 EXIT PARAMETER 1 0 FINAL APPROXIMATE SOLUTION 0 -0.5706545E+00 -0.6816283E+00 -0.7017325E+00 -0.7042129E+00 -0.7013690E+00 -0.6918656E+00 -0.6657920E+00 -0.5960342E+00 -0.4164121E+00 1 0 1 0 Page 0 Documentation for MINPACK subroutine LMDER1 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of LMDER1 is to minimize the sum of the squares of nonlinear functions in N variables by a modification of the Levenberg-Marquardt algorithm. This is done by using the more general least-squares solver LMDER. The user must provide a subroutine which calculates the functions and the Jacobian. 0 2. Subroutine and type statements. 0 SUBROUTINE LMDER1(FCN,M,N,X,FVEC,FJAC,LDFJAC,TOL, * INFO,IPVT,WA,LWA) INTEGER M,N,LDFJAC,INFO,LWA INTEGER IPVT(N) REAL TOL REAL X(N),FVEC(M),FJAC(LDFJAC,N),WA(LWA) EXTERNAL FCN 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to LMDER1 and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from LMDER1. 0 FCN is the name of the user-supplied subroutine which calculate the functions and the Jacobian. FCN must be declared in an EXTERNAL statement in the user calling program, and should be written as follows. 0 SUBROUTINE FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) INTEGER M,N,LDFJAC,IFLAG REAL X(N),FVEC(M),FJAC(LDFJAC,N) ---------- IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC. IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC. ---------- RETURN END 1 0 Page 0 The value of IFLAG should not be changed by FCN unless the user wants to terminate execution of LMDER1. In this case se IFLAG to a negative integer. 0 M is a positive integer input variable set to the number of functions. 0 N is a positive integer input variable set to the number of variables. N must not exceed M. 0 X is an array of length N. On input X must contain an initial estimate of the solution vector. On output X contains the final estimate of the solution vector. 0 FVEC is an output array of length M which contains the function evaluated at the output X. 0 FJAC is an output M by N array. The upper N by N submatrix of FJAC contains an upper triangular matrix R with diagonal ele- ments of nonincreasing magnitude such that 0 T T T P *(JAC *JAC)*P = R *R, 0 where P is a permutation matrix and JAC is the final calcu- lated Jacobian. Column j of P is column IPVT(j) (see below) of the identity matrix. The lower trapezoidal part of FJAC contains information generated during the computation of R. 0 LDFJAC is a positive integer input variable not less than M which specifies the leading dimension of the array FJAC. 0 TOL is a nonnegative input variable. Termination occurs when the algorithm estimates either that the relative error in the sum of squares is at most TOL or that the relative error between X and the solution is at most TOL. Section 4 contain more details about TOL. 0 INFO is an integer output variable. If the user has terminated execution, INFO is set to the (negative) value of IFLAG. See description of FCN. Otherwise, INFO is set as follows. 0 INFO = 0 Improper input parameters. 0 INFO = 1 Algorithm estimates that the relative error in the sum of squares is at most TOL. 0 INFO = 2 Algorithm estimates that the relative error between X and the solution is at most TOL. 0 INFO = 3 Conditions for INFO = 1 and INFO = 2 both hold. 0 INFO = 4 FVEC is orthogonal to the columns of the Jacobian t machine precision. 1 0 Page 0 INFO = 5 Number of calls to FCN with IFLAG = 1 has reached 100*(N+1). 0 INFO = 6 TOL is too small. No further reduction in the sum of squares is possible. 0 INFO = 7 TOL is too small. No further improvement in the approximate solution X is possible. 0 Sections 4 and 5 contain more details about INFO. 0 IPVT is an integer output array of length N. IPVT defines a permutation matrix P such that JAC*P = Q*R, where JAC is the final calculated Jacobian, Q is orthogonal (not stored), and is upper triangular with diagonal elements of nonincreasing magnitude. Column j of P is column IPVT(j) of the identity matrix. 0 WA is a work array of length LWA. 0 LWA is a positive integer input variable not less than 5*N+M. 0 4. Successful completion. 0 The accuracy of LMDER1 is controlled by the convergence parame- ter TOL. This parameter is used in tests which make three type of comparisons between the approximation X and a solution XSOL. LMDER1 terminates when any of the tests is satisfied. If TOL i less than the machine precision (as defined by the MINPACK func tion SPMPAR(1)), then LMDER1 only attempts to satisfy the test defined by the machine precision. Further progress is not usu- ally possible. Unless high precision solutions are required, the recommended value for TOL is the square root of the machine precision. 0 The tests assume that the functions and the Jacobian are coded consistently, and that the functions are reasonably well behaved. If these conditions are not satisfied, then LMDER1 ma incorrectly indicate convergence. The coding of the Jacobian can be checked by the MINPACK subroutine CHKDER. If the Jaco- bian is coded correctly, then the validity of the answer can be checked, for example, by rerunning LMDER1 with a tighter toler- ance. 0 First convergence test. If ENORM(Z) denotes the Euclidean norm of a vector Z, then this test attempts to guarantee that 0 ENORM(FVEC) .LE. (1+TOL)*ENORM(FVECS), 0 where FVECS denotes the functions evaluated at XSOL. If this condition is satisfied with TOL = 10**(-K), then the final residual norm ENORM(FVEC) has K significant decimal digits an INFO is set to 1 (or to 3 if the second test is also 1 0 Page 0 satisfied). 0 Second convergence test. If D is a diagonal matrix (implicitly generated by LMDER1) whose entries contain scale factors for the variables, then this test attempts to guarantee that 0 ENORM(D*(X-XSOL)) .LE. TOL*ENORM(D*XSOL). 0 If this condition is satisfied with TOL = 10**(-K), then the larger components of D*X have K significant decimal digits an INFO is set to 2 (or to 3 if the first test is also satis- fied). There is a danger that the smaller components of D*X may have large relative errors, but the choice of D is such that the accuracy of the components of X is usually related t their sensitivity. 0 Third convergence test. This test is satisfied when FVEC is orthogonal to the columns of the Jacobian to machine preci- sion. There is no clear relationship between this test and the accuracy of LMDER1, and furthermore, the test is equally well satisfied at other critical points, namely maximizers an saddle points. Therefore, termination caused by this test (INFO = 4) should be examined carefully. 0 5. Unsuccessful completion. 0 Unsuccessful termination of LMDER1 can be due to improper input parameters, arithmetic interrupts, or an excessive number of function evaluations. 0 Improper input parameters. INFO is set to 0 if N .LE. 0, or M .LT. N, or LDFJAC .LT. M, or TOL .LT. 0.E0, or LWA .LT. 5*N+M. 0 Arithmetic interrupts. If these interrupts occur in the FCN subroutine during an early stage of the computation, they may be caused by an unacceptable choice of X by LMDER1. In this case, it may be possible to remedy the situation by not evalu ating the functions here, but instead setting the components of FVEC to numbers that exceed those in the initial FVEC, thereby indirectly reducing the step length. The step length can be more directly controlled by using instead LMDER, which includes in its calling sequence the step-length- governing parameter FACTOR. 0 Excessive number of function evaluations. If the number of calls to FCN with IFLAG = 1 reaches 100*(N+1), then this indi cates that the routine is converging very slowly as measured by the progress of FVEC, and INFO is set to 5. In this case, it may be helpful to restart LMDER1, thereby forcing it to disregard old (and possibly harmful) information. 0 1 0 Page 0 6. Characteristics of the algorithm. 0 LMDER1 is a modification of the Levenberg-Marquardt algorithm. Two of its main characteristics involve the proper use of implicitly scaled variables and an optimal choice for the cor- rection. The use of implicitly scaled variables achieves scale invariance of LMDER1 and limits the size of the correction in any direction where the functions are changing rapidly. The optimal choice of the correction guarantees (under reasonable conditions) global convergence from starting points far from th solution and a fast rate of convergence for problems with small residuals. 0 Timing. The time required by LMDER1 to solve a given problem depends on M and N, the behavior of the functions, the accu- racy requested, and the starting point. The number of arith- metic operations needed by LMDER1 is about N**3 to process each evaluation of the functions (call to FCN with IFLAG = 1) and M*(N**2) to process each evaluation of the Jacobian (call to FCN with IFLAG = 2). Unless FCN can be evaluated quickly, the timing of LMDER1 will be strongly influenced by the time spent in FCN. 0 Storage. LMDER1 requires M*N + 2*M + 6*N single precision sto- rage locations and N integer storage locations, in addition t the storage required by the program. There are no internally declared storage arrays. 0 7. Subprograms required. 0 USER-supplied ...... FCN 0 MINPACK-supplied ... SPMPAR,ENORM,LMDER,LMPAR,QRFAC,QRSOLV 0 FORTRAN-supplied ... ABS,AMAX1,AMIN1,SQRT,MOD 0 8. References. 0 Jorge J. More, The Levenberg-Marquardt Algorithm, Implementatio and Theory. Numerical Analysis, G. A. Watson, editor. Lecture Notes in Mathematics 630, Springer-Verlag, 1977. 0 9. Example. 0 The problem is to determine the values of x(1), x(2), and x(3) which provide the best fit (in the least squares sense) of 0 x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)), i = 1, 15 0 to the data 1 0 Page 0 y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39, 0.37,0.58,0.73,0.96,1.34,2.10,4.39), 0 where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)). The i-th component of FVEC is thus defined by 0 y(i) - (x(1) + u(i)/(v(i)*x(2) + w(i)*x(3))). 0 C ********** C C DRIVER FOR LMDER1 EXAMPLE. C SINGLE PRECISION VERSION C C ********** INTEGER J,M,N,LDFJAC,INFO,LWA,NWRITE INTEGER IPVT(3) REAL TOL,FNORM REAL X(3),FVEC(15),FJAC(15,3),WA(30) REAL ENORM,SPMPAR EXTERNAL FCN C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C M = 15 N = 3 C C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH FIT. C X(1) = 1.E0 X(2) = 1.E0 X(3) = 1.E0 C LDFJAC = 15 LWA = 30 C C SET TOL TO THE SQUARE ROOT OF THE MACHINE PRECISION. C UNLESS HIGH PRECISION SOLUTIONS ARE REQUIRED, C THIS IS THE RECOMMENDED SETTING. C TOL = SQRT(SPMPAR(1)) C CALL LMDER1(FCN,M,N,X,FVEC,FJAC,LDFJAC,TOL, * INFO,IPVT,WA,LWA) FNORM = ENORM(M,FVEC) WRITE (NWRITE,1000) FNORM,INFO,(X(J),J=1,N) STOP 1000 FORMAT (5X,31H FINAL L2 NORM OF THE RESIDUALS,E15.7 // * 5X,15H EXIT PARAMETER,16X,I10 // * 5X,27H FINAL APPROXIMATE SOLUTION // 5X,3E15.7) C C LAST CARD OF DRIVER FOR LMDER1 EXAMPLE. C 1 0 Page 0 END SUBROUTINE FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) INTEGER M,N,LDFJAC,IFLAG REAL X(N),FVEC(M),FJAC(LDFJAC,N) C C SUBROUTINE FCN FOR LMDER1 EXAMPLE. C INTEGER I REAL TMP1,TMP2,TMP3,TMP4 REAL Y(15) DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8), * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15) * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1, * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/ C IF (IFLAG .EQ. 2) GO TO 20 DO 10 I = 1, 15 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3)) 10 CONTINUE GO TO 40 20 CONTINUE DO 30 I = 1, 15 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2 FJAC(I,1) = -1.E0 FJAC(I,2) = TMP1*TMP2/TMP4 FJAC(I,3) = TMP1*TMP3/TMP4 30 CONTINUE 40 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 0 Results obtained with different compilers or machines may be slightly different. 0 FINAL L2 NORM OF THE RESIDUALS 0.9063596E-01 0 EXIT PARAMETER 1 0 FINAL APPROXIMATE SOLUTION 0 0.8241058E-01 0.1133037E+01 0.2343695E+01 1 0 1 0 Page 0 Documentation for MINPACK subroutine LMDER 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of LMDER is to minimize the sum of the squares of M nonlinear functions in N variables by a modification of the Levenberg-Marquardt algorithm. The user must provide a subrou- tine which calculates the functions and the Jacobian. 0 2. Subroutine and type statements. 0 SUBROUTINE LMDER(FCN,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL, * MAXFEV,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,NJEV, * IPVT,QTF,WA1,WA2,WA3,WA4) INTEGER M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV INTEGER IPVT(N) REAL FTOL,XTOL,GTOL,FACTOR REAL X(N),FVEC(M),FJAC(LDFJAC,N),DIAG(N),QTF(N), * WA1(N),WA2(N),WA3(N),WA4(M) 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to LMDER and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from LMDER. 0 FCN is the name of the user-supplied subroutine which calculate the functions and the Jacobian. FCN must be declared in an EXTERNAL statement in the user calling program, and should be written as follows. 0 SUBROUTINE FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) INTEGER M,N,LDFJAC,IFLAG REAL X(N),FVEC(M),FJAC(LDFJAC,N) ---------- IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC. IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC. ---------- RETURN END 1 0 Page 0 The value of IFLAG should not be changed by FCN unless the user wants to terminate execution of LMDER. In this case set IFLAG to a negative integer. 0 M is a positive integer input variable set to the number of functions. 0 N is a positive integer input variable set to the number of variables. N must not exceed M. 0 X is an array of length N. On input X must contain an initial estimate of the solution vector. On output X contains the final estimate of the solution vector. 0 FVEC is an output array of length M which contains the function evaluated at the output X. 0 FJAC is an output M by N array. The upper N by N submatrix of FJAC contains an upper triangular matrix R with diagonal ele- ments of nonincreasing magnitude such that 0 T T T P *(JAC *JAC)*P = R *R, 0 where P is a permutation matrix and JAC is the final calcu- lated Jacobian. Column j of P is column IPVT(j) (see below) of the identity matrix. The lower trapezoidal part of FJAC contains information generated during the computation of R. 0 LDFJAC is a positive integer input variable not less than M which specifies the leading dimension of the array FJAC. 0 FTOL is a nonnegative input variable. Termination occurs when both the actual and predicted relative reductions in the sum of squares are at most FTOL. Therefore, FTOL measures the relative error desired in the sum of squares. Section 4 con- tains more details about FTOL. 0 XTOL is a nonnegative input variable. Termination occurs when the relative error between two consecutive iterates is at mos XTOL. Therefore, XTOL measures the relative error desired in the approximate solution. Section 4 contains more details about XTOL. 0 GTOL is a nonnegative input variable. Termination occurs when the cosine of the angle between FVEC and any column of the Jacobian is at most GTOL in absolute value. Therefore, GTOL measures the orthogonality desired between the function vecto and the columns of the Jacobian. Section 4 contains more details about GTOL. 0 MAXFEV is a positive integer input variable. Termination occur when the number of calls to FCN with IFLAG = 1 has reached MAXFEV. 1 0 Page 0 DIAG is an array of length N. If MODE = 1 (see below), DIAG is internally set. If MODE = 2, DIAG must contain positive entries that serve as multiplicative scale factors for the variables. 0 MODE is an integer input variable. If MODE = 1, the variables will be scaled internally. If MODE = 2, the scaling is speci fied by the input DIAG. Other values of MODE are equivalent to MODE = 1. 0 FACTOR is a positive input variable used in determining the ini tial step bound. This bound is set to the product of FACTOR and the Euclidean norm of DIAG*X if nonzero, or else to FACTO itself. In most cases FACTOR should lie in the interval (.1,100.). 100. is a generally recommended value. 0 NPRINT is an integer input variable that enables controlled printing of iterates if it is positive. In this case, FCN is called with IFLAG = 0 at the beginning of the first iteration and every NPRINT iterations thereafter and immediately prior to return, with X, FVEC, and FJAC available for printing. FVEC and FJAC should not be altered. If NPRINT is not posi- tive, no special calls of FCN with IFLAG = 0 are made. 0 INFO is an integer output variable. If the user has terminated execution, INFO is set to the (negative) value of IFLAG. See description of FCN. Otherwise, INFO is set as follows. 0 INFO = 0 Improper input parameters. 0 INFO = 1 Both actual and predicted relative reductions in th sum of squares are at most FTOL. 0 INFO = 2 Relative error between two consecutive iterates is at most XTOL. 0 INFO = 3 Conditions for INFO = 1 and INFO = 2 both hold. 0 INFO = 4 The cosine of the angle between FVEC and any column of the Jacobian is at most GTOL in absolute value. 0 INFO = 5 Number of calls to FCN with IFLAG = 1 has reached MAXFEV. 0 INFO = 6 FTOL is too small. No further reduction in the sum of squares is possible. 0 INFO = 7 XTOL is too small. No further improvement in the approximate solution X is possible. 0 INFO = 8 GTOL is too small. FVEC is orthogonal to the columns of the Jacobian to machine precision. 0 Sections 4 and 5 contain more details about INFO. 1 0 Page 0 NFEV is an integer output variable set to the number of calls t FCN with IFLAG = 1. 0 NJEV is an integer output variable set to the number of calls t FCN with IFLAG = 2. 0 IPVT is an integer output array of length N. IPVT defines a permutation matrix P such that JAC*P = Q*R, where JAC is the final calculated Jacobian, Q is orthogonal (not stored), and is upper triangular with diagonal elements of nonincreasing magnitude. Column j of P is column IPVT(j) of the identity matrix. 0 QTF is an output array of length N which contains the first N elements of the vector (Q transpose)*FVEC. 0 WA1, WA2, and WA3 are work arrays of length N. 0 WA4 is a work array of length M. 0 4. Successful completion. 0 The accuracy of LMDER is controlled by the convergence parame- ters FTOL, XTOL, and GTOL. These parameters are used in tests which make three types of comparisons between the approximation X and a solution XSOL. LMDER terminates when any of the tests is satisfied. If any of the convergence parameters is less tha the machine precision (as defined by the MINPACK function SPMPAR(1)), then LMDER only attempts to satisfy the test define by the machine precision. Further progress is not usually pos- sible. 0 The tests assume that the functions and the Jacobian are coded consistently, and that the functions are reasonably well behaved. If these conditions are not satisfied, then LMDER may incorrectly indicate convergence. The coding of the Jacobian can be checked by the MINPACK subroutine CHKDER. If the Jaco- bian is coded correctly, then the validity of the answer can be checked, for example, by rerunning LMDER with tighter toler- ances. 0 First convergence test. If ENORM(Z) denotes the Euclidean norm of a vector Z, then this test attempts to guarantee that 0 ENORM(FVEC) .LE. (1+FTOL)*ENORM(FVECS), 0 where FVECS denotes the functions evaluated at XSOL. If this condition is satisfied with FTOL = 10**(-K), then the final residual norm ENORM(FVEC) has K significant decimal digits an INFO is set to 1 (or to 3 if the second test is also satis- fied). Unless high precision solutions are required, the recommended value for FTOL is the square root of the machine precision. 1 0 Page 0 Second convergence test. If D is the diagonal matrix whose entries are defined by the array DIAG, then this test attempt to guarantee that 0 ENORM(D*(X-XSOL)) .LE. XTOL*ENORM(D*XSOL). 0 If this condition is satisfied with XTOL = 10**(-K), then the larger components of D*X have K significant decimal digits an INFO is set to 2 (or to 3 if the first test is also satis- fied). There is a danger that the smaller components of D*X may have large relative errors, but if MODE = 1, then the accuracy of the components of X is usually related to their sensitivity. Unless high precision solutions are required, the recommended value for XTOL is the square root of the machine precision. 0 Third convergence test. This test is satisfied when the cosine of the angle between FVEC and any column of the Jacobian at X is at most GTOL in absolute value. There is no clear rela- tionship between this test and the accuracy of LMDER, and furthermore, the test is equally well satisfied at other crit ical points, namely maximizers and saddle points. Therefore, termination caused by this test (INFO = 4) should be examined carefully. The recommended value for GTOL is zero. 0 5. Unsuccessful completion. 0 Unsuccessful termination of LMDER can be due to improper input parameters, arithmetic interrupts, or an excessive number of function evaluations. 0 Improper input parameters. INFO is set to 0 if N .LE. 0, or M .LT. N, or LDFJAC .LT. M, or FTOL .LT. 0.E0, or XTOL .LT. 0.E0, or GTOL .LT. 0.E0, or MAXFEV .LE. 0, or FACTOR .LE. 0.E0. 0 Arithmetic interrupts. If these interrupts occur in the FCN subroutine during an early stage of the computation, they may be caused by an unacceptable choice of X by LMDER. In this case, it may be possible to remedy the situation by rerunning LMDER with a smaller value of FACTOR. 0 Excessive number of function evaluations. A reasonable value for MAXFEV is 100*(N+1). If the number of calls to FCN with IFLAG = 1 reaches MAXFEV, then this indicates that the routin is converging very slowly as measured by the progress of FVEC and INFO is set to 5. In this case, it may be helpful to restart LMDER with MODE set to 1. 0 6. Characteristics of the algorithm. 0 LMDER is a modification of the Levenberg-Marquardt algorithm. 1 0 Page 0 Two of its main characteristics involve the proper use of implicitly scaled variables (if MODE = 1) and an optimal choice for the correction. The use of implicitly scaled variables achieves scale invariance of LMDER and limits the size of the correction in any direction where the functions are changing rapidly. The optimal choice of the correction guarantees (unde reasonable conditions) global convergence from starting points far from the solution and a fast rate of convergence for prob- lems with small residuals. 0 Timing. The time required by LMDER to solve a given problem depends on M and N, the behavior of the functions, the accu- racy requested, and the starting point. The number of arith- metic operations needed by LMDER is about N**3 to process eac evaluation of the functions (call to FCN with IFLAG = 1) and M*(N**2) to process each evaluation of the Jacobian (call to FCN with IFLAG = 2). Unless FCN can be evaluated quickly, th timing of LMDER will be strongly influenced by the time spent in FCN. 0 Storage. LMDER requires M*N + 2*M + 6*N single precision sto- rage locations and N integer storage locations, in addition t the storage required by the program. There are no internally declared storage arrays. 0 7. Subprograms required. 0 USER-supplied ...... FCN 0 MINPACK-supplied ... SPMPAR,ENORM,LMPAR,QRFAC,QRSOLV 0 FORTRAN-supplied ... ABS,AMAX1,AMIN1,SQRT,MOD 0 8. References. 0 Jorge J. More, The Levenberg-Marquardt Algorithm, Implementatio and Theory. Numerical Analysis, G. A. Watson, editor. Lecture Notes in Mathematics 630, Springer-Verlag, 1977. 0 9. Example. 0 The problem is to determine the values of x(1), x(2), and x(3) which provide the best fit (in the least squares sense) of 0 x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)), i = 1, 15 0 to the data 0 y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39, 0.37,0.58,0.73,0.96,1.34,2.10,4.39), 1 0 Page 0 where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)). The i-th component of FVEC is thus defined by 0 y(i) - (x(1) + u(i)/(v(i)*x(2) + w(i)*x(3))). 0 C ********** C C DRIVER FOR LMDER EXAMPLE. C SINGLE PRECISION VERSION C C ********** INTEGER J,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,NWRITE INTEGER IPVT(3) REAL FTOL,XTOL,GTOL,FACTOR,FNORM REAL X(3),FVEC(15),FJAC(15,3),DIAG(3),QTF(3), * WA1(3),WA2(3),WA3(3),WA4(15) REAL ENORM,SPMPAR EXTERNAL FCN C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C M = 15 N = 3 C C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH FIT. C X(1) = 1.E0 X(2) = 1.E0 X(3) = 1.E0 C LDFJAC = 15 C C SET FTOL AND XTOL TO THE SQUARE ROOT OF THE MACHINE PRECISION C AND GTOL TO ZERO. UNLESS HIGH PRECISION SOLUTIONS ARE C REQUIRED, THESE ARE THE RECOMMENDED SETTINGS. C FTOL = SQRT(SPMPAR(1)) XTOL = SQRT(SPMPAR(1)) GTOL = 0.E0 C MAXFEV = 400 MODE = 1 FACTOR = 1.E2 NPRINT = 0 C CALL LMDER(FCN,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL, * MAXFEV,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,NJEV, * IPVT,QTF,WA1,WA2,WA3,WA4) FNORM = ENORM(M,FVEC) WRITE (NWRITE,1000) FNORM,NFEV,NJEV,INFO,(X(J),J=1,N) STOP 1000 FORMAT (5X,31H FINAL L2 NORM OF THE RESIDUALS,E15.7 // 1 0 Page 0 * 5X,31H NUMBER OF FUNCTION EVALUATIONS,I10 // * 5X,31H NUMBER OF JACOBIAN EVALUATIONS,I10 // * 5X,15H EXIT PARAMETER,16X,I10 // * 5X,27H FINAL APPROXIMATE SOLUTION // 5X,3E15.7) C C LAST CARD OF DRIVER FOR LMDER EXAMPLE. C END SUBROUTINE FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) INTEGER M,N,LDFJAC,IFLAG REAL X(N),FVEC(M),FJAC(LDFJAC,N) C C SUBROUTINE FCN FOR LMDER EXAMPLE. C INTEGER I REAL TMP1,TMP2,TMP3,TMP4 REAL Y(15) DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8), * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15) * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1, * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/ C IF (IFLAG .NE. 0) GO TO 5 C C INSERT PRINT STATEMENTS HERE WHEN NPRINT IS POSITIVE. C RETURN 5 CONTINUE IF (IFLAG .EQ. 2) GO TO 20 DO 10 I = 1, 15 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3)) 10 CONTINUE GO TO 40 20 CONTINUE DO 30 I = 1, 15 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2 FJAC(I,1) = -1.E0 FJAC(I,2) = TMP1*TMP2/TMP4 FJAC(I,3) = TMP1*TMP3/TMP4 30 CONTINUE 40 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 1 0 Page 0 Results obtained with different compilers or machines may be slightly different. 0 FINAL L2 NORM OF THE RESIDUALS 0.9063596E-01 0 NUMBER OF FUNCTION EVALUATIONS 6 0 NUMBER OF JACOBIAN EVALUATIONS 5 0 EXIT PARAMETER 1 0 FINAL APPROXIMATE SOLUTION 0 0.8241058E-01 0.1133037E+01 0.2343695E+01 1 0 1 0 Page 0 Documentation for MINPACK subroutine LMSTR1 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of LMSTR1 is to minimize the sum of the squares of nonlinear functions in N variables by a modification of the Levenberg-Marquardt algorithm which uses minimal storage. This is done by using the more general least-squares solver LMSTR. The user must provide a subroutine which calculates the func- tions and the rows of the Jacobian. 0 2. Subroutine and type statements. 0 SUBROUTINE LMSTR1(FCN,M,N,X,FVEC,FJAC,LDFJAC,TOL, * INFO,IPVT,WA,LWA) INTEGER M,N,LDFJAC,INFO,LWA INTEGER IPVT(N) REAL TOL REAL X(N),FVEC(M),FJAC(LDFJAC,N),WA(LWA) EXTERNAL FCN 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to LMSTR1 and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from LMSTR1. 0 FCN is the name of the user-supplied subroutine which calculate the functions and the rows of the Jacobian. FCN must be declared in an EXTERNAL statement in the user calling program and should be written as follows. 0 SUBROUTINE FCN(M,N,X,FVEC,FJROW,IFLAG) INTEGER M,N,IFLAG REAL X(N),FVEC(M),FJROW(N) ---------- IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND RETURN THIS VECTOR IN FVEC. IF IFLAG = I CALCULATE THE (I-1)-ST ROW OF THE JACOBIAN AT X AND RETURN THIS VECTOR IN FJROW. ---------- RETURN 1 0 Page 0 END 0 The value of IFLAG should not be changed by FCN unless the user wants to terminate execution of LMSTR1. In this case se IFLAG to a negative integer. 0 M is a positive integer input variable set to the number of functions. 0 N is a positive integer input variable set to the number of variables. N must not exceed M. 0 X is an array of length N. On input X must contain an initial estimate of the solution vector. On output X contains the final estimate of the solution vector. 0 FVEC is an output array of length M which contains the function evaluated at the output X. 0 FJAC is an output N by N array. The upper triangle of FJAC con tains an upper triangular matrix R such that 0 T T T P *(JAC *JAC)*P = R *R, 0 where P is a permutation matrix and JAC is the final calcu- lated Jacobian. Column j of P is column IPVT(j) (see below) of the identity matrix. The lower triangular part of FJAC contains information generated during the computation of R. 0 LDFJAC is a positive integer input variable not less than N which specifies the leading dimension of the array FJAC. 0 TOL is a nonnegative input variable. Termination occurs when the algorithm estimates either that the relative error in the sum of squares is at most TOL or that the relative error between X and the solution is at most TOL. Section 4 contain more details about TOL. 0 INFO is an integer output variable. If the user has terminated execution, INFO is set to the (negative) value of IFLAG. See description of FCN. Otherwise, INFO is set as follows. 0 INFO = 0 Improper input parameters. 0 INFO = 1 Algorithm estimates that the relative error in the sum of squares is at most TOL. 0 INFO = 2 Algorithm estimates that the relative error between X and the solution is at most TOL. 0 INFO = 3 Conditions for INFO = 1 and INFO = 2 both hold. 0 INFO = 4 FVEC is orthogonal to the columns of the Jacobian t 1 0 Page 0 machine precision. 0 INFO = 5 Number of calls to FCN with IFLAG = 1 has reached 100*(N+1). 0 INFO = 6 TOL is too small. No further reduction in the sum of squares is possible. 0 INFO = 7 TOL is too small. No further improvement in the approximate solution X is possible. 0 Sections 4 and 5 contain more details about INFO. 0 IPVT is an integer output array of length N. IPVT defines a permutation matrix P such that JAC*P = Q*R, where JAC is the final calculated Jacobian, Q is orthogonal (not stored), and is upper triangular. Column j of P is column IPVT(j) of the identity matrix. 0 WA is a work array of length LWA. 0 LWA is a positive integer input variable not less than 5*N+M. 0 4. Successful completion. 0 The accuracy of LMSTR1 is controlled by the convergence parame- ter TOL. This parameter is used in tests which make three type of comparisons between the approximation X and a solution XSOL. LMSTR1 terminates when any of the tests is satisfied. If TOL i less than the machine precision (as defined by the MINPACK func tion SPMPAR(1)), then LMSTR1 only attempts to satisfy the test defined by the machine precision. Further progress is not usu- ally possible. Unless high precision solutions are required, the recommended value for TOL is the square root of the machine precision. 0 The tests assume that the functions and the Jacobian are coded consistently, and that the functions are reasonably well behaved. If these conditions are not satisfied, then LMSTR1 ma incorrectly indicate convergence. The coding of the Jacobian can be checked by the MINPACK subroutine CHKDER. If the Jaco- bian is coded correctly, then the validity of the answer can be checked, for example, by rerunning LMSTR1 with a tighter toler- ance. 0 First convergence test. If ENORM(Z) denotes the Euclidean norm of a vector Z, then this test attempts to guarantee that 0 ENORM(FVEC) .LE. (1+TOL)*ENORM(FVECS), 0 where FVECS denotes the functions evaluated at XSOL. If this condition is satisfied with TOL = 10**(-K), then the final residual norm ENORM(FVEC) has K significant decimal digits an 1 0 Page 0 INFO is set to 1 (or to 3 if the second test is also satis- fied). 0 Second convergence test. If D is a diagonal matrix (implicitly generated by LMSTR1) whose entries contain scale factors for the variables, then this test attempts to guarantee that 0 ENORM(D*(X-XSOL)) .LE. TOL*ENORM(D*XSOL). 0 If this condition is satisfied with TOL = 10**(-K), then the larger components of D*X have K significant decimal digits an INFO is set to 2 (or to 3 if the first test is also satis- fied). There is a danger that the smaller components of D*X may have large relative errors, but the choice of D is such that the accuracy of the components of X is usually related t their sensitivity. 0 Third convergence test. This test is satisfied when FVEC is orthogonal to the columns of the Jacobian to machine preci- sion. There is no clear relationship between this test and the accuracy of LMSTR1, and furthermore, the test is equally well satisfied at other critical points, namely maximizers an saddle points. Therefore, termination caused by this test (INFO = 4) should be examined carefully. 0 5. Unsuccessful completion. 0 Unsuccessful termination of LMSTR1 can be due to improper input parameters, arithmetic interrupts, or an excessive number of function evaluations. 0 Improper input parameters. INFO is set to 0 if N .LE. 0, or M .LT. N, or LDFJAC .LT. N, or TOL .LT. 0.E0, or LWA .LT. 5*N+M. 0 Arithmetic interrupts. If these interrupts occur in the FCN subroutine during an early stage of the computation, they may be caused by an unacceptable choice of X by LMSTR1. In this case, it may be possible to remedy the situation by not evalu ating the functions here, but instead setting the components of FVEC to numbers that exceed those in the initial FVEC, thereby indirectly reducing the step length. The step length can be more directly controlled by using instead LMSTR, which includes in its calling sequence the step-length- governing parameter FACTOR. 0 Excessive number of function evaluations. If the number of calls to FCN with IFLAG = 1 reaches 100*(N+1), then this indi cates that the routine is converging very slowly as measured by the progress of FVEC, and INFO is set to 5. In this case, it may be helpful to restart LMSTR1, thereby forcing it to disregard old (and possibly harmful) information. 1 0 Page 0 6. Characteristics of the algorithm. 0 LMSTR1 is a modification of the Levenberg-Marquardt algorithm. Two of its main characteristics involve the proper use of implicitly scaled variables and an optimal choice for the cor- rection. The use of implicitly scaled variables achieves scale invariance of LMSTR1 and limits the size of the correction in any direction where the functions are changing rapidly. The optimal choice of the correction guarantees (under reasonable conditions) global convergence from starting points far from th solution and a fast rate of convergence for problems with small residuals. 0 Timing. The time required by LMSTR1 to solve a given problem depends on M and N, the behavior of the functions, the accu- racy requested, and the starting point. The number of arith- metic operations needed by LMSTR1 is about N**3 to process each evaluation of the functions (call to FCN with IFLAG = 1) and 1.5*(N**2) to process each row of the Jacobian (call to FCN with IFLAG .GE. 2). Unless FCN can be evaluated quickly, the timing of LMSTR1 will be strongly influenced by the time spent in FCN. 0 Storage. LMSTR1 requires N**2 + 2*M + 6*N single precision sto rage locations and N integer storage locations, in addition t the storage required by the program. There are no internally declared storage arrays. 0 7. Subprograms required. 0 USER-supplied ...... FCN 0 MINPACK-supplied ... SPMPAR,ENORM,LMSTR,LMPAR,QRFAC,QRSOLV, RWUPDT 0 FORTRAN-supplied ... ABS,AMAX1,AMIN1,SQRT,MOD 0 8. References. 0 Jorge J. More, The Levenberg-Marquardt Algorithm, Implementatio and Theory. Numerical Analysis, G. A. Watson, editor. Lecture Notes in Mathematics 630, Springer-Verlag, 1977. 0 9. Example. 0 The problem is to determine the values of x(1), x(2), and x(3) which provide the best fit (in the least squares sense) of 0 x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)), i = 1, 15 1 0 Page 0 to the data 0 y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39, 0.37,0.58,0.73,0.96,1.34,2.10,4.39), 0 where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)). The i-th component of FVEC is thus defined by 0 y(i) - (x(1) + u(i)/(v(i)*x(2) + w(i)*x(3))). 0 C ********** C C DRIVER FOR LMSTR1 EXAMPLE. C SINGLE PRECISION VERSION C C ********** INTEGER J,M,N,LDFJAC,INFO,LWA,NWRITE INTEGER IPVT(3) REAL TOL,FNORM REAL X(3),FVEC(15),FJAC(3,3),WA(30) REAL ENORM,SPMPAR EXTERNAL FCN C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C M = 15 N = 3 C C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH FIT. C X(1) = 1.E0 X(2) = 1.E0 X(3) = 1.E0 C LDFJAC = 3 LWA = 30 C C SET TOL TO THE SQUARE ROOT OF THE MACHINE PRECISION. C UNLESS HIGH PRECISION SOLUTIONS ARE REQUIRED, C THIS IS THE RECOMMENDED SETTING. C TOL = SQRT(SPMPAR(1)) C CALL LMSTR1(FCN,M,N,X,FVEC,FJAC,LDFJAC,TOL, * INFO,IPVT,WA,LWA) FNORM = ENORM(M,FVEC) WRITE (NWRITE,1000) FNORM,INFO,(X(J),J=1,N) STOP 1000 FORMAT (5X,31H FINAL L2 NORM OF THE RESIDUALS,E15.7 // * 5X,15H EXIT PARAMETER,16X,I10 // * 5X,27H FINAL APPROXIMATE SOLUTION // 5X,3E15.7) C 1 0 Page 0 C LAST CARD OF DRIVER FOR LMSTR1 EXAMPLE. C END SUBROUTINE FCN(M,N,X,FVEC,FJROW,IFLAG) INTEGER M,N,IFLAG REAL X(N),FVEC(M),FJROW(N) C C SUBROUTINE FCN FOR LMSTR1 EXAMPLE. C INTEGER I REAL TMP1,TMP2,TMP3,TMP4 REAL Y(15) DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8), * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15) * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1, * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/ C IF (IFLAG .GE. 2) GO TO 20 DO 10 I = 1, 15 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3)) 10 CONTINUE GO TO 40 20 CONTINUE I = IFLAG - 1 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2 FJROW(1) = -1.E0 FJROW(2) = TMP1*TMP2/TMP4 FJROW(3) = TMP1*TMP3/TMP4 30 CONTINUE 40 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 0 Results obtained with different compilers or machines may be slightly different. 0 FINAL L2 NORM OF THE RESIDUALS 0.9063596E-01 0 EXIT PARAMETER 1 0 FINAL APPROXIMATE SOLUTION 0 0.8241058E-01 0.1133037E+01 0.2343695E+01 1 0 1 0 Page 0 Documentation for MINPACK subroutine LMSTR 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of LMSTR is to minimize the sum of the squares of M nonlinear functions in N variables by a modification of the Levenberg-Marquardt algorithm which uses minimal storage. The user must provide a subroutine which calculates the functions and the rows of the Jacobian. 0 2. Subroutine and type statements. 0 SUBROUTINE LMSTR(FCN,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL, * MAXFEV,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,NJEV, * IPVT,QTF,WA1,WA2,WA3,WA4) INTEGER M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV INTEGER IPVT(N) REAL FTOL,XTOL,GTOL,FACTOR REAL X(N),FVEC(M),FJAC(LDFJAC,N),DIAG(N),QTF(N), * WA1(N),WA2(N),WA3(N),WA4(M) 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to LMSTR and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from LMSTR. 0 FCN is the name of the user-supplied subroutine which calculate the functions and the rows of the Jacobian. FCN must be declared in an EXTERNAL statement in the user calling program and should be written as follows. 0 SUBROUTINE FCN(M,N,X,FVEC,FJROW,IFLAG) INTEGER M,N,IFLAG REAL X(N),FVEC(M),FJROW(N) ---------- IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND RETURN THIS VECTOR IN FVEC. IF IFLAG = I CALCULATE THE (I-1)-ST ROW OF THE JACOBIAN AT X AND RETURN THIS VECTOR IN FJROW. ---------- RETURN 1 0 Page 0 END 0 The value of IFLAG should not be changed by FCN unless the user wants to terminate execution of LMSTR. In this case set IFLAG to a negative integer. 0 M is a positive integer input variable set to the number of functions. 0 N is a positive integer input variable set to the number of variables. N must not exceed M. 0 X is an array of length N. On input X must contain an initial estimate of the solution vector. On output X contains the final estimate of the solution vector. 0 FVEC is an output array of length M which contains the function evaluated at the output X. 0 FJAC is an output N by N array. The upper triangle of FJAC con tains an upper triangular matrix R such that 0 T T T P *(JAC *JAC)*P = R *R, 0 where P is a permutation matrix and JAC is the final calcu- lated Jacobian. Column j of P is column IPVT(j) (see below) of the identity matrix. The lower triangular part of FJAC contains information generated during the computation of R. 0 LDFJAC is a positive integer input variable not less than N which specifies the leading dimension of the array FJAC. 0 FTOL is a nonnegative input variable. Termination occurs when both the actual and predicted relative reductions in the sum of squares are at most FTOL. Therefore, FTOL measures the relative error desired in the sum of squares. Section 4 con- tains more details about FTOL. 0 XTOL is a nonnegative input variable. Termination occurs when the relative error between two consecutive iterates is at mos XTOL. Therefore, XTOL measures the relative error desired in the approximate solution. Section 4 contains more details about XTOL. 0 GTOL is a nonnegative input variable. Termination occurs when the cosine of the angle between FVEC and any column of the Jacobian is at most GTOL in absolute value. Therefore, GTOL measures the orthogonality desired between the function vecto and the columns of the Jacobian. Section 4 contains more details about GTOL. 0 MAXFEV is a positive integer input variable. Termination occur when the number of calls to FCN with IFLAG = 1 has reached 1 0 Page 0 MAXFEV. 0 DIAG is an array of length N. If MODE = 1 (see below), DIAG is internally set. If MODE = 2, DIAG must contain positive entries that serve as multiplicative scale factors for the variables. 0 MODE is an integer input variable. If MODE = 1, the variables will be scaled internally. If MODE = 2, the scaling is speci fied by the input DIAG. Other values of MODE are equivalent to MODE = 1. 0 FACTOR is a positive input variable used in determining the ini tial step bound. This bound is set to the product of FACTOR and the Euclidean norm of DIAG*X if nonzero, or else to FACTO itself. In most cases FACTOR should lie in the interval (.1,100.). 100. is a generally recommended value. 0 NPRINT is an integer input variable that enables controlled printing of iterates if it is positive. In this case, FCN is called with IFLAG = 0 at the beginning of the first iteration and every NPRINT iterations thereafter and immediately prior to return, with X and FVEC available for printing. If NPRINT is not positive, no special calls of FCN with IFLAG = 0 are made. 0 INFO is an integer output variable. If the user has terminated execution, INFO is set to the (negative) value of IFLAG. See description of FCN. Otherwise, INFO is set as follows. 0 INFO = 0 Improper input parameters. 0 INFO = 1 Both actual and predicted relative reductions in th sum of squares are at most FTOL. 0 INFO = 2 Relative error between two consecutive iterates is at most XTOL. 0 INFO = 3 Conditions for INFO = 1 and INFO = 2 both hold. 0 INFO = 4 The cosine of the angle between FVEC and any column of the Jacobian is at most GTOL in absolute value. 0 INFO = 5 Number of calls to FCN with IFLAG = 1 has reached MAXFEV. 0 INFO = 6 FTOL is too small. No further reduction in the sum of squares is possible. 0 INFO = 7 XTOL is too small. No further improvement in the approximate solution X is possible. 0 INFO = 8 GTOL is too small. FVEC is orthogonal to the columns of the Jacobian to machine precision. 1 0 Page 0 Sections 4 and 5 contain more details about INFO. 0 NFEV is an integer output variable set to the number of calls t FCN with IFLAG = 1. 0 NJEV is an integer output variable set to the number of calls t FCN with IFLAG = 2. 0 IPVT is an integer output array of length N. IPVT defines a permutation matrix P such that JAC*P = Q*R, where JAC is the final calculated Jacobian, Q is orthogonal (not stored), and is upper triangular. Column j of P is column IPVT(j) of the identity matrix. 0 QTF is an output array of length N which contains the first N elements of the vector (Q transpose)*FVEC. 0 WA1, WA2, and WA3 are work arrays of length N. 0 WA4 is a work array of length M. 0 4. Successful completion. 0 The accuracy of LMSTR is controlled by the convergence parame- ters FTOL, XTOL, and GTOL. These parameters are used in tests which make three types of comparisons between the approximation X and a solution XSOL. LMSTR terminates when any of the tests is satisfied. If any of the convergence parameters is less tha the machine precision (as defined by the MINPACK function SPMPAR(1)), then LMSTR only attempts to satisfy the test define by the machine precision. Further progress is not usually pos- sible. 0 The tests assume that the functions and the Jacobian are coded consistently, and that the functions are reasonably well behaved. If these conditions are not satisfied, then LMSTR may incorrectly indicate convergence. The coding of the Jacobian can be checked by the MINPACK subroutine CHKDER. If the Jaco- bian is coded correctly, then the validity of the answer can be checked, for example, by rerunning LMSTR with tighter toler- ances. 0 First convergence test. If ENORM(Z) denotes the Euclidean norm of a vector Z, then this test attempts to guarantee that 0 ENORM(FVEC) .LE. (1+FTOL)*ENORM(FVECS), 0 where FVECS denotes the functions evaluated at XSOL. If this condition is satisfied with FTOL = 10**(-K), then the final residual norm ENORM(FVEC) has K significant decimal digits an INFO is set to 1 (or to 3 if the second test is also satis- fied). Unless high precision solutions are required, the recommended value for FTOL is the square root of the machine 1 0 Page 0 precision. 0 Second convergence test. If D is the diagonal matrix whose entries are defined by the array DIAG, then this test attempt to guarantee that 0 ENORM(D*(X-XSOL)) .LE. XTOL*ENORM(D*XSOL). 0 If this condition is satisfied with XTOL = 10**(-K), then the larger components of D*X have K significant decimal digits an INFO is set to 2 (or to 3 if the first test is also satis- fied). There is a danger that the smaller components of D*X may have large relative errors, but if MODE = 1, then the accuracy of the components of X is usually related to their sensitivity. Unless high precision solutions are required, the recommended value for XTOL is the square root of the machine precision. 0 Third convergence test. This test is satisfied when the cosine of the angle between FVEC and any column of the Jacobian at X is at most GTOL in absolute value. There is no clear rela- tionship between this test and the accuracy of LMSTR, and furthermore, the test is equally well satisfied at other crit ical points, namely maximizers and saddle points. Therefore, termination caused by this test (INFO = 4) should be examined carefully. The recommended value for GTOL is zero. 0 5. Unsuccessful completion. 0 Unsuccessful termination of LMSTR can be due to improper input parameters, arithmetic interrupts, or an excessive number of function evaluations. 0 Improper input parameters. INFO is set to 0 if N .LE. 0, or M .LT. N, or LDFJAC .LT. N, or FTOL .LT. 0.E0, or XTOL .LT. 0.E0, or GTOL .LT. 0.E0, or MAXFEV .LE. 0, or FACTOR .LE. 0.E0. 0 Arithmetic interrupts. If these interrupts occur in the FCN subroutine during an early stage of the computation, they may be caused by an unacceptable choice of X by LMSTR. In this case, it may be possible to remedy the situation by rerunning LMSTR with a smaller value of FACTOR. 0 Excessive number of function evaluations. A reasonable value for MAXFEV is 100*(N+1). If the number of calls to FCN with IFLAG = 1 reaches MAXFEV, then this indicates that the routin is converging very slowly as measured by the progress of FVEC and INFO is set to 5. In this case, it may be helpful to restart LMSTR with MODE set to 1. 0 6. Characteristics of the algorithm. 1 0 Page 0 LMSTR is a modification of the Levenberg-Marquardt algorithm. Two of its main characteristics involve the proper use of implicitly scaled variables (if MODE = 1) and an optimal choice for the correction. The use of implicitly scaled variables achieves scale invariance of LMSTR and limits the size of the correction in any direction where the functions are changing rapidly. The optimal choice of the correction guarantees (unde reasonable conditions) global convergence from starting points far from the solution and a fast rate of convergence for prob- lems with small residuals. 0 Timing. The time required by LMSTR to solve a given problem depends on M and N, the behavior of the functions, the accu- racy requested, and the starting point. The number of arith- metic operations needed by LMSTR is about N**3 to process eac evaluation of the functions (call to FCN with IFLAG = 1) and 1.5*(N**2) to process each row of the Jacobian (call to FCN with IFLAG .GE. 2). Unless FCN can be evaluated quickly, the timing of LMSTR will be strongly influenced by the time spent in FCN. 0 Storage. LMSTR requires N**2 + 2*M + 6*N single precision sto- rage locations and N integer storage locations, in addition t the storage required by the program. There are no internally declared storage arrays. 0 7. Subprograms required. 0 USER-supplied ...... FCN 0 MINPACK-supplied ... SPMPAR,ENORM,LMPAR,QRFAC,QRSOLV,RWUPDT 0 FORTRAN-supplied ... ABS,AMAX1,AMIN1,SQRT,MOD 0 8. References. 0 Jorge J. More, The Levenberg-Marquardt Algorithm, Implementatio and Theory. Numerical Analysis, G. A. Watson, editor. Lecture Notes in Mathematics 630, Springer-Verlag, 1977. 0 9. Example. 0 The problem is to determine the values of x(1), x(2), and x(3) which provide the best fit (in the least squares sense) of 0 x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)), i = 1, 15 0 to the data 0 y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39, 0.37,0.58,0.73,0.96,1.34,2.10,4.39), 1 0 Page 0 where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)). The i-th component of FVEC is thus defined by 0 y(i) - (x(1) + u(i)/(v(i)*x(2) + w(i)*x(3))). 0 C ********** C C DRIVER FOR LMSTR EXAMPLE. C SINGLE PRECISION VERSION C C ********** INTEGER J,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,NWRITE INTEGER IPVT(3) REAL FTOL,XTOL,GTOL,FACTOR,FNORM REAL X(3),FVEC(15),FJAC(3,3),DIAG(3),QTF(3), * WA1(3),WA2(3),WA3(3),WA4(15) REAL ENORM,SPMPAR EXTERNAL FCN C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C M = 15 N = 3 C C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH FIT. C X(1) = 1.E0 X(2) = 1.E0 X(3) = 1.E0 C LDFJAC = 3 C C SET FTOL AND XTOL TO THE SQUARE ROOT OF THE MACHINE PRECISION C AND GTOL TO ZERO. UNLESS HIGH PRECISION SOLUTIONS ARE C REQUIRED, THESE ARE THE RECOMMENDED SETTINGS. C FTOL = SQRT(SPMPAR(1)) XTOL = SQRT(SPMPAR(1)) GTOL = 0.E0 C MAXFEV = 400 MODE = 1 FACTOR = 1.E2 NPRINT = 0 C CALL LMSTR(FCN,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,GTOL, * MAXFEV,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,NJEV, * IPVT,QTF,WA1,WA2,WA3,WA4) FNORM = ENORM(M,FVEC) WRITE (NWRITE,1000) FNORM,NFEV,NJEV,INFO,(X(J),J=1,N) STOP 1000 FORMAT (5X,31H FINAL L2 NORM OF THE RESIDUALS,E15.7 // 1 0 Page 0 * 5X,31H NUMBER OF FUNCTION EVALUATIONS,I10 // * 5X,31H NUMBER OF JACOBIAN EVALUATIONS,I10 // * 5X,15H EXIT PARAMETER,16X,I10 // * 5X,27H FINAL APPROXIMATE SOLUTION // 5X,3E15.7) C C LAST CARD OF DRIVER FOR LMSTR EXAMPLE. C END SUBROUTINE FCN(M,N,X,FVEC,FJROW,IFLAG) INTEGER M,N,IFLAG REAL X(N),FVEC(M),FJROW(N) C C SUBROUTINE FCN FOR LMSTR EXAMPLE. C INTEGER I REAL TMP1,TMP2,TMP3,TMP4 REAL Y(15) DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8), * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15) * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1, * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/ C IF (IFLAG .NE. 0) GO TO 5 C C INSERT PRINT STATEMENTS HERE WHEN NPRINT IS POSITIVE. C RETURN 5 CONTINUE IF (IFLAG .GE. 2) GO TO 20 DO 10 I = 1, 15 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3)) 10 CONTINUE GO TO 40 20 CONTINUE I = IFLAG - 1 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2 FJROW(1) = -1.E0 FJROW(2) = TMP1*TMP2/TMP4 FJROW(3) = TMP1*TMP3/TMP4 30 CONTINUE 40 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 1 0 Page 0 Results obtained with different compilers or machines may be slightly different. 0 FINAL L2 NORM OF THE RESIDUALS 0.9063596E-01 0 NUMBER OF FUNCTION EVALUATIONS 6 0 NUMBER OF JACOBIAN EVALUATIONS 5 0 EXIT PARAMETER 1 0 FINAL APPROXIMATE SOLUTION 0 0.8241058E-01 0.1133037E+01 0.2343695E+01 1 0 1 0 Page 0 Documentation for MINPACK subroutine LMDIF1 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of LMDIF1 is to minimize the sum of the squares of nonlinear functions in N variables by a modification of the Levenberg-Marquardt algorithm. This is done by using the more general least-squares solver LMDIF. The user must provide a subroutine which calculates the functions. The Jacobian is the calculated by a forward-difference approximation. 0 2. Subroutine and type statements. 0 SUBROUTINE LMDIF1(FCN,M,N,X,FVEC,TOL,INFO,IWA,WA,LWA) INTEGER M,N,INFO,LWA INTEGER IWA(N) REAL TOL REAL X(N),FVEC(M),WA(LWA) EXTERNAL FCN 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to LMDIF1 and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from LMDIF1. 0 FCN is the name of the user-supplied subroutine which calculate the functions. FCN must be declared in an EXTERNAL statement in the user calling program, and should be written as follows 0 SUBROUTINE FCN(M,N,X,FVEC,IFLAG) INTEGER M,N,IFLAG REAL X(N),FVEC(M) ---------- CALCULATE THE FUNCTIONS AT X AND RETURN THIS VECTOR IN FVEC. ---------- RETURN END 0 The value of IFLAG should not be changed by FCN unless the user wants to terminate execution of LMDIF1. In this case se 1 0 Page 0 IFLAG to a negative integer. 0 M is a positive integer input variable set to the number of functions. 0 N is a positive integer input variable set to the number of variables. N must not exceed M. 0 X is an array of length N. On input X must contain an initial estimate of the solution vector. On output X contains the final estimate of the solution vector. 0 FVEC is an output array of length M which contains the function evaluated at the output X. 0 TOL is a nonnegative input variable. Termination occurs when the algorithm estimates either that the relative error in the sum of squares is at most TOL or that the relative error between X and the solution is at most TOL. Section 4 contain more details about TOL. 0 INFO is an integer output variable. If the user has terminated execution, INFO is set to the (negative) value of IFLAG. See description of FCN. Otherwise, INFO is set as follows. 0 INFO = 0 Improper input parameters. 0 INFO = 1 Algorithm estimates that the relative error in the sum of squares is at most TOL. 0 INFO = 2 Algorithm estimates that the relative error between X and the solution is at most TOL. 0 INFO = 3 Conditions for INFO = 1 and INFO = 2 both hold. 0 INFO = 4 FVEC is orthogonal to the columns of the Jacobian t machine precision. 0 INFO = 5 Number of calls to FCN has reached or exceeded 200*(N+1). 0 INFO = 6 TOL is too small. No further reduction in the sum of squares is possible. 0 INFO = 7 TOL is too small. No further improvement in the approximate solution X is possible. 0 Sections 4 and 5 contain more details about INFO. 0 IWA is an integer work array of length N. 0 WA is a work array of length LWA. 0 LWA is a positive integer input variable not less than 1 0 Page 0 M*N+5*N+M. 0 4. Successful completion. 0 The accuracy of LMDIF1 is controlled by the convergence parame- ter TOL. This parameter is used in tests which make three type of comparisons between the approximation X and a solution XSOL. LMDIF1 terminates when any of the tests is satisfied. If TOL i less than the machine precision (as defined by the MINPACK func tion SPMPAR(1)), then LMDIF1 only attempts to satisfy the test defined by the machine precision. Further progress is not usu- ally possible. Unless high precision solutions are required, the recommended value for TOL is the square root of the machine precision. 0 The tests assume that the functions are reasonably well behaved If this condition is not satisfied, then LMDIF1 may incorrectly indicate convergence. The validity of the answer can be checked, for example, by rerunning LMDIF1 with a tighter toler- ance. 0 First convergence test. If ENORM(Z) denotes the Euclidean norm of a vector Z, then this test attempts to guarantee that 0 ENORM(FVEC) .LE. (1+TOL)*ENORM(FVECS), 0 where FVECS denotes the functions evaluated at XSOL. If this condition is satisfied with TOL = 10**(-K), then the final residual norm ENORM(FVEC) has K significant decimal digits an INFO is set to 1 (or to 3 if the second test is also satis- fied). 0 Second convergence test. If D is a diagonal matrix (implicitly generated by LMDIF1) whose entries contain scale factors for the variables, then this test attempts to guarantee that 0 ENORM(D*(X-XSOL)) .LE. TOL*ENORM(D*XSOL). 0 If this condition is satisfied with TOL = 10**(-K), then the larger components of D*X have K significant decimal digits an INFO is set to 2 (or to 3 if the first test is also satis- fied). There is a danger that the smaller components of D*X may have large relative errors, but the choice of D is such that the accuracy of the components of X is usually related t their sensitivity. 0 Third convergence test. This test is satisfied when FVEC is orthogonal to the columns of the Jacobian to machine preci- sion. There is no clear relationship between this test and the accuracy of LMDIF1, and furthermore, the test is equally well satisfied at other critical points, namely maximizers an saddle points. Also, errors in the functions (see below) may result in the test being satisfied at a point not close to th 1 0 Page 0 minimum. Therefore, termination caused by this test (INFO = 4) should be examined carefully. 0 5. Unsuccessful completion. 0 Unsuccessful termination of LMDIF1 can be due to improper input parameters, arithmetic interrupts, an excessive number of func- tion evaluations, or errors in the functions. 0 Improper input parameters. INFO is set to 0 if N .LE. 0, or M .LT. N, or TOL .LT. 0.E0, or LWA .LT. M*N+5*N+M. 0 Arithmetic interrupts. If these interrupts occur in the FCN subroutine during an early stage of the computation, they may be caused by an unacceptable choice of X by LMDIF1. In this case, it may be possible to remedy the situation by not evalu ating the functions here, but instead setting the components of FVEC to numbers that exceed those in the initial FVEC, thereby indirectly reducing the step length. The step length can be more directly controlled by using instead LMDIF, which includes in its calling sequence the step-length-governing parameter FACTOR. 0 Excessive number of function evaluations. If the number of calls to FCN reaches 200*(N+1), then this indicates that the routine is converging very slowly as measured by the progress of FVEC, and INFO is set to 5. In this case, it may be help- ful to restart LMDIF1, thereby forcing it to disregard old (and possibly harmful) information. 0 Errors in the functions. The choice of step length in the for- ward-difference approximation to the Jacobian assumes that th relative errors in the functions are of the order of the machine precision. If this is not the case, LMDIF1 may fail (usually with INFO = 4). The user should then use LMDIF instead, or one of the programs which require the analytic Jacobian (LMDER1 and LMDER). 0 6. Characteristics of the algorithm. 0 LMDIF1 is a modification of the Levenberg-Marquardt algorithm. Two of its main characteristics involve the proper use of implicitly scaled variables and an optimal choice for the cor- rection. The use of implicitly scaled variables achieves scale invariance of LMDIF1 and limits the size of the correction in any direction where the functions are changing rapidly. The optimal choice of the correction guarantees (under reasonable conditions) global convergence from starting points far from th solution and a fast rate of convergence for problems with small residuals. 0 Timing. The time required by LMDIF1 to solve a given problem 1 0 Page 0 depends on M and N, the behavior of the functions, the accu- racy requested, and the starting point. The number of arith- metic operations needed by LMDIF1 is about N**3 to process each evaluation of the functions (one call to FCN) and M*(N**2) to process each approximation to the Jacobian (N calls to FCN). Unless FCN can be evaluated quickly, the tim- ing of LMDIF1 will be strongly influenced by the time spent i FCN. 0 Storage. LMDIF1 requires M*N + 2*M + 6*N single precision sto- rage locations and N integer storage locations, in addition t the storage required by the program. There are no internally declared storage arrays. 0 7. Subprograms required. 0 USER-supplied ...... FCN 0 MINPACK-supplied ... SPMPAR,ENORM,FDJAC2,LMDIF,LMPAR, QRFAC,QRSOLV 0 FORTRAN-supplied ... ABS,AMAX1,AMIN1,SQRT,MOD 0 8. References. 0 Jorge J. More, The Levenberg-Marquardt Algorithm, Implementatio and Theory. Numerical Analysis, G. A. Watson, editor. Lecture Notes in Mathematics 630, Springer-Verlag, 1977. 0 9. Example. 0 The problem is to determine the values of x(1), x(2), and x(3) which provide the best fit (in the least squares sense) of 0 x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)), i = 1, 15 0 to the data 0 y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39, 0.37,0.58,0.73,0.96,1.34,2.10,4.39), 0 where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)). The i-th component of FVEC is thus defined by 0 y(i) - (x(1) + u(i)/(v(i)*x(2) + w(i)*x(3))). 0 C ********** C C DRIVER FOR LMDIF1 EXAMPLE. C SINGLE PRECISION VERSION C 1 0 Page 0 C ********** INTEGER J,M,N,INFO,LWA,NWRITE INTEGER IWA(3) REAL TOL,FNORM REAL X(3),FVEC(15),WA(75) REAL ENORM,SPMPAR EXTERNAL FCN C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C M = 15 N = 3 C C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH FIT. C X(1) = 1.E0 X(2) = 1.E0 X(3) = 1.E0 C LWA = 75 C C SET TOL TO THE SQUARE ROOT OF THE MACHINE PRECISION. C UNLESS HIGH PRECISION SOLUTIONS ARE REQUIRED, C THIS IS THE RECOMMENDED SETTING. C TOL = SQRT(SPMPAR(1)) C CALL LMDIF1(FCN,M,N,X,FVEC,TOL,INFO,IWA,WA,LWA) FNORM = ENORM(M,FVEC) WRITE (NWRITE,1000) FNORM,INFO,(X(J),J=1,N) STOP 1000 FORMAT (5X,31H FINAL L2 NORM OF THE RESIDUALS,E15.7 // * 5X,15H EXIT PARAMETER,16X,I10 // * 5X,27H FINAL APPROXIMATE SOLUTION // 5X,3E15.7) C C LAST CARD OF DRIVER FOR LMDIF1 EXAMPLE. C END SUBROUTINE FCN(M,N,X,FVEC,IFLAG) INTEGER M,N,IFLAG REAL X(N),FVEC(M) C C SUBROUTINE FCN FOR LMDIF1 EXAMPLE. C INTEGER I REAL TMP1,TMP2,TMP3 REAL Y(15) DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8), * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15) * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1, * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/ C 1 0 Page 0 DO 10 I = 1, 15 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3)) 10 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 0 Results obtained with different compilers or machines may be slightly different. 0 FINAL L2 NORM OF THE RESIDUALS 0.9063596E-01 0 EXIT PARAMETER 1 0 FINAL APPROXIMATE SOLUTION 0 0.8241057E-01 0.1133037E+01 0.2343695E+01 1 0 1 0 Page 0 Documentation for MINPACK subroutine LMDIF 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of LMDIF is to minimize the sum of the squares of M nonlinear functions in N variables by a modification of the Levenberg-Marquardt algorithm. The user must provide a subrou- tine which calculates the functions. The Jacobian is then cal- culated by a forward-difference approximation. 0 2. Subroutine and type statements. 0 SUBROUTINE LMDIF(FCN,M,N,X,FVEC,FTOL,XTOL,GTOL,MAXFEV,EPSFCN, * DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,FJAC,LDFJAC, * IPVT,QTF,WA1,WA2,WA3,WA4) INTEGER M,N,MAXFEV,MODE,NPRINT,INFO,NFEV,LDFJAC INTEGER IPVT(N) REAL FTOL,XTOL,GTOL,EPSFCN,FACTOR REAL X(N),FVEC(M),DIAG(N),FJAC(LDFJAC,N),QTF(N), * WA1(N),WA2(N),WA3(N),WA4(M) EXTERNAL FCN 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to LMDIF and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from LMDIF. 0 FCN is the name of the user-supplied subroutine which calculate the functions. FCN must be declared in an EXTERNAL statement in the user calling program, and should be written as follows 0 SUBROUTINE FCN(M,N,X,FVEC,IFLAG) INTEGER M,N,IFLAG REAL X(N),FVEC(M) ---------- CALCULATE THE FUNCTIONS AT X AND RETURN THIS VECTOR IN FVEC. ---------- RETURN END 1 0 Page 0 The value of IFLAG should not be changed by FCN unless the user wants to terminate execution of LMDIF. In this case set IFLAG to a negative integer. 0 M is a positive integer input variable set to the number of functions. 0 N is a positive integer input variable set to the number of variables. N must not exceed M. 0 X is an array of length N. On input X must contain an initial estimate of the solution vector. On output X contains the final estimate of the solution vector. 0 FVEC is an output array of length M which contains the function evaluated at the output X. 0 FTOL is a nonnegative input variable. Termination occurs when both the actual and predicted relative reductions in the sum of squares are at most FTOL. Therefore, FTOL measures the relative error desired in the sum of squares. Section 4 con- tains more details about FTOL. 0 XTOL is a nonnegative input variable. Termination occurs when the relative error between two consecutive iterates is at mos XTOL. Therefore, XTOL measures the relative error desired in the approximate solution. Section 4 contains more details about XTOL. 0 GTOL is a nonnegative input variable. Termination occurs when the cosine of the angle between FVEC and any column of the Jacobian is at most GTOL in absolute value. Therefore, GTOL measures the orthogonality desired between the function vecto and the columns of the Jacobian. Section 4 contains more details about GTOL. 0 MAXFEV is a positive integer input variable. Termination occur when the number of calls to FCN is at least MAXFEV by the end of an iteration. 0 EPSFCN is an input variable used in determining a suitable step for the forward-difference approximation. This approximation assumes that the relative errors in the functions are of the order of EPSFCN. If EPSFCN is less than the machine preci- sion, it is assumed that the relative errors in the functions are of the order of the machine precision. 0 DIAG is an array of length N. If MODE = 1 (see below), DIAG is internally set. If MODE = 2, DIAG must contain positive entries that serve as multiplicative scale factors for the variables. 0 MODE is an integer input variable. If MODE = 1, the variables will be scaled internally. If MODE = 2, the scaling is 1 0 Page 0 specified by the input DIAG. Other values of MODE are equiva lent to MODE = 1. 0 FACTOR is a positive input variable used in determining the ini tial step bound. This bound is set to the product of FACTOR and the Euclidean norm of DIAG*X if nonzero, or else to FACTO itself. In most cases FACTOR should lie in the interval (.1,100.). 100. is a generally recommended value. 0 NPRINT is an integer input variable that enables controlled printing of iterates if it is positive. In this case, FCN is called with IFLAG = 0 at the beginning of the first iteration and every NPRINT iterations thereafter and immediately prior to return, with X and FVEC available for printing. If NPRINT is not positive, no special calls of FCN with IFLAG = 0 are made. 0 INFO is an integer output variable. If the user has terminated execution, INFO is set to the (negative) value of IFLAG. See description of FCN. Otherwise, INFO is set as follows. 0 INFO = 0 Improper input parameters. 0 INFO = 1 Both actual and predicted relative reductions in th sum of squares are at most FTOL. 0 INFO = 2 Relative error between two consecutive iterates is at most XTOL. 0 INFO = 3 Conditions for INFO = 1 and INFO = 2 both hold. 0 INFO = 4 The cosine of the angle between FVEC and any column of the Jacobian is at most GTOL in absolute value. 0 INFO = 5 Number of calls to FCN has reached or exceeded MAXFEV. 0 INFO = 6 FTOL is too small. No further reduction in the sum of squares is possible. 0 INFO = 7 XTOL is too small. No further improvement in the approximate solution X is possible. 0 INFO = 8 GTOL is too small. FVEC is orthogonal to the columns of the Jacobian to machine precision. 0 Sections 4 and 5 contain more details about INFO. 0 NFEV is an integer output variable set to the number of calls t FCN. 0 FJAC is an output M by N array. The upper N by N submatrix of FJAC contains an upper triangular matrix R with diagonal ele- ments of nonincreasing magnitude such that 1 0 Page 0 T T T P *(JAC *JAC)*P = R *R, 0 where P is a permutation matrix and JAC is the final calcu- lated Jacobian. Column j of P is column IPVT(j) (see below) of the identity matrix. The lower trapezoidal part of FJAC contains information generated during the computation of R. 0 LDFJAC is a positive integer input variable not less than M which specifies the leading dimension of the array FJAC. 0 IPVT is an integer output array of length N. IPVT defines a permutation matrix P such that JAC*P = Q*R, where JAC is the final calculated Jacobian, Q is orthogonal (not stored), and is upper triangular with diagonal elements of nonincreasing magnitude. Column j of P is column IPVT(j) of the identity matrix. 0 QTF is an output array of length N which contains the first N elements of the vector (Q transpose)*FVEC. 0 WA1, WA2, and WA3 are work arrays of length N. 0 WA4 is a work array of length M. 0 4. Successful completion. 0 The accuracy of LMDIF is controlled by the convergence parame- ters FTOL, XTOL, and GTOL. These parameters are used in tests which make three types of comparisons between the approximation X and a solution XSOL. LMDIF terminates when any of the tests is satisfied. If any of the convergence parameters is less tha the machine precision (as defined by the MINPACK function SPMPAR(1)), then LMDIF only attempts to satisfy the test define by the machine precision. Further progress is not usually pos- sible. 0 The tests assume that the functions are reasonably well behaved If this condition is not satisfied, then LMDIF may incorrectly indicate convergence. The validity of the answer can be checked, for example, by rerunning LMDIF with tighter toler- ances. 0 First convergence test. If ENORM(Z) denotes the Euclidean norm of a vector Z, then this test attempts to guarantee that 0 ENORM(FVEC) .LE. (1+FTOL)*ENORM(FVECS), 0 where FVECS denotes the functions evaluated at XSOL. If this condition is satisfied with FTOL = 10**(-K), then the final residual norm ENORM(FVEC) has K significant decimal digits an INFO is set to 1 (or to 3 if the second test is also satis- fied). Unless high precision solutions are required, the 1 0 Page 0 recommended value for FTOL is the square root of the machine precision. 0 Second convergence test. If D is the diagonal matrix whose entries are defined by the array DIAG, then this test attempt to guarantee that 0 ENORM(D*(X-XSOL)) .LE. XTOL*ENORM(D*XSOL). 0 If this condition is satisfied with XTOL = 10**(-K), then the larger components of D*X have K significant decimal digits an INFO is set to 2 (or to 3 if the first test is also satis- fied). There is a danger that the smaller components of D*X may have large relative errors, but if MODE = 1, then the accuracy of the components of X is usually related to their sensitivity. Unless high precision solutions are required, the recommended value for XTOL is the square root of the machine precision. 0 Third convergence test. This test is satisfied when the cosine of the angle between FVEC and any column of the Jacobian at X is at most GTOL in absolute value. There is no clear rela- tionship between this test and the accuracy of LMDIF, and furthermore, the test is equally well satisfied at other crit ical points, namely maximizers and saddle points. Therefore, termination caused by this test (INFO = 4) should be examined carefully. The recommended value for GTOL is zero. 0 5. Unsuccessful completion. 0 Unsuccessful termination of LMDIF can be due to improper input parameters, arithmetic interrupts, or an excessive number of function evaluations. 0 Improper input parameters. INFO is set to 0 if N .LE. 0, or M .LT. N, or LDFJAC .LT. M, or FTOL .LT. 0.E0, or XTOL .LT. 0.E0, or GTOL .LT. 0.E0, or MAXFEV .LE. 0, or FACTOR .LE. 0.E0. 0 Arithmetic interrupts. If these interrupts occur in the FCN subroutine during an early stage of the computation, they may be caused by an unacceptable choice of X by LMDIF. In this case, it may be possible to remedy the situation by rerunning LMDIF with a smaller value of FACTOR. 0 Excessive number of function evaluations. A reasonable value for MAXFEV is 200*(N+1). If the number of calls to FCN reaches MAXFEV, then this indicates that the routine is con- verging very slowly as measured by the progress of FVEC, and INFO is set to 5. In this case, it may be helpful to restart LMDIF with MODE set to 1. 0 1 0 Page 0 6. Characteristics of the algorithm. 0 LMDIF is a modification of the Levenberg-Marquardt algorithm. Two of its main characteristics involve the proper use of implicitly scaled variables (if MODE = 1) and an optimal choice for the correction. The use of implicitly scaled variables achieves scale invariance of LMDIF and limits the size of the correction in any direction where the functions are changing rapidly. The optimal choice of the correction guarantees (unde reasonable conditions) global convergence from starting points far from the solution and a fast rate of convergence for prob- lems with small residuals. 0 Timing. The time required by LMDIF to solve a given problem depends on M and N, the behavior of the functions, the accu- racy requested, and the starting point. The number of arith- metic operations needed by LMDIF is about N**3 to process eac evaluation of the functions (one call to FCN) and M*(N**2) to process each approximation to the Jacobian (N calls to FCN). Unless FCN can be evaluated quickly, the timing of LMDIF will be strongly influenced by the time spent in FCN. 0 Storage. LMDIF requires M*N + 2*M + 6*N single precision sto- rage locations and N integer storage locations, in addition t the storage required by the program. There are no internally declared storage arrays. 0 7. Subprograms required. 0 USER-supplied ...... FCN 0 MINPACK-supplied ... SPMPAR,ENORM,FDJAC2,LMPAR,QRFAC,QRSOLV 0 FORTRAN-supplied ... ABS,AMAX1,AMIN1,SQRT,MOD 0 8. References. 0 Jorge J. More, The Levenberg-Marquardt Algorithm, Implementatio and Theory. Numerical Analysis, G. A. Watson, editor. Lecture Notes in Mathematics 630, Springer-Verlag, 1977. 0 9. Example. 0 The problem is to determine the values of x(1), x(2), and x(3) which provide the best fit (in the least squares sense) of 0 x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)), i = 1, 15 0 to the data 1 0 Page 0 y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39, 0.37,0.58,0.73,0.96,1.34,2.10,4.39), 0 where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)). The i-th component of FVEC is thus defined by 0 y(i) - (x(1) + u(i)/(v(i)*x(2) + w(i)*x(3))). 0 C ********** C C DRIVER FOR LMDIF EXAMPLE. C SINGLE PRECISION VERSION C C ********** INTEGER J,M,N,MAXFEV,MODE,NPRINT,INFO,NFEV,LDFJAC,NWRITE INTEGER IPVT(3) REAL FTOL,XTOL,GTOL,EPSFCN,FACTOR,FNORM REAL X(3),FVEC(15),DIAG(3),FJAC(15,3),QTF(3), * WA1(3),WA2(3),WA3(3),WA4(15) REAL ENORM,SPMPAR EXTERNAL FCN C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C M = 15 N = 3 C C THE FOLLOWING STARTING VALUES PROVIDE A ROUGH FIT. C X(1) = 1.E0 X(2) = 1.E0 X(3) = 1.E0 C LDFJAC = 15 C C SET FTOL AND XTOL TO THE SQUARE ROOT OF THE MACHINE PRECISION C AND GTOL TO ZERO. UNLESS HIGH PRECISION SOLUTIONS ARE C REQUIRED, THESE ARE THE RECOMMENDED SETTINGS. C FTOL = SQRT(SPMPAR(1)) XTOL = SQRT(SPMPAR(1)) GTOL = 0.E0 C MAXFEV = 800 EPSFCN = 0.E0 MODE = 1 FACTOR = 1.E2 NPRINT = 0 C CALL LMDIF(FCN,M,N,X,FVEC,FTOL,XTOL,GTOL,MAXFEV,EPSFCN, * DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,FJAC,LDFJAC, * IPVT,QTF,WA1,WA2,WA3,WA4) 1 0 Page 0 FNORM = ENORM(M,FVEC) WRITE (NWRITE,1000) FNORM,NFEV,INFO,(X(J),J=1,N) STOP 1000 FORMAT (5X,31H FINAL L2 NORM OF THE RESIDUALS,E15.7 // * 5X,31H NUMBER OF FUNCTION EVALUATIONS,I10 // * 5X,15H EXIT PARAMETER,16X,I10 // * 5X,27H FINAL APPROXIMATE SOLUTION // 5X,3E15.7) C C LAST CARD OF DRIVER FOR LMDIF EXAMPLE. C END SUBROUTINE FCN(M,N,X,FVEC,IFLAG) INTEGER M,N,IFLAG REAL X(N),FVEC(M) C C SUBROUTINE FCN FOR LMDIF EXAMPLE. C INTEGER I REAL TMP1,TMP2,TMP3 REAL Y(15) DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8), * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15) * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1, * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/ C IF (IFLAG .NE. 0) GO TO 5 C C INSERT PRINT STATEMENTS HERE WHEN NPRINT IS POSITIVE. C RETURN 5 CONTINUE DO 10 I = 1, 15 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3)) 10 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 0 Results obtained with different compilers or machines may be slightly different. 0 FINAL L2 NORM OF THE RESIDUALS 0.9063596E-01 0 NUMBER OF FUNCTION EVALUATIONS 21 0 EXIT PARAMETER 1 0 FINAL APPROXIMATE SOLUTION 1 0 Page 0 0.8241057E-01 0.1133037E+01 0.2343695E+01 1 0 1 0 Page 0 Documentation for MINPACK subroutine CHKDER 0 Single precision version 0 Argonne National Laboratory 0 Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More 0 March 1980 0 1. Purpose. 0 The purpose of CHKDER is to check the gradients of M nonlinear functions in N variables, evaluated at a point X, for consis- tency with the functions themselves. The user must call CHKDER twice, first with MODE = 1 and then with MODE = 2. 0 2. Subroutine and type statements. 0 SUBROUTINE CHKDER(M,N,X,FVEC,FJAC,LDFJAC,XP,FVECP,MODE,ERR) INTEGER M,N,LDFJAC,MODE REAL X(N),FVEC(M),FJAC(LDFJAC,N),XP(N),FVECP(M),ERR(M) 0 3. Parameters. 0 Parameters designated as input parameters must be specified on entry to CHKDER and are not changed on exit, while parameters designated as output parameters need not be specified on entry and are set to appropriate values on exit from CHKDER. 0 M is a positive integer input variable set to the number of functions. 0 N is a positive integer input variable set to the number of variables. 0 X is an input array of length N. 0 FVEC is an array of length M. On input when MODE = 2, FVEC mus contain the functions evaluated at X. 0 FJAC is an M by N array. On input when MODE = 2, the rows of FJAC must contain the gradients of the respective functions evaluated at X. 0 LDFJAC is a positive integer input variable not less than M which specifies the leading dimension of the array FJAC. 0 XP is an array of length N. On output when MODE = 1, XP is set to a neighboring point of X. 1 0 Page 0 FVECP is an array of length M. On input when MODE = 2, FVECP must contain the functions evaluated at XP. 0 MODE is an integer input variable set to 1 on the first call an 2 on the second. Other values of MODE are equivalent to MODE = 1. 0 ERR is an array of length M. On output when MODE = 2, ERR con- tains measures of correctness of the respective gradients. I there is no severe loss of significance, then if ERR(I) is 1. the I-th gradient is correct, while if ERR(I) is 0.0 the I-th gradient is incorrect. For values of ERR between 0.0 and 1.0 the categorization is less certain. In general, a value of ERR(I) greater than 0.5 indicates that the I-th gradient is probably correct, while a value of ERR(I) less than 0.5 indi- cates that the I-th gradient is probably incorrect. 0 4. Successful completion. 0 CHKDER usually guarantees that if ERR(I) is 1.0, then the I-th gradient at X is consistent with the I-th function. This sug- gests that the input X be such that consistency of the gradient at X implies consistency of the gradient at all points of inter est. If all the components of X are distinct and the fractiona part of each one has two nonzero digits, then X is likely to be a satisfactory choice. 0 If ERR(I) is not 1.0 but is greater than 0.5, then the I-th gra dient is probably consistent with the I-th function (the more s the larger ERR(I) is), but the conditions for ERR(I) to be 1.0 have not been completely satisfied. In this case, it is recom- mended that CHKDER be rerun with other input values of X. If ERR(I) is always greater than 0.5, then the I-th gradient is consistent with the I-th function. 0 5. Unsuccessful completion. 0 CHKDER does not perform reliably if cancellation or rounding errors cause a severe loss of significance in the evaluation of a function. Therefore, none of the components of X should be unusually small (in particular, zero) or any other value which may cause loss of significance. The relative differences between corresponding elements of FVECP and FVEC should be at least two orders of magnitude greater than the machine precisio (as defined by the MINPACK function SPMPAR(1)). If there is a severe loss of significance in the evaluation of the I-th func- tion, then ERR(I) may be 0.0 and yet the I-th gradient could be correct. 0 If ERR(I) is not 0.0 but is less than 0.5, then the I-th gra- dient is probably not consistent with the I-th function (the more so the smaller ERR(I) is), but the conditions for ERR(I) t 1 0 Page 0 be 0.0 have not been completely satisfied. In this case, it is recommended that CHKDER be rerun with other input values of X. If ERR(I) is always less than 0.5 and if there is no severe los of significance, then the I-th gradient is not consistent with the I-th function. 0 6. Characteristics of the algorithm. 0 CHKDER checks the I-th gradient for consistency with the I-th function by computing a forward-difference approximation along suitably chosen direction and comparing this approximation with the user-supplied gradient along the same direction. The prin- cipal characteristic of CHKDER is its invariance to changes in scale of the variables or functions. 0 Timing. The time required by CHKDER depends only on M and N. The number of arithmetic operations needed by CHKDER is about N when MODE = 1 and M*N when MODE = 2. 0 Storage. CHKDER requires M*N + 3*M + 2*N single precision stor age locations, in addition to the storage required by the pro gram. There are no internally declared storage arrays. 0 7. Subprograms required. 0 MINPACK-supplied ... SPMPAR 0 FORTRAN-supplied ... ABS,ALOG10,SQRT 0 8. References. 0 None. 0 9. Example. 0 This example checks the Jacobian matrix for the problem that determines the values of x(1), x(2), and x(3) which provide the best fit (in the least squares sense) of 0 x(1) + u(i)/(v(i)*x(2) + w(i)*x(3)), i = 1, 15 0 to the data 0 y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39, 0.37,0.58,0.73,0.96,1.34,2.10,4.39), 0 where u(i) = i, v(i) = 16 - i, and w(i) = min(u(i),v(i)). The i-th component of FVEC is thus defined by 0 y(i) - (x(1) + u(i)/(v(i)*x(2) + w(i)*x(3))). 1 0 Page 0 C ********** C C DRIVER FOR CHKDER EXAMPLE. C SINGLE PRECISION VERSION C C ********** INTEGER I,M,N,LDFJAC,MODE,NWRITE REAL X(3),FVEC(15),FJAC(15,3),XP(3),FVECP(15),ERR(15) C C LOGICAL OUTPUT UNIT IS ASSUMED TO BE NUMBER 6. C DATA NWRITE /6/ C M = 15 N = 3 C C THE FOLLOWING VALUES SHOULD BE SUITABLE FOR C CHECKING THE JACOBIAN MATRIX. C X(1) = 9.2E-1 X(2) = 1.3E-1 X(3) = 5.4E-1 C LDFJAC = 15 C MODE = 1 CALL CHKDER(M,N,X,FVEC,FJAC,LDFJAC,XP,FVECP,MODE,ERR) MODE = 2 CALL FCN(M,N,X,FVEC,FJAC,LDFJAC,1) CALL FCN(M,N,X,FVEC,FJAC,LDFJAC,2) CALL FCN(M,N,XP,FVECP,FJAC,LDFJAC,1) CALL CHKDER(M,N,X,FVEC,FJAC,LDFJAC,XP,FVECP,MODE,ERR) C DO 10 I = 1, M FVECP(I) = FVECP(I) - FVEC(I) 10 CONTINUE WRITE (NWRITE,1000) (FVEC(I),I=1,M) WRITE (NWRITE,2000) (FVECP(I),I=1,M) WRITE (NWRITE,3000) (ERR(I),I=1,M) STOP 1000 FORMAT (/5X,5H FVEC // (5X,3E15.7)) 2000 FORMAT (/5X,13H FVECP - FVEC // (5X,3E15.7)) 3000 FORMAT (/5X,4H ERR // (5X,3E15.7)) C C LAST CARD OF DRIVER FOR CHKDER EXAMPLE. C END SUBROUTINE FCN(M,N,X,FVEC,FJAC,LDFJAC,IFLAG) INTEGER M,N,LDFJAC,IFLAG REAL X(N),FVEC(M),FJAC(LDFJAC,N) C C SUBROUTINE FCN FOR CHKDER EXAMPLE. C INTEGER I 1 0 Page 0 REAL TMP1,TMP2,TMP3,TMP4 REAL Y(15) DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8), * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15) * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1, * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/ C IF (IFLAG .EQ. 2) GO TO 20 DO 10 I = 1, 15 TMP1 = I TMP2 = 16 - I TMP3 = TMP1 IF (I .GT. 8) TMP3 = TMP2 FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3)) 10 CONTINUE GO TO 40 20 CONTINUE DO 30 I = 1, 15 TMP1 = I TMP2 = 16 - I C C ERROR INTRODUCED INTO NEXT STATEMENT FOR ILLUSTRATION. C CORRECTED STATEMENT SHOULD READ TMP3 = TMP1 . C TMP3 = TMP2 IF (I .GT. 8) TMP3 = TMP2 TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2 FJAC(I,1) = -1.E0 FJAC(I,2) = TMP1*TMP2/TMP4 FJAC(I,3) = TMP1*TMP3/TMP4 30 CONTINUE 40 CONTINUE RETURN C C LAST CARD OF SUBROUTINE FCN. C END 0 Results obtained with different compilers or machines may be different. In particular, the differences FVECP - FVEC are machine dependent. 0 FVEC 0 -0.1181606E+01 -0.1429655E+01 -0.1606344E+01 -0.1745269E+01 -0.1840654E+01 -0.1921586E+01 -0.1984141E+01 -0.2022537E+01 -0.2468977E+01 -0.2827562E+01 -0.3473582E+01 -0.4437612E+01 -0.6047662E+01 -0.9267761E+01 -0.1891806E+02 0 FVECP - FVEC 0 -0.7724666E-08 -0.3432405E-08 -0.2034843E-09 0.2313685E-08 0.4331078E-08 0.5984096E-08 1 0 Page 0 0.7363281E-08 0.8531470E-08 0.1488591E-07 0.2335850E-07 0.3522012E-07 0.5301255E-07 0.8266660E-07 0.1419747E-06 0.3198990E-06 0 ERR 0 0.1141397E+00 0.9943516E-01 0.9674474E-01 0.9980447E-01 0.1073116E+00 0.1220445E+00 0.1526814E+00 0.1000000E+01 0.1000000E+01 0.1000000E+01 0.1000000E+01 0.1000000E+01 0.1000000E+01 0.1000000E+01 0.1000000E+01