relibc_openlibm/slatec/dsoseq.f

*DECK DSOSEQ
      SUBROUTINE DSOSEQ (FNC, N, S, RTOLX, ATOLX, TOLF, IFLAG, MXIT,
     +   NCJS, NSRRC, NSRI, IPRINT, FMAX, C, NC, B, P, TEMP, X, Y, FAC,
     +   IS)
C***BEGIN PROLOGUE  DSOSEQ
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DSOS
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (SOSEQS-S, DSOSEQ-D)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     DSOSEQ solves a system of N simultaneous nonlinear equations.
C     See the comments in the interfacing routine DSOS for a more
C     detailed description of some of the items in the calling list.
C
C **********************************************************************
C   -Input-
C
C     FNC- Function subprogram which evaluates the equations
C     N  -number of equations
C     S  -Solution vector of initial guesses
C     RTOLX-Relative error tolerance on solution components
C     ATOLX-Absolute error tolerance on solution components
C     TOLF-Residual error tolerance
C     MXIT-Maximum number of allowable iterations.
C     NCJS-Maximum number of consecutive iterative steps to perform
C          using the same triangular Jacobian matrix approximation.
C     NSRRC-Number of consecutive iterative steps for which the
C          limiting precision accuracy test must be satisfied
C          before the routine exits with IFLAG=4.
C     NSRI-Number of consecutive iterative steps for which the
C          diverging condition test must be satisfied before
C          the routine exits with IFLAG=7.
C     IPRINT-Internal printing parameter. You must set IPRINT=-1 if you
C          want the intermediate solution iterates and a residual norm
C          to be printed.
C     C   -Internal work array, dimensioned at least N*(N+1)/2.
C     NC  -Dimension of C array. NC  .GE.  N*(N+1)/2.
C     B   -Internal work array, dimensioned N.
C     P   -Internal work array, dimensioned N.
C     TEMP-Internal work array, dimensioned N.
C     X   -Internal work array, dimensioned N.
C     Y   -Internal work array, dimensioned N.
C     FAC -Internal work array, dimensioned N.
C     IS  -Internal work array, dimensioned N.
C
C   -Output-
C     S    -Solution vector
C     IFLAG-Status indicator flag
C     MXIT-The actual number of iterations performed
C     FMAX-Residual norm
C     C   -Upper unit triangular matrix which approximates the
C          forward triangularization of the full Jacobian matrix.
C          Stored in a vector with dimension at least N*(N+1)/2.
C     B   -Contains the residuals (function values) divided
C          by the corresponding components of the P vector
C     P   -Array used to store the partial derivatives. After
C          each iteration P(K) contains the maximal derivative
C          occurring in the K-th reduced equation.
C     TEMP-Array used to store the previous solution iterate.
C     X   -Solution vector. Contains the values achieved on the
C          last iteration loop upon exit from DSOS.
C     Y   -Array containing the solution increments.
C     FAC -Array containing factors used in computing numerical
C          derivatives.
C     IS  -Records the pivotal information (column interchanges)
C
C **********************************************************************
C *** Three machine dependent parameters appear in this subroutine.
C
C *** The smallest positive magnitude, zero, is defined by the function
C *** routine D1MACH(1).
C
C *** URO, the computer unit roundoff value, is defined by D1MACH(3) for
C *** machines that round or D1MACH(4) for machines that truncate.
C *** URO is the smallest positive number such that 1.+URO  .GT.  1.
C
C *** The output tape unit number, LOUN, is defined by the function
C *** I1MACH(2).
C **********************************************************************
C
C***SEE ALSO  DSOS
C***ROUTINES CALLED  D1MACH, DSOSSL, I1MACH
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C***END PROLOGUE  DSOSEQ
C
C
      INTEGER I1MACH
      DOUBLE PRECISION D1MACH
      INTEGER IC, ICR, IFLAG, IPRINT, IS(*), ISJ, ISV, IT, ITEM, ITRY,
     1     J, JK, JS, K, KD, KJ, KK, KM1, KN, KSV, L, LOUN, LS, M, MIT,
     2     MM, MXIT, N, NC, NCJS, NP1, NSRI, NSRRC
      DOUBLE PRECISION ATOLX, B(*), C(*), CSV, F, FAC(*), FACT, FDIF,
     1     FMAX, FMIN, FMXS, FN1, FN2, FNC, FP, H, HX, P(*), PMAX, RE,
     2     RTOLX, S(*), SRURO, TEMP(*), TEST, TOLF, URO, X(*), XNORM,
     3     Y(*), YJ, YN1, YN2, YN3, YNORM, YNS, ZERO
