1*DECK SNSQE 2 SUBROUTINE SNSQE (FCN, JAC, IOPT, N, X, FVEC, TOL, NPRINT, INFO, 3 + WA, LWA) 4C***BEGIN PROLOGUE SNSQE 5C***PURPOSE An easy-to-use code to find a zero of a system of N 6C nonlinear functions in N variables by a modification of 7C the Powell hybrid method. 8C***LIBRARY SLATEC 9C***CATEGORY F2A 10C***TYPE SINGLE PRECISION (SNSQE-S, DNSQE-D) 11C***KEYWORDS EASY-TO-USE, NONLINEAR SQUARE SYSTEM, 12C POWELL HYBRID METHOD, ZEROS 13C***AUTHOR Hiebert, K. L., (SNLA) 14C***DESCRIPTION 15C 16C 1. Purpose. 17C 18C 19C The purpose of SNSQE is to find a zero of a system of N non- 20C linear functions in N variables by a modification of the Powell 21C hybrid method. This is done by using the more general nonlinear 22C equation solver SNSQ. The user must provide a subroutine which 23C calculates the functions. The user has the option of either to 24C provide a subroutine which calculates the Jacobian or to let the 25C code calculate it by a forward-difference approximation. This 26C code is the combination of the MINPACK codes (Argonne) HYBRD1 27C and HYBRJ1. 28C 29C 30C 2. Subroutine and Type Statements. 31C 32C SUBROUTINE SNSQE(FCN,JAC,IOPT,N,X,FVEC,TOL,NPRINT,INFO, 33C * WA,LWA) 34C INTEGER IOPT,N,NPRINT,INFO,LWA 35C REAL TOL 36C REAL X(N),FVEC(N),WA(LWA) 37C EXTERNAL FCN,JAC 38C 39C 40C 3. Parameters. 41C 42C Parameters designated as input parameters must be specified on 43C entry to SNSQE and are not changed on exit, while parameters 44C designated as output parameters need not be specified on entry 45C and are set to appropriate values on exit from SNSQE. 46C 47C FCN is the name of the user-supplied subroutine which calculates 48C the functions. FCN must be declared in an EXTERNAL statement 49C in the user calling program, and should be written as follows. 50C 51C SUBROUTINE FCN(N,X,FVEC,IFLAG) 52C INTEGER N,IFLAG 53C REAL X(N),FVEC(N) 54C ---------- 55C Calculate the functions at X and 56C return this vector in FVEC. 57C ---------- 58C RETURN 59C END 60C 61C The value of IFLAG should not be changed by FCN unless the 62C user wants to terminate execution of SNSQE. In this case, set 63C IFLAG to a negative integer. 64C 65C JAC is the name of the user-supplied subroutine which calculates 66C the Jacobian. If IOPT=1, then JAC must be declared in an 67C EXTERNAL statement in the user calling program, and should be 68C written as follows. 69C 70C SUBROUTINE JAC(N,X,FVEC,FJAC,LDFJAC,IFLAG) 71C INTEGER N,LDFJAC,IFLAG 72C REAL X(N),FVEC(N),FJAC(LDFJAC,N) 73C ---------- 74C Calculate the Jacobian at X and return this 75C matrix in FJAC. FVEC contains the function 76C values at X and should not be altered. 77C ---------- 78C RETURN 79C END 80C 81C The value of IFLAG should not be changed by JAC unless the 82C user wants to terminate execution of SNSQE. In this case, set 83C IFLAG to a negative integer. 84C 85C If IOPT=2, JAC can be ignored (treat it as a dummy argument). 86C 87C IOPT is an input variable which specifies how the Jacobian will 88C be calculated. If IOPT=1, then the user must supply the 89C Jacobian through the subroutine JAC. If IOPT=2, then the 90C code will approximate the Jacobian by forward-differencing. 91C 92C N is a positive integer input variable set to the number of 93C functions and variables. 94C 95C X is an array of length N. On input, X must contain an initial 96C estimate of the solution vector. On output, X contains the 97C final estimate of the solution vector. 98C 99C FVEC is an output array of length N which contains the functions 100C evaluated at the output X. 101C 102C TOL is a non-negative input variable. Termination occurs when 103C the algorithm estimates that the relative error between X and 104C the solution is at most TOL. Section 4 contains more details 105C about TOL. 106C 107C NPRINT is an integer input variable that enables controlled 108C printing of iterates if it is positive. In this case, FCN is 109C called with IFLAG = 0 at the beginning of the first iteration 110C and every NPRINT iteration thereafter and immediately prior 111C to return, with X and FVEC available for printing. Appropriate 112C print statements must be added to FCN (see example). If NPRINT 113C is not positive, no special calls of FCN with IFLAG = 0 are 114C made. 115C 116C INFO is an integer output variable. If the user has terminated 117C execution, INFO is set to the (negative) value of IFLAG. See 118C description of FCN and JAC. Otherwise, INFO is set as follows. 119C 120C INFO = 0 improper input parameters. 121C 122C INFO = 1 algorithm estimates that the relative error between 123C X and the solution is at most TOL. 124C 125C INFO = 2 number of calls to FCN has reached or exceeded 126C 100*(N+1) for IOPT=1 or 200*(N+1) for IOPT=2. 127C 128C INFO = 3 TOL is too small. No further improvement in the 129C approximate solution X is possible. 130C 131C INFO = 4 iteration is not making good progress. 132C 133C Sections 4 and 5 contain more details about INFO. 134C 135C WA is a work array of length LWA. 136C 137C LWA is a positive integer input variable not less than 138C (3*N**2+13*N))/2. 139C 140C 141C 4. Successful Completion. 142C 143C The accuracy of SNSQE is controlled by the convergence parame- 144C ter TOL. This parameter is used in a test which makes a compar- 145C ison between the approximation X and a solution XSOL. SNSQE 146C terminates when the test is satisfied. If TOL is less than the 147C machine precision (as defined by the function R1MACH(4)), then 148C SNSQE attempts only to satisfy the test defined by the machine 149C precision. Further progress is not usually possible. Unless 150C high precision solutions are required, the recommended value 151C for TOL is the square root of the machine precision. 152C 153C The test assumes that the functions are reasonably well behaved, 154C and, if the Jacobian is supplied by the user, that the functions 155C and the Jacobian coded consistently. If these conditions 156C are not satisfied, SNSQE may incorrectly indicate convergence. 157C The coding of the Jacobian can be checked by the subroutine 158C CHKDER. If the Jacobian is coded correctly or IOPT=2, then 159C the validity of the answer can be checked, for example, by 160C rerunning SNSQE with a tighter tolerance. 161C 162C Convergence Test. If ENORM(Z) denotes the Euclidean norm of a 163C vector Z, then this test attempts to guarantee that 164C 165C ENORM(X-XSOL) .LE. TOL*ENORM(XSOL). 166C 167C If this condition is satisfied with TOL = 10**(-K), then the 168C larger components of X have K significant decimal digits and 169C INFO is set to 1. There is a danger that the smaller compo- 170C nents of X may have large relative errors, but the fast rate 171C of convergence of SNSQE usually avoids this possibility. 172C 173C 174C 5. Unsuccessful Completion. 175C 176C Unsuccessful termination of SNSQE can be due to improper input 177C parameters, arithmetic interrupts, an excessive number of func- 178C tion evaluations, errors in the functions, or lack of good prog- 179C ress. 180C 181C Improper Input Parameters. INFO is set to 0 if IOPT .LT. 1, or 182C IOPT .GT. 2, or N .LE. 0, or TOL .LT. 0.E0, or 183C LWA .LT. (3*N**2+13*N)/2. 184C 185C Arithmetic Interrupts. If these interrupts occur in the FCN 186C subroutine during an early stage of the computation, they may 187C be caused by an unacceptable choice of X by SNSQE. In this 188C case, it may be possible to remedy the situation by not evalu- 189C ating the functions here, but instead setting the components 190C of FVEC to numbers that exceed those in the initial FVEC. 191C 192C Excessive Number of Function Evaluations. If the number of 193C calls to FCN reaches 100*(N+1) for IOPT=1 or 200*(N+1) for 194C IOPT=2, then this indicates that the routine is converging 195C very slowly as measured by the progress of FVEC, and INFO is 196C set to 2. This situation should be unusual because, as 197C indicated below, lack of good progress is usually diagnosed 198C earlier by SNSQE, causing termination with INFO = 4. 199C 200C Errors in the Functions. When IOPT=2, the choice of step length 201C in the forward-difference approximation to the Jacobian 202C assumes that the relative errors in the functions are of the 203C order of the machine precision. If this is not the case, 204C SNSQE may fail (usually with INFO = 4). The user should 205C then either use SNSQ and set the step length or use IOPT=1 206C and supply the Jacobian. 207C 208C Lack of Good Progress. SNSQE searches for a zero of the system 209C by minimizing the sum of the squares of the functions. In so 210C doing, it can become trapped in a region where the minimum 211C does not correspond to a zero of the system and, in this situ- 212C ation, the iteration eventually fails to make good progress. 213C In particular, this will happen if the system does not have a 214C zero. If the system has a zero, rerunning SNSQE from a dif- 215C ferent starting point may be helpful. 216C 217C 218C 6. Characteristics of the Algorithm. 219C 220C SNSQE is a modification of the Powell hybrid method. Two of 221C its main characteristics involve the choice of the correction as 222C a convex combination of the Newton and scaled gradient direc- 223C tions, and the updating of the Jacobian by the rank-1 method of 224C Broyden. The choice of the correction guarantees (under reason- 225C able conditions) global convergence for starting points far from 226C the solution and a fast rate of convergence. The Jacobian is 227C calculated at the starting point by either the user-supplied 228C subroutine or a forward-difference approximation, but it is not 229C recalculated until the rank-1 method fails to produce satis- 230C factory progress. 231C 232C Timing. The time required by SNSQE to solve a given problem 233C depends on N, the behavior of the functions, the accuracy 234C requested, and the starting point. The number of arithmetic 235C operations needed by SNSQE is about 11.5*(N**2) to process 236C each evaluation of the functions (call to FCN) and 1.3*(N**3) 237C to process each evaluation of the Jacobian (call to JAC, 238C if IOPT = 1). Unless FCN and JAC can be evaluated quickly, 239C the timing of SNSQE will be strongly influenced by the time 240C spent in FCN and JAC. 241C 242C Storage. SNSQE requires (3*N**2 + 17*N)/2 single precision 243C storage locations, in addition to the storage required by the 244C program. There are no internally declared storage arrays. 245C 246C 247C 7. Example. 248C 249C The problem is to determine the values of X(1), X(2), ..., X(9), 250C which solve the system of tridiagonal equations 251C 252C (3-2*X(1))*X(1) -2*X(2) = -1 253C -X(I-1) + (3-2*X(I))*X(I) -2*X(I+1) = -1, I=2-8 254C -X(8) + (3-2*X(9))*X(9) = -1 255C 256C ********** 257C 258C PROGRAM TEST 259C C 260C C Driver for SNSQE example. 261C C 262C INTEGER J,N,IOPT,NPRINT,INFO,LWA,NWRITE 263C REAL TOL,FNORM 264C REAL X(9),FVEC(9),WA(180) 265C REAL ENORM,R1MACH 266C EXTERNAL FCN 267C DATA NWRITE /6/ 268C C 269C IOPT = 2 270C N = 9 271C C 272C C The following starting values provide a rough solution. 273C C 274C DO 10 J = 1, 9 275C X(J) = -1.E0 276C 10 CONTINUE 277C 278C LWA = 180 279C NPRINT = 0 280C C 281C C Set TOL to the square root of the machine precision. 282C C Unless high precision solutions are required, 283C C this is the recommended setting. 284C C 285C TOL = SQRT(R1MACH(4)) 286C C 287C CALL SNSQE(FCN,JAC,IOPT,N,X,FVEC,TOL,NPRINT,INFO,WA,LWA) 288C FNORM = ENORM(N,FVEC) 289C WRITE (NWRITE,1000) FNORM,INFO,(X(J),J=1,N) 290C STOP 291C 1000 FORMAT (5X,' FINAL L2 NORM OF THE RESIDUALS',E15.7 // 292C * 5X,' EXIT PARAMETER',16X,I10 // 293C * 5X,' FINAL APPROXIMATE SOLUTION' // (5X,3E15.7)) 294C END 295C SUBROUTINE FCN(N,X,FVEC,IFLAG) 296C INTEGER N,IFLAG 297C REAL X(N),FVEC(N) 298C INTEGER K 299C REAL ONE,TEMP,TEMP1,TEMP2,THREE,TWO,ZERO 300C DATA ZERO,ONE,TWO,THREE /0.E0,1.E0,2.E0,3.E0/ 301C C 302C DO 10 K = 1, N 303C TEMP = (THREE - TWO*X(K))*X(K) 304C TEMP1 = ZERO 305C IF (K .NE. 1) TEMP1 = X(K-1) 306C TEMP2 = ZERO 307C IF (K .NE. N) TEMP2 = X(K+1) 308C FVEC(K) = TEMP - TEMP1 - TWO*TEMP2 + ONE 309C 10 CONTINUE 310C RETURN 311C END 312C 313C Results obtained with different compilers or machines 314C may be slightly different. 315C 316C FINAL L2 NORM OF THE RESIDUALS 0.1192636E-07 317C 318C EXIT PARAMETER 1 319C 320C FINAL APPROXIMATE SOLUTION 321C 322C -0.5706545E+00 -0.6816283E+00 -0.7017325E+00 323C -0.7042129E+00 -0.7013690E+00 -0.6918656E+00 324C -0.6657920E+00 -0.5960342E+00 -0.4164121E+00 325C 326C***REFERENCES M. J. D. Powell, A hybrid method for nonlinear equa- 327C tions. In Numerical Methods for Nonlinear Algebraic 328C Equations, P. Rabinowitz, Editor. Gordon and Breach, 329C 1988. 330C***ROUTINES CALLED SNSQ, XERMSG 331C***REVISION HISTORY (YYMMDD) 332C 800301 DATE WRITTEN 333C 890831 Modified array declarations. (WRB) 334C 890831 REVISION DATE from Version 3.2 335C 891214 Prologue converted to Version 4.0 format. (BAB) 336C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 337C 920501 Reformatted the REFERENCES section. (WRB) 338C***END PROLOGUE SNSQE 339 INTEGER IOPT,N,NPRINT,INFO,LWA 340 REAL TOL 341 REAL X(*),FVEC(*),WA(LWA) 342 EXTERNAL FCN, JAC 343 INTEGER INDEX,J,LR,MAXFEV,ML,MODE,MU,NFEV,NJEV 344 REAL EPSFCN,FACTOR,ONE,XTOL,ZERO 345 SAVE FACTOR, ONE, ZERO 346 DATA FACTOR,ONE,ZERO /1.0E2,1.0E0,0.0E0/ 347C***FIRST EXECUTABLE STATEMENT SNSQE 348 INFO = 0 349C 350C CHECK THE INPUT PARAMETERS FOR ERRORS. 351C 352 IF (IOPT .LT. 1 .OR. IOPT .GT. 2 .OR. N .LE. 0 353 1 .OR. TOL .LT. ZERO .OR. LWA .LT. (3*N**2 +13*N)/2) 354 2 GO TO 20 355C 356C CALL SNSQ. 357C 358 MAXFEV = 100*(N + 1) 359 IF (IOPT .EQ. 2) MAXFEV = 2*MAXFEV 360 XTOL = TOL 361 ML = N - 1 362 MU = N - 1 363 EPSFCN = ZERO 364 MODE = 2 365 DO 10 J = 1, N 366 WA(J) = ONE 367 10 CONTINUE 368 LR = (N*(N + 1))/2 369 INDEX=6*N+LR 370 CALL SNSQ(FCN,JAC,IOPT,N,X,FVEC,WA(INDEX+1),N,XTOL,MAXFEV,ML,MU, 371 1 EPSFCN,WA(1),MODE,FACTOR,NPRINT,INFO,NFEV,NJEV, 372 2 WA(6*N+1),LR,WA(N+1),WA(2*N+1),WA(3*N+1),WA(4*N+1), 373 3 WA(5*N+1)) 374 IF (INFO .EQ. 5) INFO = 4 375 20 CONTINUE 376 IF (INFO .EQ. 0) CALL XERMSG ('SLATEC', 'SNSQE', 377 + 'INVALID INPUT PARAMETER.', 2, 1) 378 RETURN 379C 380C LAST CARD OF SUBROUTINE SNSQE. 381C 382 END 383