SUBROUTINE SNSQE(FCN,JAC,IOPT,N,X,FVEC,TOL,NPRINT,INFO,WA,LWA)
C***BEGIN PROLOGUE SNSQE
C***DATE WRITTEN 800301 (YYMMDD)
C***REVISION DATE 880222 (YYMMDD)
C***CATEGORY NO. F2A
C***KEYWORDS EASY-TO-USE,NONLINEAR SQUARE SYSTEM,POWELL HYBRID METHOD,
C ZERO
C***AUTHOR HIEBERT, K. L., (SNLA)
C***PURPOSE SNSQE is the easy-to-use version of SNSQ which finds a zero
C of a system of N nonlinear functions in N variables by a
C modification of Powell hybrid method. This code is the
C combination of the MINPACK codes(Argonne) HYBRD1 and HYBRJ1
C***DESCRIPTION
C
C 1. Purpose.
C
C
C The purpose of SNSQE is to find a zero of a system of N non-
C linear functions in N variables by a modification of the Powell
C hybrid method. This is done by using the more general nonlinear
C equation solver SNSQ. The user must provide a subroutine which
C calculates the functions. The user has the option of either to
C provide a subroutine which calculates the Jacobian or to let the
C code calculate it by a forward-difference approximation. This
C code is the combination of the MINPACK codes (Argonne) HYBRD1
C and HYBRJ1.
C
C
C 2. Subroutine and Type Statements.
C
C SUBROUTINE SNSQE(FCN,JAC,IOPT,N,X,FVEC,TOL,NPRINT,INFO,
C * WA,LWA)
C INTEGER IOPT,N,NPRINT,INFO,LWA
C REAL TOL
C REAL X(N),FVEC(N),WA(LWA)
C EXTERNAL FCN,JAC
C
C
C 3. Parameters.
C
C Parameters designated as input parameters must be specified on
C entry to SNSQE 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 SNSQE.
C
C FCN is the name of the user-supplied subroutine which calculates
C the functions. FCN must be declared in an EXTERNAL statement
C in the user calling program, and should be written as follows.
C
C SUBROUTINE FCN(N,X,FVEC,IFLAG)
C INTEGER N,IFLAG
C REAL X(N),FVEC(N)
C ----------
C Calculate the functions at X and
C return this vector in FVEC.
C ----------
C RETURN
C END
C
C The value of IFLAG should not be changed by FCN unless the
C user wants to terminate execution of SNSQE. In this case, set
C IFLAG to a negative integer.
C
C JAC is the name of the user-supplied subroutine which calculates
C the Jacobian. If IOPT=1, then JAC must be declared in an
C EXTERNAL statement in the user calling program, and should be
C written as follows.
C
C SUBROUTINE JAC(N,X,FVEC,FJAC,LDFJAC,IFLAG)
C INTEGER N,LDFJAC,IFLAG
C REAL X(N),FVEC(N),FJAC(LDFJAC,N)
C ----------
C Calculate the Jacobian at X and return this
C matrix in FJAC. FVEC contains the function
C values at X and should not be altered.
C ----------
C RETURN
C END
C
C The value of IFLAG should not be changed by JAC unless the
C user wants to terminate execution of SNSQE. In this case, set
C IFLAG to a negative integer.
C
C If IOPT=2, JAC can be ignored (treat it as a dummy argument).
C
C IOPT is an input variable which specifies how the Jacobian will
C be calculated. If IOPT=1, then the user must supply the
C Jacobian through the subroutine JAC. If IOPT=2, then the
C code will approximate the Jacobian by forward-differencing.
C
C N is a positive integer input variable set to the number of
C functions and variables.
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 N which contains the functions
C evaluated at the output X.
C
C TOL is a non-negative input variable. Termination occurs when
C the algorithm estimates that the relative error between X and
C the solution is at most TOL. Section 4 contains more details
C about TOL.
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 iteration 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). If NPRINT
C is not positive, no special calls of FCN with IFLAG = 0 are
C 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 algorithm estimates that the relative error between
C X and the solution is at most TOL.
C
C INFO = 2 number of calls to FCN has reached or exceeded
C 100*(N+1) for IOPT=1 or 200*(N+1) for IOPT=2.
C
C INFO = 3 TOL is too small. No further improvement in the
C approximate solution X is possible.
C
C INFO = 4 iteration is not making good progress.
C
C Sections 4 and 5 contain more details about INFO.
C
C WA is a work array of length LWA.
C
C LWA is a positive integer input variable not less than
C (3*N**2+13*N))/2.
C
C
C 4. Successful Completion.
C
C The accuracy of SNSQE is controlled by the convergence parame-
C ter TOL. This parameter is used in a test which makes a compar-
C ison between the approximation X and a solution XSOL. SNSQE
C terminates when the test is satisfied. If TOL is less than the
C machine precision (as defined by the function R1MACH(4)), then
C SNSQE attemps only to satisfy the test defined by the machine
C precision. Further progress is not usually possible. Unless
C high precision solutions are required, the recommended value
C for TOL is the square root of the machine precision.
C
C The test assumes that the functions are reasonably well behaved,
C and, if the Jacobian is supplied by the user, that the functions
C and the Jacobian coded consistently. If these conditions
C are not satisfied, SNSQE may incorrectly indicate convergence.
C The coding of the Jacobian can be checked by the subroutine
C CHKDER. If the Jacobian is coded correctly or IOPT=2, then
C the validity of the answer can be checked, for example, by
C rerunning SNSQE with a tighter tolerance.
C
C Convergence Test. If SNRM2(Z) denotes the Euclidean norm of a
C vector Z, then this test attempts to guarantee that
C
C SNRM2(X-XSOL) .LE. TOL*SNRM2(XSOL).
C
C If this condition is satisfied with TOL = 10**(-K), then the
C larger components of X have K significant decimal digits and
C INFO is set to 1. There is a danger that the smaller compo-
C nents of X may have large relative errors, but the fast rate
C of convergence of SNSQE usually avoids this possibility.
C
C
C 5. Unsuccessful Completion.
C
C Unsuccessful termination of SNSQE can be due to improper input
C parameters, arithmetic interrupts, an excessive number of func-
C tion evaluations, errors in the functions, or lack of good prog-
C ress.
C
C Improper Input Parameters. INFO is set to 0 if IOPT .LT. 1, or
C IOPT .GT. 2, or N .LE. 0, or TOL .LT. 0.E0, or
C LWA .LT. (3*N**2+13*N)/2.
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 SNSQE. In this
C case, it may be possible to remedy the situation by not evalu-
C ating the functions here, but instead setting the components
C of FVEC to numbers that exceed those in the initial FVEC.
C
C Excessive Number of Function Evaluations. If the number of
C calls to FCN reaches 100*(N+1) for IOPT=1 or 200*(N+1) for
C IOPT=2, then this indicates that the routine is converging
C very slowly as measured by the progress of FVEC, and INFO is
C set to 2. This situation should be unusual because, as
C indicated below, lack of good progress is usually diagnosed
C earlier by SNSQE, causing termination with INFO = 4.
C
C Errors in the Functions. When IOPT=2, the choice of step length
C in the forward-difference approximation to the Jacobian
C assumes that the relative errors in the functions are of the
C order of the machine precision. If this is not the case,
C SNSQE may fail (usually with INFO = 4). The user should
C then either use SNSQ and set the step length or use IOPT=1
C and supply the Jacobian.
C
C Lack of Good Progress. SNSQE searches for a zero of the system
C by minimizing the sum of the squares of the functions. In so
C doing, it can become trapped in a region where the minimum
C does not correspond to a zero of the system and, in this situ-
C ation, the iteration eventually fails to make good progress.
C In particular, this will happen if the system does not have a
C zero. If the system has a zero, rerunning SNSQE from a dif-
C ferent starting point may be helpful.
C
C
C 6. Characteristics of the Algorithm.
C
C SNSQE is a modification of the Powell hybrid method. Two of
C its main characteristics involve the choice of the correction as
C a convex combination of the Newton and scaled gradient direc-
C tions, and the updating of the Jacobian by the rank-1 method of
C Broyden. The choice of the correction guarantees (under reason-
C able conditions) global convergence for starting points far from
C the solution and a fast rate of convergence. The Jacobian is
C calculated at the starting point by either the user-supplied
C subroutine or a forward-difference approximation, but it is not
C recalculated until the rank-1 method fails to produce satis-
C factory progress.
C
C Timing. The time required by SNSQE to solve a given problem
C depends on N, the behavior of the functions, the accuracy
C requested, and the starting point. The number of arithmetic
C operations needed by SNSQE is about 11.5*(N**2) to process
C each evaluation of the functions (call to FCN) and 1.3*(N**3)
C to process each evaluation of the Jacobian (call to JAC,
C if IOPT = 1). Unless FCN and JAC can be evaluated quickly,
C the timing of SNSQE will be strongly influenced by the time
C spent in FCN and JAC.
C
C Storage. SNSQE requires (3*N**2 + 17*N)/2 single precision
C storage locations, in addition to the storage required by the
C program. There are no internally declared storage arrays.
C
C
C 7. Example.
C
C The problem is to determine the values of X(1), X(2), ..., X(9),
C which solve the system of tridiagonal equations
C
C (3-2*X(1))*X(1) -2*X(2) = -1
C -X(I-1) + (3-2*X(I))*X(I) -2*X(I+1) = -1, I=2-8
C -X(8) + (3-2*X(9))*X(9) = -1
C
C **********
C
C PROGRAM TEST(INPUT,OUTPUT,TAPE6=OUTPUT)
C C
C C Driver for SNSQE example.
C C
C INTEGER J,N,IOPT,NPRINT,INFO,LWA,NWRITE
C REAL TOL,FNORM
C REAL X(9),FVEC(9),WA(180)
C REAL SNRM2,R1MACH
C EXTERNAL FCN
C DATA NWRITE /6/
C C
C IOPT = 2
C N = 9
C C
C C The following starting values provide a rough solution.
C C
C DO 10 J = 1, 9
C X(J) = -1.E0
C 10 CONTINUE
C
C LWA = 180
C NPRINT = 0
C C
C C Set TOL to the square root of the machine precision.
C C Unless high precision solutions are required,
C C this is the recommended setting.
C C
C TOL = SQRT(R1MACH(4))
C C
C CALL SNSQE(FCN,JAC,IOPT,N,X,FVEC,TOL,NPRINT,INFO,WA,LWA)
C FNORM = SNRM2(N,FVEC)
C WRITE (NWRITE,1000) FNORM,INFO,(X(J),J=1,N)
C STOP
C 1000 FORMAT (5X,' FINAL L2 NORM OF THE RESIDUALS',E15.7 //
C * 5X,' EXIT PARAMETER',16X,I10 //
C * 5X,' FINAL APPROXIMATE SOLUTION' // (5X,3E15.7))
C END
C SUBROUTINE FCN(N,X,FVEC,IFLAG)
C INTEGER N,IFLAG
C REAL X(N),FVEC(N)
C INTEGER K
C REAL ONE,TEMP,TEMP1,TEMP2,THREE,TWO,ZERO
C DATA ZERO,ONE,TWO,THREE /0.E0,1.E0,2.E0,3.E0/
C C
C DO 10 K = 1, N
C TEMP = (THREE - TWO*X(K))*X(K)
C TEMP1 = ZERO
C IF (K .NE. 1) TEMP1 = X(K-1)
C TEMP2 = ZERO
C IF (K .NE. N) TEMP2 = X(K+1)
C FVEC(K) = TEMP - TEMP1 - TWO*TEMP2 + ONE
C 10 CONTINUE
C RETURN
C END
C
C Results obtained with different compilers or machines
C may be slightly different.
C
C FINAL L2 NORM OF THE RESIDUALS 0.1192636E-07
C
C EXIT PARAMETER 1
C
C FINAL APPROXIMATE SOLUTION
C
C -0.5706545E+00 -0.6816283E+00 -0.7017325E+00
C -0.7042129E+00 -0.7013690E+00 -0.6918656E+00
C -0.6657920E+00 -0.5960342E+00 -0.4164121E+00
C***REFERENCES POWELL, M. J. D.
C A HYBRID METHOD FOR NONLINEAR EQUATIONS.
C NUMERICAL METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS,
C P. RABINOWITZ, EDITOR. GORDON AND BREACH, 1970.
C***ROUTINES CALLED SNSQ,XERROR
C***END PROLOGUE SNSQE
INTEGER IOPT,N,NPRINT,INFO,LWA
REAL TOL
REAL X(N),FVEC(N),WA(LWA)
EXTERNAL FCN,JAC
INTEGER INDEX,J,LR,MAXFEV,ML,MODE,MU,NFEV,NJEV
REAL EPSFCN,FACTOR,ONE,XTOL,ZERO
DATA FACTOR,ONE,ZERO /1.0E2,1.0E0,0.0E0/
C***FIRST EXECUTABLE STATEMENT SNSQE
INFO = 0
C
C CHECK THE INPUT PARAMETERS FOR ERRORS.
C
IF (IOPT .LT. 1 .OR. IOPT .GT. 2 .OR. N .LE. 0
1 .OR. TOL .LT. ZERO .OR. LWA .LT. (3*N**2 +13*N)/2)
2 GO TO 20
C
C CALL SNSQ.
C
MAXFEV = 100*(N + 1)
IF (IOPT .EQ. 2) MAXFEV = 2*MAXFEV
XTOL = TOL
ML = N - 1
MU = N - 1
EPSFCN = ZERO
MODE = 2
DO 10 J = 1, N
WA(J) = ONE
10 CONTINUE
LR = (N*(N + 1))/2
INDEX=6*N+LR
CALL SNSQ(FCN,JAC,IOPT,N,X,FVEC,WA(INDEX+1),N,XTOL,MAXFEV,ML,MU,
1 EPSFCN,WA(1),MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,
2 WA(6*N+1),LR,WA(N+1),WA(2*N+1),WA(3*N+1),WA(4*N+1),
3 WA(5*N+1))
IF (INFO .EQ. 5) INFO = 4
20 CONTINUE
IF (INFO .EQ. 0) CALL XERROR( 'SNSQE -- INVALID INPUT PARAMETER.'
1,34,2,1)
RETURN
C
C LAST CARD OF SUBROUTINE SNSQE.
C
END
SUBROUTINE DOGLEG(N,R,LR,DIAG,QTB,DELTA,X,WA1,WA2)
INTEGER N,LR
REAL DELTA
REAL R(LR),DIAG(N),QTB(N),X(N),WA1(N),WA2(N)
INTEGER I,J,JJ,JP1,K,L
REAL ALPHA,BNORM,EPSMCH,GNORM,ONE,QNORM,SGNORM,SUM,TEMP,ZERO
REAL R1MACH,SNRM2
DATA ONE,ZERO /1.0E0,0.0E0/
EPSMCH = R1MACH(4)
JJ = (N*(N + 1))/2 + 1
DO 50 K = 1, N
J = N - K + 1
JP1 = J + 1
JJ = JJ - K
L = JJ + 1
SUM = ZERO
IF (N .LT. JP1) GO TO 20
DO 10 I = JP1, N
SUM = SUM + R(L)*X(I)
L = L + 1
10 CONTINUE
20 CONTINUE
TEMP = R(JJ)
IF (TEMP .NE. ZERO) GO TO 40
L = J
DO 30 I = 1, J
TEMP = AMAX1(TEMP,ABS(R(L)))
L = L + N - I
30 CONTINUE
TEMP = EPSMCH*TEMP
IF (TEMP .EQ. ZERO) TEMP = EPSMCH
40 CONTINUE
X(J) = (QTB(J) - SUM)/TEMP
50 CONTINUE
DO 60 J = 1, N
WA1(J) = ZERO
WA2(J) = DIAG(J)*X(J)
60 CONTINUE
QNORM = SNRM2(N,WA2,1)
IF (QNORM .LE. DELTA) GO TO 140
L = 1
DO 80 J = 1, N
TEMP = QTB(J)
DO 70 I = J, N
WA1(I) = WA1(I) + R(L)*TEMP
L = L + 1
70 CONTINUE
WA1(J) = WA1(J)/DIAG(J)
80 CONTINUE
GNORM = SNRM2(N,WA1,1)
SGNORM = ZERO
ALPHA = DELTA/QNORM
IF (GNORM .EQ. ZERO) GO TO 120
DO 90 J = 1, N
WA1(J) = (WA1(J)/GNORM)/DIAG(J)
90 CONTINUE
L = 1
DO 110 J = 1, N
SUM = ZERO
DO 100 I = J, N
SUM = SUM + R(L)*WA1(I)
L = L + 1
100 CONTINUE
WA2(J) = SUM
110 CONTINUE
TEMP = SNRM2(N,WA2,1)
SGNORM = (GNORM/TEMP)/TEMP
ALPHA = ZERO
IF (SGNORM .GE. DELTA) GO TO 120
BNORM = SNRM2(N,QTB,1)
TEMP = (BNORM/GNORM)*(BNORM/QNORM)*(SGNORM/DELTA)
TEMP = TEMP - (DELTA/QNORM)*(SGNORM/DELTA)**2
1 + SQRT((TEMP-(DELTA/QNORM))**2
2 +(ONE-(DELTA/QNORM)**2)*(ONE-(SGNORM/DELTA)**2))
ALPHA = ((DELTA/QNORM)*(ONE - (SGNORM/DELTA)**2))/TEMP
120 CONTINUE
TEMP = (ONE - ALPHA)*AMIN1(SGNORM,DELTA)
DO 130 J = 1, N
X(J) = TEMP*WA1(J) + ALPHA*X(J)
130 CONTINUE
140 CONTINUE
RETURN
END
SUBROUTINE FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN,
1 WA1,WA2)
INTEGER N,LDFJAC,IFLAG,ML,MU
REAL EPSFCN
REAL X(N),FVEC(N),FJAC(LDFJAC,N),WA1(N),WA2(N)
INTEGER I,J,K,MSUM
REAL EPS,EPSMCH,H,TEMP,ZERO
REAL R1MACH
DATA ZERO /0.0E0/
EPSMCH = R1MACH(4)
EPS = SQRT(AMAX1(EPSFCN,EPSMCH))
MSUM = ML + MU + 1
IF (MSUM .LT. N) GO TO 40
DO 20 J = 1, N
TEMP = X(J)
H = EPS*ABS(TEMP)
IF (H .EQ. ZERO) H = EPS
X(J) = TEMP + H
CALL FCN(N,X,WA1,IFLAG)
IF (IFLAG .LT. 0) GO TO 30
X(J) = TEMP
DO 10 I = 1, N
FJAC(I,J) = (WA1(I) - FVEC(I))/H
10 CONTINUE
20 CONTINUE
30 CONTINUE
GO TO 110
40 CONTINUE
DO 90 K = 1, MSUM
DO 60 J = K, N, MSUM
WA2(J) = X(J)
H = EPS*ABS(WA2(J))
IF (H .EQ. ZERO) H = EPS
X(J) = WA2(J) + H
60 CONTINUE
CALL FCN(N,X,WA1,IFLAG)
IF (IFLAG .LT. 0) GO TO 100
DO 80 J = K, N, MSUM
X(J) = WA2(J)
H = EPS*ABS(WA2(J))
IF (H .EQ. ZERO) H = EPS
DO 70 I = 1, N
FJAC(I,J) = ZERO
IF (I .GE. J - MU .AND. I .LE. J + ML)
1 FJAC(I,J) = (WA1(I) - FVEC(I))/H
70 CONTINUE
80 CONTINUE
90 CONTINUE
100 CONTINUE
110 CONTINUE
RETURN
END
SUBROUTINE QFORM(M,N,Q,LDQ,WA)
INTEGER M,N,LDQ
REAL Q(LDQ,M),WA(M)
INTEGER I,J,JM1,K,L,MINMN,NP1
REAL ONE,SUM,TEMP,ZERO
DATA ONE,ZERO /1.0E0,0.0E0/
MINMN = MIN0(M,N)
IF (MINMN .LT. 2) GO TO 30
DO 20 J = 2, MINMN
JM1 = J - 1
DO 10 I = 1, JM1
Q(I,J) = ZERO
10 CONTINUE
20 CONTINUE
30 CONTINUE
NP1 = N + 1
IF (M .LT. NP1) GO TO 60
DO 50 J = NP1, M
DO 40 I = 1, M
Q(I,J) = ZERO
40 CONTINUE
Q(J,J) = ONE
50 CONTINUE
60 CONTINUE
DO 120 L = 1, MINMN
K = MINMN - L + 1
DO 70 I = K, M
WA(I) = Q(I,K)
Q(I,K) = ZERO
70 CONTINUE
Q(K,K) = ONE
IF (WA(K) .EQ. ZERO) GO TO 110
DO 100 J = K, M
SUM = ZERO
DO 80 I = K, M
SUM = SUM + Q(I,J)*WA(I)
80 CONTINUE
TEMP = SUM/WA(K)
DO 90 I = K, M
Q(I,J) = Q(I,J) - TEMP*WA(I)
90 CONTINUE
100 CONTINUE
110 CONTINUE
120 CONTINUE
RETURN
END
SUBROUTINE QRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,SIGMA,ACNORM,WA)
INTEGER M,N,LDA,LIPVT
INTEGER IPVT(LIPVT)
LOGICAL PIVOT
REAL A(LDA,N),SIGMA(N),ACNORM(N),WA(N)
INTEGER I,J,JP1,K,KMAX,MINMN
REAL AJNORM,EPSMCH,ONE,P05,SUM,TEMP,ZERO
REAL R1MACH,SNRM2
DATA ONE,P05,ZERO /1.0E0,5.0E-2,0.0E0/
EPSMCH = R1MACH(4)
DO 10 J = 1, N
ACNORM(J) = SNRM2(M,A(1,J),1)
SIGMA(J) = ACNORM(J)
WA(J) = SIGMA(J)
IF (PIVOT) IPVT(J) = J
10 CONTINUE
MINMN = MIN0(M,N)
DO 110 J = 1, MINMN
IF (.NOT.PIVOT) GO TO 40
KMAX = J
DO 20 K = J, N
IF (SIGMA(K) .GT. SIGMA(KMAX)) KMAX = K
20 CONTINUE
IF (KMAX .EQ. J) GO TO 40
DO 30 I = 1, M
TEMP = A(I,J)
A(I,J) = A(I,KMAX)
A(I,KMAX) = TEMP
30 CONTINUE
SIGMA(KMAX) = SIGMA(J)
WA(KMAX) = WA(J)
K = IPVT(J)
IPVT(J) = IPVT(KMAX)
IPVT(KMAX) = K
40 CONTINUE
AJNORM = SNRM2(M-J+1,A(J,J),1)
IF (AJNORM .EQ. ZERO) GO TO 100
IF (A(J,J) .LT. ZERO) AJNORM = -AJNORM
DO 50 I = J, M
A(I,J) = A(I,J)/AJNORM
50 CONTINUE
A(J,J) = A(J,J) + ONE
JP1 = J + 1
IF (N .LT. JP1) GO TO 100
DO 90 K = JP1, N
SUM = ZERO
DO 60 I = J, M
SUM = SUM + A(I,J)*A(I,K)
60 CONTINUE
TEMP = SUM/A(J,J)
DO 70 I = J, M
A(I,K) = A(I,K) - TEMP*A(I,J)
70 CONTINUE
IF (.NOT.PIVOT .OR. SIGMA(K) .EQ. ZERO) GO TO 80
TEMP = A(J,K)/SIGMA(K)
SIGMA(K) = SIGMA(K)*SQRT(AMAX1(ZERO,ONE-TEMP**2))
IF (P05*(SIGMA(K)/WA(K))**2 .GT. EPSMCH) GO TO 80
SIGMA(K) = SNRM2(M-J,A(JP1,K),1)
WA(K) = SIGMA(K)
80 CONTINUE
90 CONTINUE
100 CONTINUE
SIGMA(J) = -AJNORM
110 CONTINUE
RETURN
END
SUBROUTINE R1MPYQ(M,N,A,LDA,V,W)
INTEGER M,N,LDA
REAL A(LDA,N),V(N),W(N)
INTEGER I,J,NMJ,NM1
REAL COS,ONE,SIN,TEMP
DATA ONE /1.0E0/
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 50
DO 20 NMJ = 1, NM1
J = N - NMJ
IF (ABS(V(J)) .GT. ONE) COS = ONE/V(J)
IF (ABS(V(J)) .GT. ONE) SIN = SQRT(ONE-COS**2)
IF (ABS(V(J)) .LE. ONE) SIN = V(J)
IF (ABS(V(J)) .LE. ONE) COS = SQRT(ONE-SIN**2)
DO 10 I = 1, M
TEMP = COS*A(I,J) - SIN*A(I,N)
A(I,N) = SIN*A(I,J) + COS*A(I,N)
A(I,J) = TEMP
10 CONTINUE
20 CONTINUE
DO 40 J = 1, NM1
IF (ABS(W(J)) .GT. ONE) COS = ONE/W(J)
IF (ABS(W(J)) .GT. ONE) SIN = SQRT(ONE-COS**2)
IF (ABS(W(J)) .LE. ONE) SIN = W(J)
IF (ABS(W(J)) .LE. ONE) COS = SQRT(ONE-SIN**2)
DO 30 I = 1, M
TEMP = COS*A(I,J) + SIN*A(I,N)
A(I,N) = -SIN*A(I,J) + COS*A(I,N)
A(I,J) = TEMP
30 CONTINUE
40 CONTINUE
50 CONTINUE
RETURN
END
SUBROUTINE R1UPDT(M,N,S,LS,U,V,W,SING)
INTEGER M,N,LS
LOGICAL SING
REAL S(LS),U(M),V(N),W(M)
INTEGER I,J,JJ,L,NMJ,NM1
REAL COS,COTAN,GIANT,ONE,P5,P25,SIN,TAN,TAU,TEMP,ZERO
REAL R1MACH
DATA ONE,P5,P25,ZERO /1.0E0,5.0E-1,2.5E-1,0.0E0/
GIANT = R1MACH(2)
JJ = (N*(2*M - N + 1))/2 - (M - N)
L = JJ
DO 10 I = N, M
W(I) = S(L)
L = L + 1
10 CONTINUE
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 70
DO 60 NMJ = 1, NM1
J = N - NMJ
JJ = JJ - (M - J + 1)
W(J) = ZERO
IF (V(J) .EQ. ZERO) GO TO 50
IF (ABS(V(N)) .GE. ABS(V(J))) GO TO 20
COTAN = V(N)/V(J)
SIN = P5/SQRT(P25+P25*COTAN**2)
COS = SIN*COTAN
TAU = ONE
IF (ABS(COS)*GIANT .GT. ONE) TAU = ONE/COS
GO TO 30
20 CONTINUE
TAN = V(J)/V(N)
COS = P5/SQRT(P25+P25*TAN**2)
SIN = COS*TAN
TAU = SIN
30 CONTINUE
V(N) = SIN*V(J) + COS*V(N)
V(J) = TAU
L = JJ
DO 40 I = J, M
TEMP = COS*S(L) - SIN*W(I)
W(I) = SIN*S(L) + COS*W(I)
S(L) = TEMP
L = L + 1
40 CONTINUE
50 CONTINUE
60 CONTINUE
70 CONTINUE
DO 80 I = 1, M
W(I) = W(I) + V(N)*U(I)
80 CONTINUE
SING = .FALSE.
IF (NM1 .LT. 1) GO TO 140
DO 130 J = 1, NM1
IF (W(J) .EQ. ZERO) GO TO 120
IF (ABS(S(JJ)) .GE. ABS(W(J))) GO TO 90
COTAN = S(JJ)/W(J)
SIN = P5/SQRT(P25+P25*COTAN**2)
COS = SIN*COTAN
TAU = ONE
IF (ABS(COS)*GIANT .GT. ONE) TAU = ONE/COS
GO TO 100
90 CONTINUE
TAN = W(J)/S(JJ)
COS = P5/SQRT(P25+P25*TAN**2)
SIN = COS*TAN
TAU = SIN
100 CONTINUE
L = JJ
DO 110 I = J, M
TEMP = COS*S(L) + SIN*W(I)
W(I) = -SIN*S(L) + COS*W(I)
S(L) = TEMP
L = L + 1
110 CONTINUE
W(J) = TAU
120 CONTINUE
IF (S(JJ) .EQ. ZERO) SING = .TRUE.
JJ = JJ + (M - J + 1)
130 CONTINUE
140 CONTINUE
L = JJ
DO 150 I = N, M
S(L) = W(I)
L = L + 1
150 CONTINUE
IF (S(JJ) .EQ. ZERO) SING = .TRUE.
RETURN
END
SUBROUTINE SNSQ(FCN,JAC,IOPT,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,ML,
1 MU,EPSFCN,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF,WA1,
2 WA2,WA3,WA4)
INTEGER IOPT,N,MAXFEV,ML,MU,MODE,NPRINT,INFO,NFEV,LDFJAC,LR,NJEV
REAL XTOL,EPSFCN,FACTOR
REAL X(N),FVEC(N),DIAG(N),FJAC(LDFJAC,N),R(LR),QTF(N),WA1(N),
1 WA2(N),WA3(N),WA4(N)
EXTERNAL FCN
INTEGER I,IFLAG,ITER,J,JM1,L,NCFAIL,NCSUC,NSLOW1,NSLOW2
INTEGER IWA(1)
LOGICAL JEVAL,SING
REAL ACTRED,DELTA,EPSMCH,FNORM,FNORM1,ONE,PNORM,PRERED,P1,P5,
1 P001,P0001,RATIO,SUM,TEMP,XNORM,ZERO
REAL R1MACH,SNRM2
DATA ONE,P1,P5,P001,P0001,ZERO
1 /1.0E0,1.0E-1,5.0E-1,1.0E-3,1.0E-4,0.0E0/
EPSMCH = R1MACH(4)
INFO = 0
IFLAG = 0
NFEV = 0
NJEV = 0
IF (IOPT .LT. 1 .OR. IOPT .GT. 2 .OR.
1 N .LE. 0 .OR. XTOL .LT. ZERO .OR. MAXFEV .LE. 0
2 .OR. ML .LT. 0 .OR. MU .LT. 0 .OR. FACTOR .LE. ZERO
3 .OR. LDFJAC .LT. N .OR. LR .LT. (N*(N + 1))/2) 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
IFLAG = 1
CALL FCN(N,X,FVEC,IFLAG)
NFEV = 1
IF (IFLAG .LT. 0) GO TO 300
FNORM = SNRM2(N,FVEC,1)
ITER = 1
NCSUC = 0
NCFAIL = 0
NSLOW1 = 0
NSLOW2 = 0
30 CONTINUE
JEVAL = .TRUE.
IF (IOPT .EQ. 2) GO TO 31
CALL JAC(N,X,FVEC,FJAC,LDFJAC,IFLAG)
NJEV = NJEV+1
GO TO 32
31 IFLAG = 2
CALL FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN,WA1,
1 WA2)
NFEV = NFEV + MIN0(ML+MU+1,N)
32 IF (IFLAG .LT. 0) GO TO 300
CALL QRFAC(N,N,FJAC,LDFJAC,.FALSE.,IWA,1,WA1,WA2,WA3)
IF (ITER .NE. 1) GO TO 70
IF (MODE .EQ. 2) GO TO 50
DO 40 J = 1, N
DIAG(J) = WA2(J)
IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE
40 CONTINUE
50 CONTINUE
DO 60 J = 1, N
WA3(J) = DIAG(J)*X(J)
60 CONTINUE
XNORM = SNRM2(N,WA3,1)
DELTA = FACTOR*XNORM
IF (DELTA .EQ. ZERO) DELTA = FACTOR
70 CONTINUE
DO 80 I = 1, N
QTF(I) = FVEC(I)
80 CONTINUE
DO 120 J = 1, N
IF (FJAC(J,J) .EQ. ZERO) GO TO 110
SUM = ZERO
DO 90 I = J, N
SUM = SUM + FJAC(I,J)*QTF(I)
90 CONTINUE
TEMP = -SUM/FJAC(J,J)
DO 100 I = J, N
QTF(I) = QTF(I) + FJAC(I,J)*TEMP
100 CONTINUE
110 CONTINUE
120 CONTINUE
SING = .FALSE.
DO 150 J = 1, N
L = J
JM1 = J - 1
IF (JM1 .LT. 1) GO TO 140
DO 130 I = 1, JM1
R(L) = FJAC(I,J)
L = L + N - I
130 CONTINUE
140 CONTINUE
R(L) = WA1(J)
IF (WA1(J) .EQ. ZERO) SING = .TRUE.
150 CONTINUE
CALL QFORM(N,N,FJAC,LDFJAC,WA1)
IF (MODE .EQ. 2) GO TO 170
DO 160 J = 1, N
DIAG(J) = AMAX1(DIAG(J),WA2(J))
160 CONTINUE
170 CONTINUE
180 CONTINUE
IF (NPRINT .LE. 0) GO TO 190
IFLAG = 0
IF (MOD(ITER-1,NPRINT) .EQ. 0) CALL FCN(N,X,FVEC,IFLAG)
IF (IFLAG .LT. 0) GO TO 300
190 CONTINUE
CALL DOGLEG(N,R,LR,DIAG,QTF,DELTA,WA1,WA2,WA3)
DO 200 J = 1, N
WA1(J) = -WA1(J)
WA2(J) = X(J) + WA1(J)
WA3(J) = DIAG(J)*WA1(J)
200 CONTINUE
PNORM = SNRM2(N,WA3,1)
IF (ITER .EQ. 1) DELTA = AMIN1(DELTA,PNORM)
IFLAG = 1
CALL FCN(N,WA2,WA4,IFLAG)
NFEV = NFEV + 1
IF (IFLAG .LT. 0) GO TO 300
FNORM1 = SNRM2(N,WA4,1)
ACTRED = -ONE
IF (FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2
L = 1
DO 220 I = 1, N
SUM = ZERO
DO 210 J = I, N
SUM = SUM + R(L)*WA1(J)
L = L + 1
210 CONTINUE
WA3(I) = QTF(I) + SUM
220 CONTINUE
TEMP = SNRM2(N,WA3,1)
PRERED = ZERO
IF (TEMP .LT. FNORM) PRERED = ONE - (TEMP/FNORM)**2
RATIO = ZERO
IF (PRERED .GT. ZERO) RATIO = ACTRED/PRERED
IF (RATIO .GE. P1) GO TO 230
NCSUC = 0
NCFAIL = NCFAIL + 1
DELTA = P5*DELTA
GO TO 240
230 CONTINUE
NCFAIL = 0
NCSUC = NCSUC + 1
IF (RATIO .GE. P5 .OR. NCSUC .GT. 1)
1 DELTA = AMAX1(DELTA,PNORM/P5)
IF (ABS(RATIO-ONE) .LE. P1) DELTA = PNORM/P5
240 CONTINUE
IF (RATIO .LT. P0001) GO TO 260
DO 250 J = 1, N
X(J) = WA2(J)
WA2(J) = DIAG(J)*X(J)
FVEC(J) = WA4(J)
250 CONTINUE
XNORM = SNRM2(N,WA2,1)
FNORM = FNORM1
ITER = ITER + 1
260 CONTINUE
NSLOW1 = NSLOW1 + 1
IF (ACTRED .GE. P001) NSLOW1 = 0
IF (JEVAL) NSLOW2 = NSLOW2 + 1
IF (ACTRED .GE. P1) NSLOW2 = 0
IF (DELTA .LE. XTOL*XNORM .OR. FNORM .EQ. ZERO) INFO = 1
IF (INFO .NE. 0) GO TO 300
IF (NFEV .GE. MAXFEV) INFO = 2
IF (P1*AMAX1(P1*DELTA,PNORM) .LE. EPSMCH*XNORM) INFO = 3
IF (NSLOW2 .EQ. 5) INFO = 4
IF (NSLOW1 .EQ. 10) INFO = 5
IF (INFO .NE. 0) GO TO 300
IF (NCFAIL .EQ. 2) GO TO 290
DO 280 J = 1, N
SUM = ZERO
DO 270 I = 1, N
SUM = SUM + FJAC(I,J)*WA4(I)
270 CONTINUE
WA2(J) = (SUM - WA3(J))/PNORM
WA1(J) = DIAG(J)*((DIAG(J)*WA1(J))/PNORM)
IF (RATIO .GE. P0001) QTF(J) = SUM
280 CONTINUE
CALL R1UPDT(N,N,R,LR,WA1,WA2,WA3,SING)
CALL R1MPYQ(N,N,FJAC,LDFJAC,WA2,WA3)
CALL R1MPYQ(1,N,QTF,1,WA2,WA3)
JEVAL = .FALSE.
GO TO 180
290 CONTINUE
GO TO 30
300 CONTINUE
IF (IFLAG .LT. 0) INFO = IFLAG
IFLAG = 0
IF (NPRINT .GT. 0) CALL FCN(N,X,FVEC,IFLAG)
IF (INFO .LT. 0) CALL XERROR( 'SNSQ -- EXECUTION TERMINATED BECA
1USE USER SET IFLAG NEGATIVE.',63,1,1)
IF (INFO .EQ. 0) CALL XERROR( 'SNSQ -- INVALID INPUT PARAMETER.'
1 ,34,2,1)
IF (INFO .EQ. 2) CALL XERROR( 'SNSQ -- TOO MANY FUNCTION EVALUAT
1IONS.',40,9,1)
IF (INFO .EQ. 3) CALL XERROR( 'SNSQ -- XTOL TOO SMALL. NO FURTHE
1R IMPROVEMENT POSSIBLE.',58,3,1)
IF (INFO .GT. 4) CALL XERROR( 'SNSQ -- ITERATION NOT MAKING GOOD
1 PROGRESS.',45,1,1)
RETURN
END