      SUBROUTINE SOSEQS (FNC, N, S, RTOLX, ATOLX, TOLF, IFLAG, MXIT,
     +   NCJS, NSRRC, NSRI, IPRINT, FMAX, C, NC, B, P, TEMP, X, Y, FAC,
     +   IS)
C
C
      DIMENSION S(*), C(NC), B(*), IS(*), P(*), TEMP(*), X(*), Y(*),
     1          FAC(*)
C
C***FIRST EXECUTABLE STATEMENT  SOSEQS
      URO = R1MACH(4)
      LOUN = I1MACH(2)
      ZERO = R1MACH(1)
      RE = MAX(RTOLX,URO)
      SRURO = SQRT(URO)
C
      IFLAG = 0
      NP1 = N + 1
      ICR = 0
      IC = 0
      ITRY = NCJS
      YN1 = 0.
      YN2 = 0.
      YN3 = 0.
      YNS = 0.
      MIT = 0
      FN1 = 0.
      FN2 = 0.
      FMXS = 0.
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 330 M=1,MXIT
C
        DO 20 K=1,N
          FAC(K) = SRURO
   20   CONTINUE
C
   30   KN = 1
        FMAX = 0.
C
C
C    ******** BEGIN SUBITERATION LOOP DEFINING THE LINEARIZATION OF EACH
C    ******** EQUATION WHICH RESULTS IN THE CONSTRUCTION OF AN UPPER
C    ******** TRIANGULAR MATRIX APPROXIMATING THE FORWARD
C    ******** TRIANGULARIZATION OF THE FULL JACOBIAN MATRIX
C
        DO 170 K=1,N
          KM1 = K - 1
C
C     BACK-SOLVE A TRIANGULAR LINEAR SYSTEM OBTAINING
C     IMPROVED SOLUTION VALUES FOR K-1 OF THE VARIABLES
C     FROM THE FIRST K-1 EQUATIONS. THESE VARIABLES ARE THEN
C     ELIMINATED FROM THE K-TH EQUATION.
C
          IF (KM1 .EQ. 0) GO TO 50
          CALL SOSSOL(K, N, KM1, Y, C, B, KN)
          DO 40 J=1,KM1
            JS = IS(J)
            X(JS) = TEMP(JS) + Y(J)
   40     CONTINUE
C
C
C     EVALUATE THE K-TH EQUATION AND THE INTERMEDIATE COMPUTATION
C     FOR THE MAX NORM OF THE RESIDUAL VECTOR.
C
   50     F = FNC(X,K)
          FMAX = MAX(FMAX,ABS(F))
C
C     IF WE WISH TO PERFORM SEVERAL ITERATIONS USING A FIXED
C     FACTORIZATION OF AN APPROXIMATE JACOBIAN,WE NEED ONLY
C     UPDATE THE CONSTANT VECTOR.
C
          IF (ITRY .LT. NCJS) GO TO 160
C
C
          IT = 0
C
C     COMPUTE PARTIAL DERIVATIVES THAT ARE REQUIRED IN THE LINEARIZATION
C     OF THE 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) H = FAC(ITEM)
            X(ITEM) = HX + H
            IF (KM1 .EQ. 0) GO TO 70
            Y(J) = H
            CALL SOSSOL(K, N, J, Y, C, B, KN)
            DO 60 L=1,KM1
              LS = IS(L)
              X(LS) = TEMP(LS) + Y(L)
   60       CONTINUE
   70       FP = FNC(X,K)
            X(ITEM) = HX
            FDIF = FP - F
            IF (ABS(FDIF) .GT. URO*ABS(F)) GO TO 80
            FDIF = 0.
            IT = IT + 1
   80       P(J) = FDIF/H
   90     CONTINUE
C
          IF (IT .LE. (N-K)) GO TO 110
C
C     ALL COMPUTED PARTIAL DERIVATIVES OF THE K-TH EQUATION
C     ARE EFFECTIVELY ZERO.TRY LARGER PERTURBATIONS OF THE
C     INDEPENDENT VARIABLES.
C
          DO 100 J=K,N
            ISJ = IS(J)
            FACT = 100.*FAC(ISJ)
            IF (FACT .GT. 1.E+10) GO TO 340
            FAC(ISJ) = FACT
  100     CONTINUE
          GO TO 30
C
  110     IF (K .EQ. N) GO TO 160
C
C     ACHIEVE A PIVOTING EFFECT BY CHOOSING THE MAXIMAL DERIVATIVE
C     ELEMENT
C
          PMAX = 0.
          DO 120 J=K,N
            TEST = ABS(P(J))
            IF (TEST .LE. PMAX) GO TO 120
            PMAX = TEST
            ISV = J
  120     CONTINUE
          IF (PMAX .EQ. 0.) GO TO 340
C
C     SET UP THE COEFFICIENTS FOR THE K-TH ROW OF THE TRIANGULAR
C     LINEAR SYSTEM AND SAVE THE PARTIAL DERIVATIVE OF
C     LARGEST MAGNITUDE
C
          PMAX = P(ISV)
          KK = KN
          DO 140 J=K,N
            IF (J .EQ. ISV) GO TO 130
            C(KK) = -P(J)/PMAX
  130       KK = KK + 1
  140     CONTINUE
          P(K) = PMAX
C
C
          IF (ISV .EQ. K) GO TO 160
