*DECK DSOS SUBROUTINE DSOS (FNC, NEQ, X, RTOLX, ATOLX, TOLF, IFLAG, RW, LRW, + IW, LIW) C***BEGIN PROLOGUE DSOS C***PURPOSE Solve a square system of nonlinear equations. C***LIBRARY SLATEC C***CATEGORY F2A C***TYPE DOUBLE PRECISION (SOS-S, DSOS-D) C***KEYWORDS BROWN'S METHOD, NEWTON'S METHOD, NONLINEAR EQUATIONS, C ROOTS, SOLUTIONS C***AUTHOR Watts, H. A., (SNLA) C***DESCRIPTION C C DSOS solves a system of NEQ simultaneous nonlinear equations in C NEQ unknowns. That is, it solves the problem F(X)=0 C where X is a vector with components X(1),...,X(NEQ) and F C is a vector of nonlinear functions. Each equation is of the form C C F (X(1),...,X(NEQ))=0 for K=1,...,NEQ. C K C C The algorithm is based on an iterative method which is a C variation of Newton's method using Gaussian elimination C in a manner similar to the Gauss-Seidel process. Convergence C is roughly quadratic. All partial derivatives required by C the algorithm are approximated by first difference quotients. C The convergence behavior of this code is affected by the C ordering of the equations, and it is advantageous to place linear C and mildly nonlinear equations first in the ordering. C C Actually, DSOS is merely an interfacing routine for C calling subroutine DSOSEQ which embodies the solution C algorithm. The purpose of this is to add greater C flexibility and ease of use for the prospective user. C C DSOSEQ calls the accompanying routine DSOSSL which solves special C triangular linear systems by back-substitution. C C The user must supply a function subprogram which evaluates the C K-th equation only (K specified by DSOSEQ) for each call C to the subprogram. C C DSOS represents an implementation of the mathematical algorithm C described in the references below. It is a modification of the C code SOSNLE written by H. A. Watts in 1973. C C ********************************************************************** C -Input- C C FNC -Name of the function program which evaluates the equations. C This name must be in an EXTERNAL statement in the calling C program. The user must supply FNC in the form FNC(X,K), C where X is the solution vector (which must be dimensioned C in FNC) and FNC returns the value of the K-th function. C C NEQ -Number of equations to be solved. C C X -Solution vector. Initial guesses must be supplied. C C RTOLX -Relative error tolerance used in the convergence criteria. C Each solution component X(I) is checked by an accuracy test C of the form ABS(X(I)-XOLD(I)) .LE. RTOLX*ABS(X(I))+ATOLX, C where XOLD(I) represents the previous iteration value. C RTOLX must be non-negative. C C ATOLX -Absolute error tolerance used in the convergence criteria. C ATOLX must be non-negative. If the user suspects some C solution component may be zero, he should set ATOLX to an C appropriate (depends on the scale of the remaining variables) C positive value for better efficiency. C C TOLF -Residual error tolerance used in the convergence criteria. C Convergence will be indicated if all residuals (values of the C functions or equations) are not bigger than TOLF in C magnitude. Note that extreme care must be given in assigning C an appropriate value for TOLF because this convergence test C is dependent on the scaling of the equations. An C inappropriate value can cause premature termination of the C iteration process. C C IFLAG -Optional input indicator. You must set IFLAG=-1 if you C want to use any of the optional input items listed below. C Otherwise set it to zero. C C RW -A DOUBLE PRECISION work array which is split apart by DSOS C and used internally by DSOSEQ. C C LRW -Dimension of the RW array. LRW must be at least C 1 + 6*NEQ + NEQ*(NEQ+1)/2 C C IW -An INTEGER work array which is split apart by DSOS and used C internally by DSOSEQ. C C LIW -Dimension of the IW array. LIW must be at least 3 + NEQ. C C -Optional Input- C C IW(1) -Internal printing parameter. You must set IW(1)=-1 if C you want the intermediate solution iterates to be printed. C C IW(2) -Iteration limit. The maximum number of allowable C iterations can be specified, if desired. To override the C default value of 50, set IW(2) to the number wanted. C C Remember, if you tell the code that you are using one of the C options (by setting IFLAG=-1), you must supply values C for both IW(1) and IW(2). C C ********************************************************************** C -Output- C C X -Solution vector. C C IFLAG -Status indicator C C *** Convergence to a Solution *** C C 1 Means satisfactory convergence to a solution was achieved. C Each solution component X(I) satisfies the error tolerance C test ABS(X(I)-XOLD(I)) .LE. RTOLX*ABS(X(I))+ATOLX. C C 2 Means procedure converged to a solution such that all C residuals are at most TOLF in magnitude, C ABS(FNC(X,I)) .LE. TOLF. C C 3 Means that conditions for both IFLAG=1 and IFLAG=2 hold. C C 4 Means possible numerical convergence. Behavior indicates C limiting precision calculations as a result of user asking C for too much accuracy or else convergence is very slow. C Residual norms and solution increment norms have C remained roughly constant over several consecutive C iterations. C C *** Task Interrupted *** C C 5 Means the allowable number of iterations has been met C without obtaining a solution to the specified accuracy. C Very slow convergence may be indicated. Examine the C approximate solution returned and see if the error C tolerances seem appropriate. C C 6 Means the allowable number of iterations has been met and C the iterative process does not appear to be converging. C A local minimum may have been encountered or there may be C limiting precision difficulties. C C 7 Means that the iterative scheme appears to be diverging. C Residual norms and solution increment norms have C increased over several consecutive iterations. C C *** Task Cannot Be Continued *** C C 8 Means that a Jacobian-related matrix was singular. C C 9 Means improper input parameters. C C *** IFLAG should be examined after each call to *** C *** DSOS with the appropriate action being taken. *** C C C RW(1) -Contains a norm of the residual. C C IW(3) -Contains the number of iterations used by the process. C C ********************************************************************** C C***REFERENCES K. M. Brown, Solution of simultaneous nonlinear C equations, Algorithm 316, Communications of the C A.C.M. 10, (1967), pp. 728-729. C K. M. Brown, A quadratically convergent Newton-like C method based upon Gaussian elimination, SIAM Journal C on Numerical Analysis 6, (1969), pp. 560-569. C***ROUTINES CALLED DSOSEQ, XERMSG C***REVISION HISTORY (YYMMDD) C 801001 DATE WRITTEN C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900510 Convert XERRWV calls to XERMSG calls, change Prologue C comments to agree with SOS. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DSOS INTEGER IFLAG, INPFLG, IPRINT, IW(*), K1, K2, K3, K4, K5, K6, 1 LIW, LRW, MXIT, NC, NCJS, NEQ, NSRI, NSRRC DOUBLE PRECISION ATOLX, FNC, RTOLX, RW(*), TOLF, X(*) CHARACTER*8 XERN1 CHARACTER*16 XERN3, XERN4 EXTERNAL FNC C***FIRST EXECUTABLE STATEMENT DSOS INPFLG = IFLAG C C CHECK FOR VALID INPUT C IF (NEQ .LE. 0) THEN WRITE (XERN1, '(I8)') NEQ CALL XERMSG ('SLATEC', 'DSOS', 'THE NUMBER OF EQUATIONS ' // * 'MUST BE A POSITIVE INTEGER. YOU HAVE CALLED THE ' // * 'CODE WITH NEQ = ' // XERN1, 1, 1) IFLAG = 9 ENDIF C IF (RTOLX .LT. 0.0D0 .OR. ATOLX .LT. 0.0D0) THEN WRITE (XERN3, '(1PE15.6)') ATOLX WRITE (XERN4, '(1PE15.6)') RTOLX CALL XERMSG ('SLATEC', 'DSOS', 'THE ERROR TOLERANCES FOR ' // * 'THE SOLUTION ITERATES CANNOT BE NEGATIVE. YOU HAVE ' // * 'CALLED THE CODE WITH RTOLX = ' // XERN3 // * ' AND ATOLX = ' // XERN4,2, 1) IFLAG = 9 ENDIF C IF (TOLF .LT. 0.0D0) THEN WRITE (XERN3, '(1PE15.6)') TOLF CALL XERMSG ('SLATEC', 'DSOS', 'THE RESIDUAL ERROR ' // * 'TOLERANCE MUST BE NON-NEGATIVE. YOU HAVE CALLED THE ' // * 'CODE WITH TOLF = ' // XERN3, 3, 1) IFLAG = 9 ENDIF C IPRINT = 0 MXIT = 50 IF (INPFLG .EQ. (-1)) THEN IF (IW(1) .EQ. (-1)) IPRINT = -1 MXIT = IW(2) IF (MXIT .LE. 0) THEN WRITE (XERN1, '(I8)') MXIT CALL XERMSG ('SLATEC', 'DSOS', 'YOU HAVE TOLD THE CODE ' // * 'TO USE OPTIONAL INPUT ITEMS BY SETTING IFLAG=-1. ' // * 'HOWEVER YOU HAVE CALLED THE CODE WITH THE MAXIMUM ' // * 'ALLOWABLE NUMBER OF ITERATIONS SET TO IW(2) = ' // * XERN1, 4, 1) IFLAG = 9 ENDIF ENDIF C NC = (NEQ*(NEQ+1))/2 IF (LRW .LT. 1 + 6*NEQ + NC) THEN WRITE (XERN1, '(I8)') LRW CALL XERMSG ('SLATEC', 'DSOS', 'DIMENSION OF THE RW ARRAY ' // * 'MUST BE AT LEAST 1 + 6*NEQ + NEQ*(NEQ+1)/2 . YOU HAVE ' // * 'CALLED THE CODE WITH LRW = ' // XERN1, 5, 1) IFLAG = 9 ENDIF C IF (LIW .LT. 3 + NEQ) THEN WRITE (XERN1, '(I8)') LIW CALL XERMSG ('SLATEC', 'DSOS', 'DIMENSION OF THE IW ARRAY ' // * 'MUST BE AT LEAST 3 + NEQ. YOU HAVE CALLED THE CODE ' // * 'WITH LIW = ' // XERN1, 6, 1) IFLAG = 9 ENDIF C IF (IFLAG .NE. 9) THEN NCJS = 6 NSRRC = 4 NSRI = 5 C K1 = NC + 2 K2 = K1 + NEQ K3 = K2 + NEQ K4 = K3 + NEQ K5 = K4 + NEQ K6 = K5 + NEQ C CALL DSOSEQ(FNC, NEQ, X, RTOLX, ATOLX, TOLF, IFLAG, MXIT, NCJS, 1 NSRRC, NSRI, IPRINT, RW(1), RW(2), NC, RW(K1), 2 RW(K2), RW(K3), RW(K4), RW(K5), RW(K6), IW(4)) C IW(3) = MXIT ENDIF RETURN END