1      SUBROUTINE SDCOR (DFDY,EL,FA,H,IMPL,IPVT,MATDIM,MITER,ML,MU,N,
2     8   NDE,NQ,T,Y,YH,YWT,EVALFA,SAVE1,SAVE2,A,D)
3C***BEGIN PROLOGUE  SDCOR
4C***REFER TO  SDRIV3
5C  Subroutine SDCOR is called to compute corrections to the Y array.
6C  In the case of functional iteration, update Y directly from the
7C  result of the last call to F.
8C  In the case of the chord method, compute the corrector error and
9C  solve the linear system with that as right hand side and DFDY as
10C  coefficient matrix, using the LU decomposition if MITER is 1, 2, 4,
11C  or 5.
12C***ROUTINES CALLED  SGESL,SGBSL,SNRM2
13C***DATE WRITTEN   790601   (YYMMDD)
14C***REVISION DATE  841119   (YYMMDD)
15C***CATEGORY NO.  I1A2,I1A1B
16C***AUTHOR  KAHANER, D. K., NATIONAL BUREAU OF STANDARDS,
17C           SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY
18C***END PROLOGUE  SDCOR
19      REAL A(MATDIM,*), D, DFDY(MATDIM,*), EL(13,12), H,
20     8     SAVE1(*), SAVE2(*), SNRM2, T, Y(*), YH(N,*), YWT(*)
21      INTEGER IPVT(*)
22      LOGICAL EVALFA
23C***FIRST EXECUTABLE STATEMENT  SDCOR
24      IF (MITER .EQ. 0) THEN
25        DO 100 I = 1,N
26 100      SAVE1(I) = (H*SAVE2(I) - YH(I,2) - SAVE1(I))/YWT(I)
27        D = SNRM2(N, SAVE1, 1)/SQRT(REAL(N))
28        DO 105 I = 1,N
29 105      SAVE1(I) = H*SAVE2(I) - YH(I,2)
30      ELSE IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN
31        IF (IMPL .EQ. 0) THEN
32          DO 130 I = 1,N
33 130        SAVE2(I) = H*SAVE2(I) - YH(I,2) - SAVE1(I)
34        ELSE IF (IMPL .EQ. 1) THEN
35          IF (EVALFA) THEN
36            CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
37          ELSE
38            EVALFA = .TRUE.
39          END IF
40          DO 150 I = 1,N
41 150        SAVE2(I) = H*SAVE2(I)
42          DO 160 J = 1,N
43            DO 160 I = 1,N
44 160          SAVE2(I) = SAVE2(I) - A(I,J)*(YH(J,2) + SAVE1(J))
45        ELSE IF (IMPL .EQ. 2) THEN
46          IF (EVALFA) THEN
47            CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
48          ELSE
49            EVALFA = .TRUE.
50          END IF
51          DO 180 I = 1,N
52 180        SAVE2(I) = H*SAVE2(I) - A(I,1)*(YH(I,2) + SAVE1(I))
53        END IF
54        CALL SGESL (DFDY, MATDIM, N, IPVT, SAVE2, 0)
55        DO 200 I = 1,N
56          SAVE1(I) = SAVE1(I) + SAVE2(I)
57 200      SAVE2(I) = SAVE2(I)/YWT(I)
58        D = SNRM2(N, SAVE2, 1)/SQRT(REAL(N))
59      ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN
60        IF (IMPL .EQ. 0) THEN
61          DO 230 I = 1,N
62 230        SAVE2(I) = H*SAVE2(I) - YH(I,2) - SAVE1(I)
63        ELSE IF (IMPL .EQ. 1) THEN
64          IF (EVALFA) THEN
65            CALL FA (N, T, Y, A(ML+1,1), MATDIM, ML, MU, NDE)
66          ELSE
67            EVALFA = .TRUE.
68          END IF
69          DO 250 I = 1,N
70 250        SAVE2(I) = H*SAVE2(I)
71          MW = ML + 1 + MU
72          DO 260 J = 1,N
73            I1 = MAX(ML+1, MW+1-J)
74            I2 = MIN(MW+N-J, MW+ML)
75            DO 260 I = I1,I2
76              I3 = I + J - MW
77 260          SAVE2(I3) = SAVE2(I3) - A(I,J)*(YH(J,2) + SAVE1(J))
78        ELSE IF (IMPL .EQ. 2) THEN
79          IF (EVALFA) THEN
80            CALL FA (N, T, Y, A, MATDIM, ML, MU, NDE)
81          ELSE
82            EVALFA = .TRUE.
83          END IF
84          DO 280 I = 1,N
85 280        SAVE2(I) = H*SAVE2(I) - A(I,1)*(YH(I,2) + SAVE1(I))
86        END IF
87        CALL SGBSL (DFDY, MATDIM, N, ML, MU, IPVT, SAVE2, 0)
88        DO 300 I = 1,N
89          SAVE1(I) = SAVE1(I) + SAVE2(I)
90 300      SAVE2(I) = SAVE2(I)/YWT(I)
91        D = SNRM2(N, SAVE2, 1)/SQRT(REAL(N))
92      ELSE IF (MITER .EQ. 3) THEN
93        IFLAG = 2
94        CALL USERS (Y, YH(1,2), YWT, SAVE1, SAVE2, T, H, EL(1,NQ), IMPL,
95     8              N, NDE, IFLAG)
96        DO 320 I = 1,N
97          SAVE1(I) = SAVE1(I) + SAVE2(I)
98 320      SAVE2(I) = SAVE2(I)/YWT(I)
99        D = SNRM2(N, SAVE2, 1)/SQRT(REAL(N))
100      END IF
101      END
102