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