C
C     BEGIN BLOCK PERMITTING ...EXITS TO 430
C        BEGIN BLOCK PERMITTING ...EXITS TO 410
C           BEGIN BLOCK PERMITTING ...EXITS TO 390
C***FIRST EXECUTABLE STATEMENT  DSOSEQ
               URO = D1MACH(4)
               LOUN = I1MACH(2)
               ZERO = D1MACH(1)
               RE = MAX(RTOLX,URO)
               SRURO = SQRT(URO)
C
               IFLAG = 0
               NP1 = N + 1
               ICR = 0
               IC = 0
               ITRY = NCJS
               YN1 = 0.0D0
               YN2 = 0.0D0
               YN3 = 0.0D0
               YNS = 0.0D0
               MIT = 0
               FN1 = 0.0D0
               FN2 = 0.0D0
               FMXS = 0.0D0
C
C              INITIALIZE THE INTERCHANGE (PIVOTING) VECTOR AND
C              SAVE THE CURRENT SOLUTION APPROXIMATION FOR FUTURE USE.
C
               DO 10 K = 1, N
                  IS(K) = K
                  X(K) = S(K)
                  TEMP(K) = X(K)
   10          CONTINUE
C
C
C              *********************************************************
C              **** BEGIN PRINCIPAL ITERATION LOOP  ****
C              *********************************************************
C
               DO 380 M = 1, MXIT
C                 BEGIN BLOCK PERMITTING ...EXITS TO 350
C                    BEGIN BLOCK PERMITTING ...EXITS TO 240
C
                        DO 20 K = 1, N
                           FAC(K) = SRURO
   20                   CONTINUE
C
   30                   CONTINUE
C                          BEGIN BLOCK PERMITTING ...EXITS TO 180
                              KN = 1
                              FMAX = 0.0D0
C
C
C                             ******** BEGIN SUBITERATION LOOP DEFINING
C                             THE LINEARIZATION OF EACH ********
C                             EQUATION WHICH RESULTS IN THE CONSTRUCTION
C                             OF AN UPPER ******** TRIANGULAR MATRIX
C                             APPROXIMATING THE FORWARD ********
C                             TRIANGULARIZATION OF THE FULL JACOBIAN
C                             MATRIX
C
                              DO 170 K = 1, N
C                                BEGIN BLOCK PERMITTING ...EXITS TO 160
                                    KM1 = K - 1
C
C                                   BACK-SOLVE A TRIANGULAR LINEAR
C                                   SYSTEM OBTAINING IMPROVED SOLUTION
C                                   VALUES FOR K-1 OF THE VARIABLES FROM
C                                   THE FIRST K-1 EQUATIONS. THESE
C                                   VARIABLES ARE THEN ELIMINATED FROM
C                                   THE K-TH EQUATION.
C
                                    IF (KM1 .EQ. 0) GO TO 50
                                       CALL DSOSSL(K,N,KM1,Y,C,B,KN)
                                       DO 40 J = 1, KM1
                                          JS = IS(J)
                                          X(JS) = TEMP(JS) + Y(J)
   40                                  CONTINUE
   50                               CONTINUE
