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