+++ /dev/null
-*DECK SNLS1
- SUBROUTINE SNLS1 (FCN, IOPT, M, N, X, FVEC, FJAC, LDFJAC, FTOL,
- + XTOL, GTOL, MAXFEV, EPSFCN, DIAG, MODE, FACTOR, NPRINT, INFO,
- + NFEV, NJEV, IPVT, QTF, WA1, WA2, WA3, WA4)
-C***BEGIN PROLOGUE SNLS1
-C***PURPOSE Minimize the sum of the squares of M nonlinear functions
-C in N variables by a modification of the Levenberg-Marquardt
-C algorithm.
-C***LIBRARY SLATEC
-C***CATEGORY K1B1A1, K1B1A2
-C***TYPE SINGLE PRECISION (SNLS1-S, DNLS1-D)
-C***KEYWORDS LEVENBERG-MARQUARDT, NONLINEAR DATA FITTING,
-C NONLINEAR LEAST SQUARES
-C***AUTHOR Hiebert, K. L., (SNLA)
-C***DESCRIPTION
-C
-C 1. Purpose.
-C
-C The purpose of SNLS1 is to minimize the sum of the squares of M
-C nonlinear functions in N variables by a modification of the
-C Levenberg-Marquardt algorithm. The user must provide a subrou-
-C tine which calculates the functions. The user has the option
-C of how the Jacobian will be supplied. The user can supply the
-C full Jacobian, or the rows of the Jacobian (to avoid storing
-C the full Jacobian), or let the code approximate the Jacobian by
-C forward-differencing. This code is the combination of the
-C MINPACK codes (Argonne) LMDER, LMDIF, and LMSTR.
-C
-C
-C 2. Subroutine and Type Statements.
-C
-C SUBROUTINE SNLS1(FCN,IOPT,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,
-C * GTOL,MAXFEV,EPSFCN,DIAG,MODE,FACTOR,NPRINT,INFO
-C * ,NFEV,NJEV,IPVT,QTF,WA1,WA2,WA3,WA4)
-C INTEGER IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV
-C INTEGER IPVT(N)
-C REAL FTOL,XTOL,GTOL,EPSFCN,FACTOR
-C REAL X(N),FVEC(M),FJAC(LDFJAC,N),DIAG(N),QTF(N),
-C * WA1(N),WA2(N),WA3(N),WA4(M)
-C
-C
-C 3. Parameters.
-C
-C Parameters designated as input parameters must be specified on
-C entry to SNLS1 and are not changed on exit, while parameters
-C designated as output parameters need not be specified on entry
-C and are set to appropriate values on exit from SNLS1.
-C
-C FCN is the name of the user-supplied subroutine which calculates
-C the functions. If the user wants to supply the Jacobian
-C (IOPT=2 or 3), then FCN must be written to calculate the
-C Jacobian, as well as the functions. See the explanation
-C of the IOPT argument below.
-C If the user wants the iterates printed (NPRINT positive), then
-C FCN must do the printing. See the explanation of NPRINT
-C below. FCN must be declared in an EXTERNAL statement in the
-C calling program and should be written as follows.
-C
-C
-C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
-C INTEGER IFLAG,LDFJAC,M,N
-C REAL X(N),FVEC(M)
-C ----------
-C FJAC and LDFJAC may be ignored , if IOPT=1.
-C REAL FJAC(LDFJAC,N) , if IOPT=2.
-C REAL FJAC(N) , if IOPT=3.
-C ----------
-C If IFLAG=0, the values in X and FVEC are available
-C for printing. See the explanation of NPRINT below.
-C IFLAG will never be zero unless NPRINT is positive.
-C The values of X and FVEC must not be changed.
-C RETURN
-C ----------
-C If IFLAG=1, calculate the functions at X and return
-C this vector in FVEC.
-C RETURN
-C ----------
-C If IFLAG=2, calculate the full Jacobian at X and return
-C this matrix in FJAC. Note that IFLAG will never be 2 unless
-C IOPT=2. FVEC contains the function values at X and must
-C not be altered. FJAC(I,J) must be set to the derivative
-C of FVEC(I) with respect to X(J).
-C RETURN
-C ----------
-C If IFLAG=3, calculate the LDFJAC-th row of the Jacobian
-C and return this vector in FJAC. Note that IFLAG will
-C never be 3 unless IOPT=3. FVEC contains the function
-C values at X and must not be altered. FJAC(J) must be
-C set to the derivative of FVEC(LDFJAC) with respect to X(J).
-C RETURN
-C ----------
-C END
-C
-C
-C The value of IFLAG should not be changed by FCN unless the
-C user wants to terminate execution of SNLS1. In this case, set
-C IFLAG to a negative integer.
-C
-C
-C IOPT is an input variable which specifies how the Jacobian will
-C be calculated. If IOPT=2 or 3, then the user must supply the
-C Jacobian, as well as the function values, through the
-C subroutine FCN. If IOPT=2, the user supplies the full
-C Jacobian with one call to FCN. If IOPT=3, the user supplies
-C one row of the Jacobian with each call. (In this manner,
-C storage can be saved because the full Jacobian is not stored.)
-C If IOPT=1, the code will approximate the Jacobian by forward
-C differencing.
-C
-C M is a positive integer input variable set to the number of
-C functions.
-C
-C N is a positive integer input variable set to the number of
-C variables. N must not exceed M.
-C
-C X is an array of length N. On input, X must contain an initial
-C estimate of the solution vector. On output, X contains the
-C final estimate of the solution vector.
-C
-C FVEC is an output array of length M which contains the functions
-C evaluated at the output X.
-C
-C FJAC is an output array. For IOPT=1 and 2, FJAC is an M by N
-C array. For IOPT=3, FJAC is an N by N array. The upper N by N
-C submatrix of FJAC contains an upper triangular matrix R with
-C diagonal elements of nonincreasing magnitude such that
-C
-C T T T
-C P *(JAC *JAC)*P = R *R,
-C
-C where P is a permutation matrix and JAC is the final calcu-
-C lated Jacobian. Column J of P is column IPVT(J) (see below)
-C of the identity matrix. The lower part of FJAC contains
-C information generated during the computation of R.
-C
-C LDFJAC is a positive integer input variable which specifies
-C the leading dimension of the array FJAC. For IOPT=1 and 2,
-C LDFJAC must not be less than M. For IOPT=3, LDFJAC must not
-C be less than N.
-C
-C FTOL is a non-negative input variable. Termination occurs when
-C both the actual and predicted relative reductions in the sum
-C of squares are at most FTOL. Therefore, FTOL measures the
-C relative error desired in the sum of squares. Section 4 con-
-C tains more details about FTOL.
-C
-C XTOL is a non-negative input variable. Termination occurs when
-C the relative error between two consecutive iterates is at most
-C XTOL. Therefore, XTOL measures the relative error desired in
-C the approximate solution. Section 4 contains more details
-C about XTOL.
-C
-C GTOL is a non-negative input variable. Termination occurs when
-C the cosine of the angle between FVEC and any column of the
-C Jacobian is at most GTOL in absolute value. Therefore, GTOL
-C measures the orthogonality desired between the function vector
-C and the columns of the Jacobian. Section 4 contains more
-C details about GTOL.
-C
-C MAXFEV is a positive integer input variable. Termination occurs
-C when the number of calls to FCN to evaluate the functions
-C has reached MAXFEV.
-C
-C EPSFCN is an input variable used in determining a suitable step
-C for the forward-difference approximation. This approximation
-C assumes that the relative errors in the functions are of the
-C order of EPSFCN. If EPSFCN is less than the machine preci-
-C sion, it is assumed that the relative errors in the functions
-C are of the order of the machine precision. If IOPT=2 or 3,
-C then EPSFCN can be ignored (treat it as a dummy argument).
-C
-C DIAG is an array of length N. If MODE = 1 (see below), DIAG is
-C internally set. If MODE = 2, DIAG must contain positive
-C entries that serve as implicit (multiplicative) scale factors
-C for the variables.
-C
-C MODE is an integer input variable. If MODE = 1, the variables
-C will be scaled internally. If MODE = 2, the scaling is speci-
-C fied by the input DIAG. Other values of MODE are equivalent
-C to MODE = 1.
-C
-C FACTOR is a positive input variable used in determining the ini-
-C tial step bound. This bound is set to the product of FACTOR
-C and the Euclidean norm of DIAG*X if nonzero, or else to FACTOR
-C itself. In most cases FACTOR should lie in the interval
-C (.1,100.). 100. is a generally recommended value.
-C
-C NPRINT is an integer input variable that enables controlled
-C printing of iterates if it is positive. In this case, FCN is
-C called with IFLAG = 0 at the beginning of the first iteration
-C and every NPRINT iterations thereafter and immediately prior
-C to return, with X and FVEC available for printing. Appropriate
-C print statements must be added to FCN (see example) and
-C FVEC should not be altered. If NPRINT is not positive, no
-C special calls to FCN with IFLAG = 0 are made.
-C
-C INFO is an integer output variable. If the user has terminated
-C execution, INFO is set to the (negative) value of IFLAG. See
-C description of FCN and JAC. Otherwise, INFO is set as follows.
-C
-C INFO = 0 improper input parameters.
-C
-C INFO = 1 both actual and predicted relative reductions in the
-C sum of squares are at most FTOL.
-C
-C INFO = 2 relative error between two consecutive iterates is
-C at most XTOL.
-C
-C INFO = 3 conditions for INFO = 1 and INFO = 2 both hold.
-C
-C INFO = 4 the cosine of the angle between FVEC and any column
-C of the Jacobian is at most GTOL in absolute value.
-C
-C INFO = 5 number of calls to FCN for function evaluation
-C has reached MAXFEV.
-C
-C INFO = 6 FTOL is too small. No further reduction in the sum
-C of squares is possible.
-C
-C INFO = 7 XTOL is too small. No further improvement in the
-C approximate solution X is possible.
-C
-C INFO = 8 GTOL is too small. FVEC is orthogonal to the
-C columns of the Jacobian to machine precision.
-C
-C Sections 4 and 5 contain more details about INFO.
-C
-C NFEV is an integer output variable set to the number of calls to
-C FCN for function evaluation.
-C
-C NJEV is an integer output variable set to the number of
-C evaluations of the full Jacobian. If IOPT=2, only one call to
-C FCN is required for each evaluation of the full Jacobian.
-C If IOPT=3, the M calls to FCN are required.
-C If IOPT=1, then NJEV is set to zero.
-C
-C IPVT is an integer output array of length N. IPVT defines a
-C permutation matrix P such that JAC*P = Q*R, where JAC is the
-C final calculated Jacobian, Q is orthogonal (not stored), and R
-C is upper triangular with diagonal elements of nonincreasing
-C magnitude. Column J of P is column IPVT(J) of the identity
-C matrix.
-C
-C QTF is an output array of length N which contains the first N
-C elements of the vector (Q transpose)*FVEC.
-C
-C WA1, WA2, and WA3 are work arrays of length N.
-C
-C WA4 is a work array of length M.
-C
-C
-C 4. Successful Completion.
-C
-C The accuracy of SNLS1 is controlled by the convergence parame-
-C ters FTOL, XTOL, and GTOL. These parameters are used in tests
-C which make three types of comparisons between the approximation
-C X and a solution XSOL. SNLS1 terminates when any of the tests
-C is satisfied. If any of the convergence parameters is less than
-C the machine precision (as defined by the function R1MACH(4)),
-C then SNLS1 only attempts to satisfy the test defined by the
-C machine precision. Further progress is not usually possible.
-C
-C The tests assume that the functions are reasonably well behaved,
-C and, if the Jacobian is supplied by the user, that the functions
-C and the Jacobian are coded consistently. If these conditions
-C are not satisfied, then SNLS1 may incorrectly indicate conver-
-C gence. If the Jacobian is coded correctly or IOPT=1,
-C then the validity of the answer can be checked, for example, by
-C rerunning SNLS1 with tighter tolerances.
-C
-C First Convergence Test. If ENORM(Z) denotes the Euclidean norm
-C of a vector Z, then this test attempts to guarantee that
-C
-C ENORM(FVEC) .LE. (1+FTOL)*ENORM(FVECS),
-C
-C where FVECS denotes the functions evaluated at XSOL. If this
-C condition is satisfied with FTOL = 10**(-K), then the final
-C residual norm ENORM(FVEC) has K significant decimal digits and
-C INFO is set to 1 (or to 3 if the second test is also satis-
-C fied). Unless high precision solutions are required, the
-C recommended value for FTOL is the square root of the machine
-C precision.
-C
-C Second Convergence Test. If D is the diagonal matrix whose
-C entries are defined by the array DIAG, then this test attempts
-C to guarantee that
-C
-C ENORM(D*(X-XSOL)) .LE. XTOL*ENORM(D*XSOL).
-C
-C If this condition is satisfied with XTOL = 10**(-K), then the
-C larger components of D*X have K significant decimal digits and
-C INFO is set to 2 (or to 3 if the first test is also satis-
-C fied). There is a danger that the smaller components of D*X
-C may have large relative errors, but if MODE = 1, then the
-C accuracy of the components of X is usually related to their
-C sensitivity. Unless high precision solutions are required,
-C the recommended value for XTOL is the square root of the
-C machine precision.
-C
-C Third Convergence Test. This test is satisfied when the cosine
-C of the angle between FVEC and any column of the Jacobian at X
-C is at most GTOL in absolute value. There is no clear rela-
-C tionship between this test and the accuracy of SNLS1, and
-C furthermore, the test is equally well satisfied at other crit-
-C ical points, namely maximizers and saddle points. Therefore,
-C termination caused by this test (INFO = 4) should be examined
-C carefully. The recommended value for GTOL is zero.
-C
-C
-C 5. Unsuccessful Completion.
-C
-C Unsuccessful termination of SNLS1 can be due to improper input
-C parameters, arithmetic interrupts, or an excessive number of
-C function evaluations.
-C
-C Improper Input Parameters. INFO is set to 0 if IOPT .LT. 1
-C or IOPT .GT. 3, or N .LE. 0, or M .LT. N, or for IOPT=1 or 2
-C LDFJAC .LT. M, or for IOPT=3 LDFJAC .LT. N, or FTOL .LT. 0.E0,
-C or XTOL .LT. 0.E0, or GTOL .LT. 0.E0, or MAXFEV .LE. 0, or
-C FACTOR .LE. 0.E0.
-C
-C Arithmetic Interrupts. If these interrupts occur in the FCN
-C subroutine during an early stage of the computation, they may
-C be caused by an unacceptable choice of X by SNLS1. In this
-C case, it may be possible to remedy the situation by rerunning
-C SNLS1 with a smaller value of FACTOR.
-C
-C Excessive Number of Function Evaluations. A reasonable value
-C for MAXFEV is 100*(N+1) for IOPT=2 or 3 and 200*(N+1) for
-C IOPT=1. If the number of calls to FCN reaches MAXFEV, then
-C this indicates that the routine is converging very slowly
-C as measured by the progress of FVEC, and INFO is set to 5.
-C In this case, it may be helpful to restart SNLS1 with MODE
-C set to 1.
-C
-C
-C 6. Characteristics of the Algorithm.
-C
-C SNLS1 is a modification of the Levenberg-Marquardt algorithm.
-C Two of its main characteristics involve the proper use of
-C implicitly scaled variables (if MODE = 1) and an optimal choice
-C for the correction. The use of implicitly scaled variables
-C achieves scale invariance of SNLS1 and limits the size of the
-C correction in any direction where the functions are changing
-C rapidly. The optimal choice of the correction guarantees (under
-C reasonable conditions) global convergence from starting points
-C far from the solution and a fast rate of convergence for
-C problems with small residuals.
-C
-C Timing. The time required by SNLS1 to solve a given problem
-C depends on M and N, the behavior of the functions, the accu-
-C racy requested, and the starting point. The number of arith-
-C metic operations needed by SNLS1 is about N**3 to process each
-C evaluation of the functions (call to FCN) and to process each
-C evaluation of the Jacobian it takes M*N**2 for IOPT=2 (one
-C call to FCN), M*N**2 for IOPT=1 (N calls to FCN) and
-C 1.5*M*N**2 for IOPT=3 (M calls to FCN). Unless FCN
-C can be evaluated quickly, the timing of SNLS1 will be
-C strongly influenced by the time spent in FCN.
-C
-C Storage. SNLS1 requires (M*N + 2*M + 6*N) for IOPT=1 or 2 and
-C (N**2 + 2*M + 6*N) for IOPT=3 single precision storage
-C locations and N integer storage locations, in addition to
-C the storage required by the program. There are no internally
-C declared storage arrays.
-C
-C *Long Description:
-C
-C 7. Example.
-C
-C The problem is to determine the values of X(1), X(2), and X(3)
-C which provide the best fit (in the least squares sense) of
-C
-C X(1) + U(I)/(V(I)*X(2) + W(I)*X(3)), I = 1, 15
-C
-C to the data
-C
-C Y = (0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,
-C 0.37,0.58,0.73,0.96,1.34,2.10,4.39),
-C
-C where U(I) = I, V(I) = 16 - I, and W(I) = MIN(U(I),V(I)). The
-C I-th component of FVEC is thus defined by
-C
-C Y(I) - (X(1) + U(I)/(V(I)*X(2) + W(I)*X(3))).
-C
-C **********
-C
-C PROGRAM TEST
-C C
-C C Driver for SNLS1 example.
-C C
-C INTEGER J,IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,
-C * NWRITE
-C INTEGER IPVT(3)
-C REAL FTOL,XTOL,GTOL,FACTOR,FNORM,EPSFCN
-C REAL X(3),FVEC(15),FJAC(15,3),DIAG(3),QTF(3),
-C * WA1(3),WA2(3),WA3(3),WA4(15)
-C REAL ENORM,R1MACH
-C EXTERNAL FCN
-C DATA NWRITE /6/
-C C
-C IOPT = 1
-C M = 15
-C N = 3
-C C
-C C The following starting values provide a rough fit.
-C C
-C X(1) = 1.E0
-C X(2) = 1.E0
-C X(3) = 1.E0
-C C
-C LDFJAC = 15
-C C
-C C Set FTOL and XTOL to the square root of the machine precision
-C C and GTOL to zero. Unless high precision solutions are
-C C required, these are the recommended settings.
-C C
-C FTOL = SQRT(R1MACH(4))
-C XTOL = SQRT(R1MACH(4))
-C GTOL = 0.E0
-C C
-C MAXFEV = 400
-C EPSFCN = 0.0
-C MODE = 1
-C FACTOR = 1.E2
-C NPRINT = 0
-C C
-C CALL SNLS1(FCN,IOPT,M,N,X,FVEC,FJAC,LDFJAC,FTOL,XTOL,
-C * GTOL,MAXFEV,EPSFCN,DIAG,MODE,FACTOR,NPRINT,
-C * INFO,NFEV,NJEV,IPVT,QTF,WA1,WA2,WA3,WA4)
-C FNORM = ENORM(M,FVEC)
-C WRITE (NWRITE,1000) FNORM,NFEV,NJEV,INFO,(X(J),J=1,N)
-C STOP
-C 1000 FORMAT (5X,' FINAL L2 NORM OF THE RESIDUALS',E15.7 //
-C * 5X,' NUMBER OF FUNCTION EVALUATIONS',I10 //
-C * 5X,' NUMBER OF JACOBIAN EVALUATIONS',I10 //
-C * 5X,' EXIT PARAMETER',16X,I10 //
-C * 5X,' FINAL APPROXIMATE SOLUTION' // 5X,3E15.7)
-C END
-C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,DUM,IDUM)
-C C This is the form of the FCN routine if IOPT=1,
-C C that is, if the user does not calculate the Jacobian.
-C INTEGER M,N,IFLAG
-C REAL X(N),FVEC(M)
-C INTEGER I
-C REAL TMP1,TMP2,TMP3,TMP4
-C REAL Y(15)
-C DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
-C * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
-C * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1,
-C * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
-C C
-C IF (IFLAG .NE. 0) GO TO 5
-C C
-C C Insert print statements here when NPRINT is positive.
-C C
-C RETURN
-C 5 CONTINUE
-C DO 10 I = 1, M
-C TMP1 = I
-C TMP2 = 16 - I
-C TMP3 = TMP1
-C IF (I .GT. 8) TMP3 = TMP2
-C FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
-C 10 CONTINUE
-C RETURN
-C END
-C
-C
-C Results obtained with different compilers or machines
-C may be slightly different.
-C
-C FINAL L2 NORM OF THE RESIDUALS 0.9063596E-01
-C
-C NUMBER OF FUNCTION EVALUATIONS 25
-C
-C NUMBER OF JACOBIAN EVALUATIONS 0
-C
-C EXIT PARAMETER 1
-C
-C FINAL APPROXIMATE SOLUTION
-C
-C 0.8241058E-01 0.1133037E+01 0.2343695E+01
-C
-C
-C For IOPT=2, FCN would be modified as follows to also
-C calculate the full Jacobian when IFLAG=2.
-C
-C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
-C C
-C C This is the form of the FCN routine if IOPT=2,
-C C that is, if the user calculates the full Jacobian.
-C C
-C INTEGER LDFJAC,M,N,IFLAG
-C REAL X(N),FVEC(M)
-C REAL FJAC(LDFJAC,N)
-C INTEGER I
-C REAL TMP1,TMP2,TMP3,TMP4
-C REAL Y(15)
-C DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
-C * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
-C * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1,
-C * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
-C C
-C IF (IFLAG .NE. 0) GO TO 5
-C C
-C C Insert print statements here when NPRINT is positive.
-C C
-C RETURN
-C 5 CONTINUE
-C IF(IFLAG.NE.1) GO TO 20
-C DO 10 I = 1, M
-C TMP1 = I
-C TMP2 = 16 - I
-C TMP3 = TMP1
-C IF (I .GT. 8) TMP3 = TMP2
-C FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
-C 10 CONTINUE
-C RETURN
-C C
-C C Below, calculate the full Jacobian.
-C C
-C 20 CONTINUE
-C C
-C DO 30 I = 1, M
-C TMP1 = I
-C TMP2 = 16 - I
-C TMP3 = TMP1
-C IF (I .GT. 8) TMP3 = TMP2
-C TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2
-C FJAC(I,1) = -1.E0
-C FJAC(I,2) = TMP1*TMP2/TMP4
-C FJAC(I,3) = TMP1*TMP3/TMP4
-C 30 CONTINUE
-C RETURN
-C END
-C
-C
-C For IOPT = 3, FJAC would be dimensioned as FJAC(3,3),
-C LDFJAC would be set to 3, and FCN would be written as
-C follows to calculate a row of the Jacobian when IFLAG=3.
-C
-C SUBROUTINE FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
-C C This is the form of the FCN routine if IOPT=3,
-C C that is, if the user calculates the Jacobian row by row.
-C INTEGER M,N,IFLAG
-C REAL X(N),FVEC(M)
-C REAL FJAC(N)
-C INTEGER I
-C REAL TMP1,TMP2,TMP3,TMP4
-C REAL Y(15)
-C DATA Y(1),Y(2),Y(3),Y(4),Y(5),Y(6),Y(7),Y(8),
-C * Y(9),Y(10),Y(11),Y(12),Y(13),Y(14),Y(15)
-C * /1.4E-1,1.8E-1,2.2E-1,2.5E-1,2.9E-1,3.2E-1,3.5E-1,3.9E-1,
-C * 3.7E-1,5.8E-1,7.3E-1,9.6E-1,1.34E0,2.1E0,4.39E0/
-C C
-C IF (IFLAG .NE. 0) GO TO 5
-C C
-C C Insert print statements here when NPRINT is positive.
-C C
-C RETURN
-C 5 CONTINUE
-C IF( IFLAG.NE.1) GO TO 20
-C DO 10 I = 1, M
-C TMP1 = I
-C TMP2 = 16 - I
-C TMP3 = TMP1
-C IF (I .GT. 8) TMP3 = TMP2
-C FVEC(I) = Y(I) - (X(1) + TMP1/(X(2)*TMP2 + X(3)*TMP3))
-C 10 CONTINUE
-C RETURN
-C C
-C C Below, calculate the LDFJAC-th row of the Jacobian.
-C C
-C 20 CONTINUE
-C
-C I = LDFJAC
-C TMP1 = I
-C TMP2 = 16 - I
-C TMP3 = TMP1
-C IF (I .GT. 8) TMP3 = TMP2
-C TMP4 = (X(2)*TMP2 + X(3)*TMP3)**2
-C FJAC(1) = -1.E0
-C FJAC(2) = TMP1*TMP2/TMP4
-C FJAC(3) = TMP1*TMP3/TMP4
-C RETURN
-C END
-C
-C***REFERENCES Jorge J. More, The Levenberg-Marquardt algorithm:
-C implementation and theory. In Numerical Analysis
-C Proceedings (Dundee, June 28 - July 1, 1977, G. A.
-C Watson, Editor), Lecture Notes in Mathematics 630,
-C Springer-Verlag, 1978.
-C***ROUTINES CALLED CHKDER, ENORM, FDJAC3, LMPAR, QRFAC, R1MACH,
-C RWUPDT, XERMSG
-C***REVISION HISTORY (YYMMDD)
-C 800301 DATE WRITTEN
-C 890531 Changed all specific intrinsics to generic. (WRB)
-C 890531 REVISION DATE from Version 3.2
-C 891214 Prologue converted to Version 4.0 format. (BAB)
-C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
-C 900510 Convert XERRWV calls to XERMSG calls. (RWC)
-C 920501 Reformatted the REFERENCES section. (WRB)
-C***END PROLOGUE SNLS1
- INTEGER IOPT,M,N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV
- INTEGER IJUNK,NROW,IPVT(*)
- REAL FTOL,XTOL,GTOL,FACTOR,EPSFCN
- REAL X(*),FVEC(*),FJAC(LDFJAC,*),DIAG(*),QTF(*),WA1(*),WA2(*),
- 1 WA3(*),WA4(*)
- LOGICAL SING
- EXTERNAL FCN
- INTEGER I,IFLAG,ITER,J,L,MODECH
- REAL ACTRED,DELTA,DIRDER,EPSMCH,FNORM,FNORM1,GNORM,ONE,PAR,
- 1 PNORM,PRERED,P1,P5,P25,P75,P0001,RATIO,SUM,TEMP,TEMP1,
- 2 TEMP2,XNORM,ZERO
- REAL R1MACH,ENORM,ERR,CHKLIM
- CHARACTER*8 XERN1
- CHARACTER*16 XERN3
-C
- SAVE CHKLIM, ONE, P1, P5, P25, P75, P0001, ZERO
- DATA CHKLIM/.1E0/
- DATA ONE,P1,P5,P25,P75,P0001,ZERO
- 1 /1.0E0,1.0E-1,5.0E-1,2.5E-1,7.5E-1,1.0E-4,0.0E0/
-C
-C***FIRST EXECUTABLE STATEMENT SNLS1
- EPSMCH = R1MACH(4)
-C
- INFO = 0
- IFLAG = 0
- NFEV = 0
- NJEV = 0
-C
-C CHECK THE INPUT PARAMETERS FOR ERRORS.
-C
- IF (IOPT .LT. 1 .OR. IOPT .GT. 3 .OR. N .LE. 0 .OR.
- 1 M .LT. N .OR. LDFJAC .LT. N .OR. FTOL .LT. ZERO
- 2 .OR. XTOL .LT. ZERO .OR. GTOL .LT. ZERO
- 3 .OR. MAXFEV .LE. 0 .OR. FACTOR .LE. ZERO) GO TO 300
- IF (IOPT .LT. 3 .AND. LDFJAC .LT. M) GO TO 300
- IF (MODE .NE. 2) GO TO 20
- DO 10 J = 1, N
- IF (DIAG(J) .LE. ZERO) GO TO 300
- 10 CONTINUE
- 20 CONTINUE
-C
-C EVALUATE THE FUNCTION AT THE STARTING POINT
-C AND CALCULATE ITS NORM.
-C
- IFLAG = 1
- IJUNK = 1
- CALL FCN(IFLAG,M,N,X,FVEC,FJAC,IJUNK)
- NFEV = 1
- IF (IFLAG .LT. 0) GO TO 300
- FNORM = ENORM(M,FVEC)
-C
-C INITIALIZE LEVENBERG-MARQUARDT PARAMETER AND ITERATION COUNTER.
-C
- PAR = ZERO
- ITER = 1
-C
-C BEGINNING OF THE OUTER LOOP.
-C
- 30 CONTINUE
-C
-C IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.
-C
- IF (NPRINT .LE. 0) GO TO 40
- IFLAG = 0
- IF (MOD(ITER-1,NPRINT) .EQ. 0)
- 1 CALL FCN(IFLAG,M,N,X,FVEC,FJAC,IJUNK)
- IF (IFLAG .LT. 0) GO TO 300
- 40 CONTINUE
-C
-C CALCULATE THE JACOBIAN MATRIX.
-C
- IF (IOPT .EQ. 3) GO TO 475
-C
-C STORE THE FULL JACOBIAN USING M*N STORAGE
-C
- IF (IOPT .EQ. 1) GO TO 410
-C
-C THE USER SUPPLIES THE JACOBIAN
-C
- IFLAG = 2
- CALL FCN(IFLAG,M,N,X,FVEC,FJAC,LDFJAC)
- NJEV = NJEV + 1
-C
-C ON THE FIRST ITERATION, CHECK THE USER SUPPLIED JACOBIAN
-C
- IF (ITER .LE. 1) THEN
- IF (IFLAG .LT. 0) GO TO 300
-C
-C GET THE INCREMENTED X-VALUES INTO WA1(*).
-C
- MODECH = 1
- CALL CHKDER(M,N,X,FVEC,FJAC,LDFJAC,WA1,WA4,MODECH,ERR)
-C
-C EVALUATE FUNCTION AT INCREMENTED VALUE AND PUT IN WA4(*).
-C
- IFLAG = 1
- CALL FCN(IFLAG,M,N,WA1,WA4,FJAC,LDFJAC)
- NFEV = NFEV + 1
- IF(IFLAG .LT. 0) GO TO 300
- DO 350 I = 1, M
- MODECH = 2
- CALL CHKDER(1,N,X,FVEC(I),FJAC(I,1),LDFJAC,WA1,
- 1 WA4(I),MODECH,ERR)
- IF (ERR .LT. CHKLIM) THEN
- WRITE (XERN1, '(I8)') I
- WRITE (XERN3, '(1PE15.6)') ERR
- CALL XERMSG ('SLATEC', 'SNLS1', 'DERIVATIVE OF ' //
- * 'FUNCTION ' // XERN1 // ' MAY BE WRONG, ERR = ' //
- * XERN3 // ' TOO CLOSE TO 0.', 7, 0)
- ENDIF
- 350 CONTINUE
- ENDIF
-C
- GO TO 420
-C
-C THE CODE APPROXIMATES THE JACOBIAN
-C
-410 IFLAG = 1
- CALL FDJAC3(FCN,M,N,X,FVEC,FJAC,LDFJAC,IFLAG,EPSFCN,WA4)
- NFEV = NFEV + N
- 420 IF (IFLAG .LT. 0) GO TO 300
-C
-C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.
-C
- CALL QRFAC(M,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3)
-C
-C FORM (Q TRANSPOSE)*FVEC AND STORE THE FIRST N COMPONENTS IN
-C QTF.
-C
- DO 430 I = 1, M
- WA4(I) = FVEC(I)
- 430 CONTINUE
- DO 470 J = 1, N
- IF (FJAC(J,J) .EQ. ZERO) GO TO 460
- SUM = ZERO
- DO 440 I = J, M
- SUM = SUM + FJAC(I,J)*WA4(I)
- 440 CONTINUE
- TEMP = -SUM/FJAC(J,J)
- DO 450 I = J, M
- WA4(I) = WA4(I) + FJAC(I,J)*TEMP
- 450 CONTINUE
- 460 CONTINUE
- FJAC(J,J) = WA1(J)
- QTF(J) = WA4(J)
- 470 CONTINUE
- GO TO 560
-C
-C ACCUMULATE THE JACOBIAN BY ROWS IN ORDER TO SAVE STORAGE.
-C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN MATRIX
-C CALCULATED ONE ROW AT A TIME, WHILE SIMULTANEOUSLY
-C FORMING (Q TRANSPOSE)*FVEC AND STORING THE FIRST
-C N COMPONENTS IN QTF.
-C
- 475 DO 490 J = 1, N
- QTF(J) = ZERO
- DO 480 I = 1, N
- FJAC(I,J) = ZERO
- 480 CONTINUE
- 490 CONTINUE
- DO 500 I = 1, M
- NROW = I
- IFLAG = 3
- CALL FCN(IFLAG,M,N,X,FVEC,WA3,NROW)
- IF (IFLAG .LT. 0) GO TO 300
-C
-C ON THE FIRST ITERATION, CHECK THE USER SUPPLIED JACOBIAN.
-C
- IF(ITER .GT. 1) GO TO 498
-C
-C GET THE INCREMENTED X-VALUES INTO WA1(*).
-C
- MODECH = 1
- CALL CHKDER(M,N,X,FVEC,FJAC,LDFJAC,WA1,WA4,MODECH,ERR)
-C
-C EVALUATE AT INCREMENTED VALUES, IF NOT ALREADY EVALUATED.
-C
- IF(I .NE. 1) GO TO 495
-C
-C EVALUATE FUNCTION AT INCREMENTED VALUE AND PUT INTO WA4(*).
-C
- IFLAG = 1
- CALL FCN(IFLAG,M,N,WA1,WA4,FJAC,NROW)
- NFEV = NFEV + 1
- IF(IFLAG .LT. 0) GO TO 300
-495 CONTINUE
- MODECH = 2
- CALL CHKDER(1,N,X,FVEC(I),WA3,1,WA1,WA4(I),MODECH,ERR)
- IF (ERR .LT. CHKLIM) THEN
- WRITE (XERN1, '(I8)') I
- WRITE (XERN3, '(1PE15.6)') ERR
- CALL XERMSG ('SLATEC', 'SNLS1', 'DERIVATIVE OF FUNCTION '
- * // XERN1 // ' MAY BE WRONG, ERR = ' // XERN3 //
- * ' TOO CLOSE TO 0.', 7, 0)
- ENDIF
-498 CONTINUE
-C
- TEMP = FVEC(I)
- CALL RWUPDT(N,FJAC,LDFJAC,WA3,QTF,TEMP,WA1,WA2)
- 500 CONTINUE
- NJEV = NJEV + 1
-C
-C IF THE JACOBIAN IS RANK DEFICIENT, CALL QRFAC TO
-C REORDER ITS COLUMNS AND UPDATE THE COMPONENTS OF QTF.
-C
- SING = .FALSE.
- DO 510 J = 1, N
- IF (FJAC(J,J) .EQ. ZERO) SING = .TRUE.
- IPVT(J) = J
- WA2(J) = ENORM(J,FJAC(1,J))
- 510 CONTINUE
- IF (.NOT.SING) GO TO 560
- CALL QRFAC(N,N,FJAC,LDFJAC,.TRUE.,IPVT,N,WA1,WA2,WA3)
- DO 550 J = 1, N
- IF (FJAC(J,J) .EQ. ZERO) GO TO 540
- SUM = ZERO
- DO 520 I = J, N
- SUM = SUM + FJAC(I,J)*QTF(I)
- 520 CONTINUE
- TEMP = -SUM/FJAC(J,J)
- DO 530 I = J, N
- QTF(I) = QTF(I) + FJAC(I,J)*TEMP
- 530 CONTINUE
- 540 CONTINUE
- FJAC(J,J) = WA1(J)
- 550 CONTINUE
- 560 CONTINUE
-C
-C ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING
-C TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.
-C
- IF (ITER .NE. 1) GO TO 80
- IF (MODE .EQ. 2) GO TO 60
- DO 50 J = 1, N
- DIAG(J) = WA2(J)
- IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE
- 50 CONTINUE
- 60 CONTINUE
-C
-C ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X
-C AND INITIALIZE THE STEP BOUND DELTA.
-C
- DO 70 J = 1, N
- WA3(J) = DIAG(J)*X(J)
- 70 CONTINUE
- XNORM = ENORM(N,WA3)
- DELTA = FACTOR*XNORM
- IF (DELTA .EQ. ZERO) DELTA = FACTOR
- 80 CONTINUE
-C
-C COMPUTE THE NORM OF THE SCALED GRADIENT.
-C
- GNORM = ZERO
- IF (FNORM .EQ. ZERO) GO TO 170
- DO 160 J = 1, N
- L = IPVT(J)
- IF (WA2(L) .EQ. ZERO) GO TO 150
- SUM = ZERO
- DO 140 I = 1, J
- SUM = SUM + FJAC(I,J)*(QTF(I)/FNORM)
- 140 CONTINUE
- GNORM = MAX(GNORM,ABS(SUM/WA2(L)))
- 150 CONTINUE
- 160 CONTINUE
- 170 CONTINUE
-C
-C TEST FOR CONVERGENCE OF THE GRADIENT NORM.
-C
- IF (GNORM .LE. GTOL) INFO = 4
- IF (INFO .NE. 0) GO TO 300
-C
-C RESCALE IF NECESSARY.
-C
- IF (MODE .EQ. 2) GO TO 190
- DO 180 J = 1, N
- DIAG(J) = MAX(DIAG(J),WA2(J))
- 180 CONTINUE
- 190 CONTINUE
-C
-C BEGINNING OF THE INNER LOOP.
-C
- 200 CONTINUE
-C
-C DETERMINE THE LEVENBERG-MARQUARDT PARAMETER.
-C
- CALL LMPAR(N,FJAC,LDFJAC,IPVT,DIAG,QTF,DELTA,PAR,WA1,WA2,
- 1 WA3,WA4)
-C
-C STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.
-C
- DO 210 J = 1, N
- WA1(J) = -WA1(J)
- WA2(J) = X(J) + WA1(J)
- WA3(J) = DIAG(J)*WA1(J)
- 210 CONTINUE
- PNORM = ENORM(N,WA3)
-C
-C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.
-C
- IF (ITER .EQ. 1) DELTA = MIN(DELTA,PNORM)
-C
-C EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.
-C
- IFLAG = 1
- CALL FCN(IFLAG,M,N,WA2,WA4,FJAC,IJUNK)
- NFEV = NFEV + 1
- IF (IFLAG .LT. 0) GO TO 300
- FNORM1 = ENORM(M,WA4)
-C
-C COMPUTE THE SCALED ACTUAL REDUCTION.
-C
- ACTRED = -ONE
- IF (P1*FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2
-C
-C COMPUTE THE SCALED PREDICTED REDUCTION AND
-C THE SCALED DIRECTIONAL DERIVATIVE.
-C
- DO 230 J = 1, N
- WA3(J) = ZERO
- L = IPVT(J)
- TEMP = WA1(L)
- DO 220 I = 1, J
- WA3(I) = WA3(I) + FJAC(I,J)*TEMP
- 220 CONTINUE
- 230 CONTINUE
- TEMP1 = ENORM(N,WA3)/FNORM
- TEMP2 = (SQRT(PAR)*PNORM)/FNORM
- PRERED = TEMP1**2 + TEMP2**2/P5
- DIRDER = -(TEMP1**2 + TEMP2**2)
-C
-C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED
-C REDUCTION.
-C
- RATIO = ZERO
- IF (PRERED .NE. ZERO) RATIO = ACTRED/PRERED
-C
-C UPDATE THE STEP BOUND.
-C
- IF (RATIO .GT. P25) GO TO 240
- IF (ACTRED .GE. ZERO) TEMP = P5
- IF (ACTRED .LT. ZERO)
- 1 TEMP = P5*DIRDER/(DIRDER + P5*ACTRED)
- IF (P1*FNORM1 .GE. FNORM .OR. TEMP .LT. P1) TEMP = P1
- DELTA = TEMP*MIN(DELTA,PNORM/P1)
- PAR = PAR/TEMP
- GO TO 260
- 240 CONTINUE
- IF (PAR .NE. ZERO .AND. RATIO .LT. P75) GO TO 250
- DELTA = PNORM/P5
- PAR = P5*PAR
- 250 CONTINUE
- 260 CONTINUE
-C
-C TEST FOR SUCCESSFUL ITERATION.
-C
- IF (RATIO .LT. P0001) GO TO 290
-C
-C SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.
-C
- DO 270 J = 1, N
- X(J) = WA2(J)
- WA2(J) = DIAG(J)*X(J)
- 270 CONTINUE
- DO 280 I = 1, M
- FVEC(I) = WA4(I)
- 280 CONTINUE
- XNORM = ENORM(N,WA2)
- FNORM = FNORM1
- ITER = ITER + 1
- 290 CONTINUE
-C
-C TESTS FOR CONVERGENCE.
-C
- IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL
- 1 .AND. P5*RATIO .LE. ONE) INFO = 1
- IF (DELTA .LE. XTOL*XNORM) INFO = 2
- IF (ABS(ACTRED) .LE. FTOL .AND. PRERED .LE. FTOL
- 1 .AND. P5*RATIO .LE. ONE .AND. INFO .EQ. 2) INFO = 3
- IF (INFO .NE. 0) GO TO 300
-C
-C TESTS FOR TERMINATION AND STRINGENT TOLERANCES.
-C
- IF (NFEV .GE. MAXFEV) INFO = 5
- IF (ABS(ACTRED) .LE. EPSMCH .AND. PRERED .LE. EPSMCH
- 1 .AND. P5*RATIO .LE. ONE) INFO = 6
- IF (DELTA .LE. EPSMCH*XNORM) INFO = 7
- IF (GNORM .LE. EPSMCH) INFO = 8
- IF (INFO .NE. 0) GO TO 300
-C
-C END OF THE INNER LOOP. REPEAT IF ITERATION UNSUCCESSFUL.
-C
- IF (RATIO .LT. P0001) GO TO 200
-C
-C END OF THE OUTER LOOP.
-C
- GO TO 30
- 300 CONTINUE
-C
-C TERMINATION, EITHER NORMAL OR USER IMPOSED.
-C
- IF (IFLAG .LT. 0) INFO = IFLAG
- IFLAG = 0
- IF (NPRINT .GT. 0) CALL FCN(IFLAG,M,N,X,FVEC,FJAC,IJUNK)
- IF (INFO .LT. 0) CALL XERMSG ('SLATEC', 'SNLS1',
- + 'EXECUTION TERMINATED BECAUSE USER SET IFLAG NEGATIVE.', 1, 1)
- IF (INFO .EQ. 0) CALL XERMSG ('SLATEC', 'SNLS1',
- + 'INVALID INPUT PARAMETER.', 2, 1)
- IF (INFO .EQ. 4) CALL XERMSG ('SLATEC', 'SNLS1',
- + 'THIRD CONVERGENCE CONDITION, CHECK RESULTS BEFORE ACCEPTING.',
- + 1, 1)
- IF (INFO .EQ. 5) CALL XERMSG ('SLATEC', 'SNLS1',
- + 'TOO MANY FUNCTION EVALUATIONS.', 9, 1)
- IF (INFO .GE. 6) CALL XERMSG ('SLATEC', 'SNLS1',
- + 'TOLERANCES TOO SMALL, NO FURTHER IMPROVEMENT POSSIBLE.', 3, 1)
- RETURN
-C
-C LAST CARD OF SUBROUTINE SNLS1.
-C
- END
git://source.jalview.org / jabaws.git / blobdiff