C
C
C                                   EVALUATE THE K-TH EQUATION AND THE
C                                   INTERMEDIATE COMPUTATION FOR THE MAX
C                                   NORM OF THE RESIDUAL VECTOR.
C
                                    F = FNC(X,K)
                                    FMAX = MAX(FMAX,ABS(F))
C
C                                   IF WE WISH TO PERFORM SEVERAL
C                                   ITERATIONS USING A FIXED
C                                   FACTORIZATION OF AN APPROXIMATE
C                                   JACOBIAN,WE NEED ONLY UPDATE THE
C                                   CONSTANT VECTOR.
C
C                                ...EXIT
                                    IF (ITRY .LT. NCJS) GO TO 160
C
C
                                    IT = 0
C
C                                   COMPUTE PARTIAL DERIVATIVES THAT ARE
C                                   REQUIRED IN THE LINEARIZATION OF THE
C                                   K-TH REDUCED EQUATION
C
                                    DO 90 J = K, N
                                       ITEM = IS(J)
                                       HX = X(ITEM)
                                       H = FAC(ITEM)*HX
                                       IF (ABS(H) .LE. ZERO)
     1                                    H = FAC(ITEM)
                                       X(ITEM) = HX + H
                                       IF (KM1 .EQ. 0) GO TO 70
                                          Y(J) = H
                                          CALL DSOSSL(K,N,J,Y,C,B,KN)
                                          DO 60 L = 1, KM1
                                             LS = IS(L)
                                             X(LS) = TEMP(LS) + Y(L)
   60                                     CONTINUE
   70                                  CONTINUE
                                       FP = FNC(X,K)
                                       X(ITEM) = HX
                                       FDIF = FP - F
                                       IF (ABS(FDIF) .GT. URO*ABS(F))
     1                                    GO TO 80
                                          FDIF = 0.0D0
                                          IT = IT + 1
   80                                  CONTINUE
                                       P(J) = FDIF/H
   90                               CONTINUE
C
                                    IF (IT .LE. (N - K)) GO TO 110
C
C                                      ALL COMPUTED PARTIAL DERIVATIVES
C                                      OF THE K-TH EQUATION ARE
C                                      EFFECTIVELY ZERO.TRY LARGER
C                                      PERTURBATIONS OF THE INDEPENDENT
C                                      VARIABLES.
C
                                       DO 100 J = K, N
                                          ISJ = IS(J)
                                          FACT = 100.0D0*FAC(ISJ)
C           ..............................EXIT
                                          IF (FACT .GT. 1.0D10)
     1                                       GO TO 390
                                          FAC(ISJ) = FACT
  100                                  CONTINUE
C                          ............EXIT
                                       GO TO 180
  110                               CONTINUE
C
C                                ...EXIT
                                    IF (K .EQ. N) GO TO 160
C
C                                   ACHIEVE A PIVOTING EFFECT BY
C                                   CHOOSING THE MAXIMAL DERIVATIVE
C                                   ELEMENT
C
                                    PMAX = 0.0D0
                                    DO 130 J = K, N
                                       TEST = ABS(P(J))
                                       IF (TEST .LE. PMAX) GO TO 120
                                          PMAX = TEST
                                          ISV = J
  120                                  CONTINUE
  130                               CONTINUE
C           ........................EXIT
                                    IF (PMAX .EQ. 0.0D0) GO TO 390
C
C                                   SET UP THE COEFFICIENTS FOR THE K-TH
C                                   ROW OF THE TRIANGULAR LINEAR SYSTEM
C                                   AND SAVE THE PARTIAL DERIVATIVE OF
C                                   LARGEST MAGNITUDE
C
                                    PMAX = P(ISV)
                                    KK = KN
                                    DO 140 J = K, N
                                       IF (J .NE. ISV)
     1                                    C(KK) = -P(J)/PMAX
                                       KK = KK + 1
  140                               CONTINUE
                                    P(K) = PMAX
C
C
C                                ...EXIT
                                    IF (ISV .EQ. K) GO TO 160
