dsos.f

      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