*DECK SOS SUBROUTINE SOS (FNC, NEQ, X, RTOLX, ATOLX, TOLF, IFLAG, RW, LRW, + IW, LIW) C***BEGIN PROLOGUE SOS C***PURPOSE Solve a square system of nonlinear equations. C***LIBRARY SLATEC C***CATEGORY F2A C***TYPE SINGLE 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 SOS 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, SOS is merely an interfacing routine for C calling subroutine SOSEQS 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 SOSEQS calls the accompanying routine SOSSOL, 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 SOSEQS) for each call C to the subprogram. C C SOS 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 REAL work array which is split apart by SOS and used C internally by SOSEQS. 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 SOS and used C internally by SOSEQS. 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 *** SOS 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***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 SOSEQS, 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, changed Prologue C comments to agree with DSOS. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE SOS