C
C                                   INTERCHANGE THE TWO COLUMNS OF C
C                                   DETERMINED BY THE PIVOTAL STRATEGY
C
                                    KSV = IS(K)
                                    IS(K) = IS(ISV)
                                    IS(ISV) = KSV
C
                                    KD = ISV - K
                                    KJ = K
                                    DO 150 J = 1, K
                                       CSV = C(KJ)
                                       JK = KJ + KD
                                       C(KJ) = C(JK)
                                       C(JK) = CSV
                                       KJ = KJ + N - J
  150                               CONTINUE
  160                            CONTINUE
C
                                 KN = KN + NP1 - K
C
C                                STORE THE COMPONENTS FOR THE CONSTANT
C                                VECTOR
C
                                 B(K) = -F/P(K)
C
  170                         CONTINUE
C                       ......EXIT
                              GO TO 190
  180                      CONTINUE
                        GO TO 30
  190                   CONTINUE
C
C                       ********
C                       ******** END OF LOOP CREATING THE TRIANGULAR
C                       LINEARIZATION MATRIX
C                       ********
C
C
C                        SOLVE THE RESULTING TRIANGULAR SYSTEM FOR A NEW
C                        SOLUTION APPROXIMATION AND OBTAIN THE SOLUTION
C                        INCREMENT NORM.
C
                        KN = KN - 1
                        Y(N) = B(N)
                        IF (N .GT. 1) CALL DSOSSL(N,N,N,Y,C,B,KN)
                        XNORM = 0.0D0
                        YNORM = 0.0D0
                        DO 200 J = 1, N
                           YJ = Y(J)
                           YNORM = MAX(YNORM,ABS(YJ))
                           JS = IS(J)
                           X(JS) = TEMP(JS) + YJ
                           XNORM = MAX(XNORM,ABS(X(JS)))
  200                   CONTINUE
C
C
C                       PRINT INTERMEDIATE SOLUTION ITERATES AND
C                       RESIDUAL NORM IF DESIRED
C
                        IF (IPRINT .NE. (-1)) GO TO 220
                           MM = M - 1
                           WRITE (LOUN,210) FMAX,MM,(X(J), J = 1, N)
  210                      FORMAT ('0RESIDUAL NORM =', D9.2, / 1X,
     1                             'SOLUTION ITERATE (', I3, ')', /
     2                             (1X, 5D26.14))
  220                   CONTINUE
C
C                       TEST FOR CONVERGENCE TO A SOLUTION (RELATIVE
C                       AND/OR ABSOLUTE ERROR COMPARISON ON SUCCESSIVE
C                       APPROXIMATIONS OF EACH SOLUTION VARIABLE)
C
                        DO 230 J = 1, N
                           JS = IS(J)
C                    ......EXIT
                           IF (ABS(Y(J)) .GT. RE*ABS(X(JS)) + ATOLX)
     1                        GO TO 240
  230                   CONTINUE
                        IF (FMAX .LE. FMXS) IFLAG = 1
  240                CONTINUE
C
C                    TEST FOR CONVERGENCE TO A SOLUTION BASED ON
C                    RESIDUALS
C
                     IF (FMAX .LE. TOLF) IFLAG = IFLAG + 2
C        ............EXIT
                     IF (IFLAG .GT. 0) GO TO 410
C
C
                     IF (M .GT. 1) GO TO 250
                        FMIN = FMAX
                     GO TO 330
  250                CONTINUE
C                       BEGIN BLOCK PERMITTING ...EXITS TO 320
C
C                          SAVE SOLUTION HAVING MINIMUM RESIDUAL NORM.
C
                           IF (FMAX .GE. FMIN) GO TO 270
                              MIT = M + 1
                              YN1 = YNORM
                              YN2 = YNS
                              FN1 = FMXS
                              FMIN = FMAX
                              DO 260 J = 1, N
                                 S(J) = X(J)
  260                         CONTINUE
                              IC = 0
  270                      CONTINUE
