# snsq.f

          SUBROUTINE SNSQ (FCN, JAC, IOPT, N, X, FVEC, FJAC, LDFJAC, XTOL,
         +   MAXFEV, ML, MU, EPSFCN, DIAG, MODE, FACTOR, NPRINT, INFO, NFEV,
         +   NJEV, R, LR, QTF, WA1, WA2, WA3, WA4)
    C***BEGIN PROLOGUE  SNSQ
    C***PURPOSE  Find a zero of a system of a N nonlinear functions in N
    C            variables by a modification of the Powell hybrid method.
    C***LIBRARY   SLATEC
    C***CATEGORY  F2A
    C***TYPE      SINGLE PRECISION (SNSQ-S, DNSQ-D)
    C***KEYWORDS  NONLINEAR SQUARE SYSTEM, POWELL HYBRID METHOD, ZEROS
    C***AUTHOR  Hiebert, K. L., (SNLA)
    C***DESCRIPTION
    C
    C 1. Purpose.
    C
    C       The purpose of SNSQ 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.  The user must provide a subroutine which calcu-
    C       lates 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.
    C       This code is the combination of the MINPACK codes (Argonne)
    C       HYBRD and HYBRDJ.
    C
    C
    C 2. Subroutine and Type Statements.
    C
    C       SUBROUTINE SNSQ(FCN,JAC,IOPT,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,
    C      *                 ML,MU,EPSFCN,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,
    C      *                 NJEV,R,LR,QTF,WA1,WA2,WA3,WA4)
    C       INTEGER IOPT,N,MAXFEV,ML,MU,MODE,NPRINT,INFO,NFEV,LDFJAC,NJEV,LR
    C       REAL XTOL,EPSFCN,FACTOR
    C       REAL X(N),FVEC(N),DIAG(N),FJAC(LDFJAC,N),R(LR),QTF(N),
    C      *     WA1(N),WA2(N),WA3(N),WA4(N)
    C       EXTERNAL FCN,JAC
    C
    C
    C 3. Parameters.
    C
    C       Parameters designated as input parameters must be specified on
    C       entry to SNSQ 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 SNSQ.
    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 SNSQ.  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 SNSQ.  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       FJAC is an output N by N array which contains the orthogonal
    C         matrix Q produced by the QR factorization of the final approx-
    C         imate Jacobian.
    C
    C       LDFJAC is a positive integer input variable not less than N
    C         which specifies the leading dimension of the array FJAC.
    C
    C       XTOL is a non-negative input variable.  Termination occurs when
    C         the relative error between two consecutive iterates is at most
    C         XTOL.  Therefore, XTOL measures the relative error desired in
    C         the approximate solution.  Section 4 contains more details
    C         about XTOL.
    C
    C       MAXFEV is a positive integer input variable.  Termination occurs
    C         when the number of calls to FCN is at least MAXFEV by the end
    C         of an iteration.
    C
    C       ML is a non-negative integer input variable which specifies the
    C         number of subdiagonals within the band of the Jacobian matrix.
    C         If the Jacobian is not banded or IOPT=1, set ML to at
    C         least N - 1.
    C
    C       MU is a non-negative integer input variable which specifies the
    C         number of superdiagonals within the band of the Jacobian
    C         matrix.  If the Jacobian is not banded or IOPT=1, set MU to at
    C         least N - 1.
    C
    C       EPSFCN is an input variable used in determining a suitable step
    C         for the forward-difference approximation.  This approximation
    C         assumes that the relative errors in the functions are of the
    C         order of EPSFCN.  If EPSFCN is less than the machine preci-
    C         sion, it is assumed that the relative errors in the functions
    C         are of the order of the machine precision.  If IOPT=1, then
    C         EPSFCN can be ignored (treat it as a dummy argument).
    C
    C       DIAG is an array of length N.  If MODE = 1 (see below), DIAG is
    C         internally set.  If MODE = 2, DIAG must contain positive
    C         entries that serve as implicit (multiplicative) scale factors
    C         for the variables.
    C
    C       MODE is an integer input variable.  If MODE = 1, the variables
    C         will be scaled internally.  If MODE = 2, the scaling is speci-
    C         fied by the input DIAG.  Other values of MODE are equivalent
    C         to MODE = 1.
    C
    C       FACTOR is a positive input variable used in determining the ini-
    C         tial step bound.  This bound is set to the product of FACTOR
    C         and the Euclidean norm of DIAG*X if nonzero, or else to FACTOR
    C         itself.  In most cases FACTOR should lie in the interval
    C         (.1,100.).  100. is a generally recommended value.
    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  relative error between two consecutive iterates is
    C                   at most XTOL.
    C
    C         INFO = 2  number of calls to FCN has reached or exceeded
    C                   MAXFEV.
    C
    C         INFO = 3  XTOL is too small.  No further improvement in the
    C                   approximate solution X is possible.
    C
    C         INFO = 4  iteration is not making good progress, as measured
    C                   by the improvement from the last five Jacobian eval-
    C                   uations.
    C
    C         INFO = 5  iteration is not making good progress, as measured
    C                   by the improvement from the last ten iterations.
    C
    C         Sections 4 and 5 contain more details about INFO.
    C
    C       NFEV is an integer output variable set to the number of calls to
    C         FCN.
    C
    C       NJEV is an integer output variable set to the number of calls to
    C         JAC. (If IOPT=2, then NJEV is set to zero.)
    C
    C       R is an output array of length LR which contains the upper
    C         triangular matrix produced by the QR factorization of the
    C         final approximate Jacobian, stored rowwise.
    C
    C       LR is a positive integer input variable not less than
    C         (N*(N+1))/2.
    C
    C       QTF is an output array of length N which contains the vector
    C         (Q TRANSPOSE)*FVEC.
    C
    C       WA1, WA2, WA3, and WA4 are work arrays of length N.
    C
    C
    C 4. Successful Completion.
    C
    C       The accuracy of SNSQ is controlled by the convergence parameter
    C       XTOL.  This parameter is used in a test which makes a comparison
    C       between the approximation X and a solution XSOL.  SNSQ termi-
    C       nates when the test is satisfied.  If the convergence parameter
    C       is less than the machine precision (as defined by the function
    C       R1MACH(4)), then SNSQ only attempts to satisfy the test
    C       defined by the machine precision.  Further progress is not
    C       usually possible.
    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 are coded consistently.  If these conditions
    C       are not satisfied, then SNSQ may incorrectly indicate conver-
    C       gence.  The coding of the Jacobian can be checked by the
    C       subroutine CHKDER. If the Jacobian is coded correctly or IOPT=2,
    C       then the validity of the answer can be checked, for example, by
    C       rerunning SNSQ with a tighter tolerance.
    C
    C       Convergence Test.  If ENORM(Z) denotes the Euclidean norm of a
    C         vector Z and D is the diagonal matrix whose entries are
    C         defined by the array DIAG, then this test attempts to guaran-
    C         tee that
    C
    C               ENORM(D*(X-XSOL)) .LE. XTOL*ENORM(D*XSOL).
    C
    C         If this condition is satisfied with XTOL = 10**(-K), then the
    C         larger components of D*X have K significant decimal digits and
    C         INFO is set to 1.  There is a danger that the smaller compo-
    C         nents of D*X may have large relative errors, but the fast rate
    C         of convergence of SNSQ usually avoids this possibility.
    C         Unless high precision solutions are required, the recommended
    C         value for XTOL is the square root of the machine precision.
    C
    C
    C 5. Unsuccessful Completion.
    C
    C       Unsuccessful termination of SNSQ can be due to improper input
    C       parameters, arithmetic interrupts, an excessive number of func-
    C       tion evaluations, or lack of good progress.
    C
    C       Improper Input Parameters.  INFO is set to 0 if IOPT .LT. 1,
    C         or IOPT .GT. 2, or N .LE. 0, or LDFJAC .LT. N, or
    C         XTOL .LT. 0.E0, or MAXFEV .LE. 0, or ML .LT. 0, or MU .LT. 0,
    C         or FACTOR .LE. 0.E0, or LR .LT. (N*(N+1))/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 SNSQ.  In this
    C         case, it may be possible to remedy the situation by rerunning
    C         SNSQ with a smaller value of FACTOR.
    C
    C       Excessive Number of Function Evaluations.  A reasonable value
    C         for MAXFEV is 100*(N+1) for IOPT=1 and 200*(N+1) for IOPT=2.
    C         If the number of calls to FCN reaches MAXFEV, then this
    C         indicates that the routine is converging very slowly as
    C         measured by the progress of FVEC, and INFO is set to 2.  This
    C         situation should be unusual because, as indicated below, lack
    C         of good progress is usually diagnosed earlier by SNSQ,
    C         causing termination with INFO = 4 or INFO = 5.
    C
    C       Lack of Good Progress.  SNSQ 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 SNSQ from a dif-
    C         ferent starting point may be helpful.
    C
    C
    C 6. Characteristics of the Algorithm.
    C
    C       SNSQ is a modification of the Powell hybrid method.  Two of its
    C       main characteristics involve the choice of the correction as a
    C       convex combination of the Newton and scaled gradient directions,
    C       and the updating of the Jacobian by the rank-1 method of Broy-
    C       den.  The choice of the correction guarantees (under reasonable
    C       conditions) global convergence for starting points far from the
    C       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 SNSQ 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 SNSQ 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 SNSQ will be strongly influenced by the time
    C         spent in FCN and JAC.
    C
    C       Storage.  SNSQ 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
    C C
    C C     Driver for SNSQ example.
    C C
    C       INTEGER J,IOPT,N,MAXFEV,ML,MU,MODE,NPRINT,INFO,NFEV,LDFJAC,LR,
    C      *        NWRITE
    C       REAL XTOL,EPSFCN,FACTOR,FNORM
    C       REAL X(9),FVEC(9),DIAG(9),FJAC(9,9),R(45),QTF(9),
    C      *     WA1(9),WA2(9),WA3(9),WA4(9)
    C       REAL ENORM,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
    C       LDFJAC = 9
    C       LR = 45
    C C
    C C     Set XTOL 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       XTOL = SQRT(R1MACH(4))
    C C
    C       MAXFEV = 2000
    C       ML = 1
    C       MU = 1
    C       EPSFCN = 0.E0
    C       MODE = 2
    C       DO 20 J = 1, 9
    C          DIAG(J) = 1.E0
    C    20    CONTINUE
    C       FACTOR = 1.E2
    C       NPRINT = 0
    C C
    C       CALL SNSQ(FCN,JAC,IOPT,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,ML,MU,
    C      *           EPSFCN,DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,
    C      *           R,LR,QTF,WA1,WA2,WA3,WA4)
    C       FNORM = ENORM(N,FVEC)
    C       WRITE (NWRITE,1000) FNORM,NFEV,INFO,(X(J),J=1,N)
    C       STOP
    C  1000 FORMAT (5X,' FINAL L2 NORM OF THE RESIDUALS',E15.7 //
    C      *        5X,' NUMBER OF FUNCTION EVALUATIONS',I10 //
    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       IF (IFLAG .NE. 0) GO TO 5
    C C
    C C     Insert print statements here when NPRINT is positive.
    C C
    C       RETURN
    C     5 CONTINUE
    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       NUMBER OF FUNCTION EVALUATIONS        14
    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
    C***REFERENCES  M. J. D. Powell, A hybrid method for nonlinear equa-
    C                 tions. In Numerical Methods for Nonlinear Algebraic
    C                 Equations, P. Rabinowitz, Editor.  Gordon and Breach,
    C                 1988.
    C***ROUTINES CALLED  DOGLEG, ENORM, FDJAC1, QFORM, QRFAC, R1MACH,
    C                    R1MPYQ, R1UPDT, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   800301  DATE WRITTEN
    C   890531  Changed all specific intrinsics to generic.  (WRB)
    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   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    C   920501  Reformatted the REFERENCES section.  (WRB)
    C***END PROLOGUE  SNSQ