C
C     INTERCHANGE THE TWO COLUMNS OF C DETERMINED BY THE
C     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
C
  160     KN = KN + NP1 - K
C
C     STORE THE COMPONENTS FOR THE CONSTANT VECTOR
C
          B(K) = -F/P(K)
C
  170   CONTINUE
C
C    ********
C    ******** END OF LOOP CREATING THE TRIANGULAR LINEARIZATION MATRIX
C    ********
C
C
C     SOLVE THE RESULTING TRIANGULAR SYSTEM FOR A NEW SOLUTION
C     APPROXIMATION AND OBTAIN THE SOLUTION INCREMENT NORM.
C
        KN = KN - 1
        Y(N) = B(N)
        IF (N .GT. 1) CALL SOSSOL(N, N, N, Y, C, B, KN)
        XNORM = 0.
        YNORM = 0.
        DO 180 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)))
  180   CONTINUE
C
C
C     PRINT INTERMEDIATE SOLUTION ITERATES AND RESIDUAL NORM IF DESIRED
C
        IF (IPRINT.NE.(-1)) GO TO 190
        MM = M - 1
        WRITE (LOUN,1234) FMAX, MM, (X(J),J=1,N)
 1234   FORMAT ('0RESIDUAL NORM =', E9.2, /1X, 'SOLUTION ITERATE',
     1   ' (', I3, ')', /(1X, 5E26.14))
  190   CONTINUE
C
C     TEST FOR CONVERGENCE TO A SOLUTION (RELATIVE AND/OR ABSOLUTE ERROR
C     COMPARISON ON SUCCESSIVE APPROXIMATIONS OF EACH SOLUTION VARIABLE)
C
        DO 200 J=1,N
          JS = IS(J)
          IF (ABS(Y(J)) .GT. RE*ABS(X(JS))+ATOLX) GO TO 210
  200   CONTINUE
        IF (FMAX .LE. FMXS) IFLAG = 1
C
C     TEST FOR CONVERGENCE TO A SOLUTION BASED ON RESIDUALS
C
  210   IF (FMAX .GT. TOLF) GO TO 220
        IFLAG = IFLAG + 2
  220   IF (IFLAG .GT. 0) GO TO 360
C
C
        IF (M .GT. 1) GO TO 230
        FMIN = FMAX
        GO TO 280
C
C     SAVE SOLUTION HAVING MINIMUM RESIDUAL NORM.
C
  230   IF (FMAX .GE. FMIN) GO TO 250
        MIT = M + 1
        YN1 = YNORM
        YN2 = YNS
        FN1 = FMXS
        FMIN = FMAX
        DO 240 J=1,N
          S(J) = X(J)
  240   CONTINUE
        IC = 0
C
C     TEST FOR LIMITING PRECISION CONVERGENCE.  VERY SLOWLY CONVERGENT
C     PROBLEMS MAY ALSO BE DETECTED.
C
  250   IF (YNORM .GT. SRURO*XNORM) GO TO 260
        IF ((FMAX .LT. 0.2*FMXS) .OR. (FMAX .GT. 5.*FMXS)) GO TO 260
        IF ((YNORM .LT. 0.2*YNS) .OR. (YNORM .GT. 5.*YNS)) GO TO 260
        ICR = ICR + 1
        IF (ICR .LT. NSRRC) GO TO 270
        IFLAG = 4
        FMAX = FMIN
        GO TO 380
  260   ICR = 0
C
C     TEST FOR DIVERGENCE OF THE ITERATIVE SCHEME.
C
        IF ((YNORM .LE. 2.*YNS) .AND. (FMAX .LE. 2.*FMXS)) GO TO 270
        IC = IC + 1
        IF (IC .LT. NSRI) GO TO 280
        IFLAG = 7
        GO TO 360
  270   IC = 0
C
C     CHECK TO SEE IF NEXT ITERATION CAN USE THE OLD JACOBIAN
C     FACTORIZATION
C
  280   ITRY = ITRY - 1
        IF (ITRY .EQ. 0) GO TO 290
        IF (20.*YNORM .GT. XNORM) GO TO 290
        IF (YNORM .GT. 2.*YNS) GO TO 290
        IF (FMAX .LT. 2.*FMXS) GO TO 300
  290   ITRY = NCJS
C
C     SAVE THE CURRENT SOLUTION APPROXIMATION AND THE RESIDUAL AND
C     SOLUTION INCREMENT NORMS FOR USE IN THE NEXT ITERATION.
C
  300   DO 310 J=1,N
          TEMP(J) = X(J)
  310   CONTINUE
        IF (M.NE.MIT) GO TO 320
        FN2 = FMAX
        YN3 = YNORM
  320   FMXS = FMAX
        YNS = YNORM
C
C
  330 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.0*YN2 .OR. YN3 .GT. 10.0*YN1) IFLAG = 6
      IF (FN1 .GT. 5.0*FMIN .OR. FN2 .GT. 5.0*FMIN) IFLAG = 6
      IF (FMAX .GT. 5.0*FMIN) IFLAG = 6
      GO TO 360
C
C
C     A JACOBIAN-RELATED MATRIX IS EFFECTIVELY SINGULAR.
  340 IFLAG = 8
      DO 350 J=1,N
        S(J) = TEMP(J)
  350 CONTINUE
      GO TO 380
C
C
  360 DO 370 J=1,N
        S(J) = X(J)
  370 CONTINUE
C
C
  380 MXIT = M
      RETURN
      END