C
C                          TEST FOR LIMITING PRECISION CONVERGENCE. VERY
C                          SLOWLY CONVERGENT PROBLEMS MAY ALSO BE
C                          DETECTED.
C
                           IF (YNORM .GT. SRURO*XNORM) GO TO 290
                           IF (FMAX .LT. 0.2D0*FMXS
     1                         .OR. FMAX .GT. 5.0D0*FMXS) GO TO 290
                           IF (YNORM .LT. 0.2D0*YNS
     1                         .OR. YNORM .GT. 5.0D0*YNS) GO TO 290
                              ICR = ICR + 1
                              IF (ICR .GE. NSRRC) GO TO 280
                                 IC = 0
C                       .........EXIT
                                 GO TO 320
  280                         CONTINUE
                              IFLAG = 4
                              FMAX = FMIN
C     ........................EXIT
                              GO TO 430
  290                      CONTINUE
                           ICR = 0
C
C                          TEST FOR DIVERGENCE OF THE ITERATIVE SCHEME.
C
                           IF (YNORM .GT. 2.0D0*YNS
     1                         .OR. FMAX .GT. 2.0D0*FMXS) GO TO 300
                              IC = 0
                           GO TO 310
  300                      CONTINUE
                              IC = IC + 1
C                       ......EXIT
                              IF (IC .LT. NSRI) GO TO 320
                              IFLAG = 7
C        .....................EXIT
                              GO TO 410
  310                      CONTINUE
  320                   CONTINUE
  330                CONTINUE
C
C                    CHECK TO SEE IF NEXT ITERATION CAN USE THE OLD
C                    JACOBIAN FACTORIZATION
C
                     ITRY = ITRY - 1
                     IF (ITRY .EQ. 0) GO TO 340
                     IF (20.0D0*YNORM .GT. XNORM) GO TO 340
                     IF (YNORM .GT. 2.0D0*YNS) GO TO 340
C                 ......EXIT
                        IF (FMAX .LT. 2.0D0*FMXS) GO TO 350
  340                CONTINUE
                     ITRY = NCJS
  350             CONTINUE
C
C                 SAVE THE CURRENT SOLUTION APPROXIMATION AND THE
C                 RESIDUAL AND SOLUTION INCREMENT NORMS FOR USE IN THE
C                 NEXT ITERATION.
C
                  DO 360 J = 1, N
                     TEMP(J) = X(J)
  360             CONTINUE
                  IF (M .NE. MIT) GO TO 370
                     FN2 = FMAX
                     YN3 = YNORM
  370             CONTINUE
                  FMXS = FMAX
                  YNS = YNORM
C
C
  380          CONTINUE
C
C              *********************************************************
C              **** END OF PRINCIPAL ITERATION LOOP ****
C              *********************************************************
C
C
C               TOO MANY ITERATIONS, CONVERGENCE WAS NOT ACHIEVED.
               M = MXIT
               IFLAG = 5
               IF (YN1 .GT. 10.0D0*YN2 .OR. YN3 .GT. 10.0D0*YN1)
     1            IFLAG = 6
               IF (FN1 .GT. 5.0D0*FMIN .OR. FN2 .GT. 5.0D0*FMIN)
     1            IFLAG = 6
               IF (FMAX .GT. 5.0D0*FMIN) IFLAG = 6
C        ......EXIT
               GO TO 410
  390       CONTINUE
C
C
C           A JACOBIAN-RELATED MATRIX IS EFFECTIVELY SINGULAR.
            IFLAG = 8
            DO 400 J = 1, N
               S(J) = TEMP(J)
  400       CONTINUE
C     ......EXIT
            GO TO 430
  410    CONTINUE
C
C
         DO 420 J = 1, N
            S(J) = X(J)
  420    CONTINUE
  430 CONTINUE
C
C
      MXIT = M
      RETURN
      END