1*DECK PJAC
2      SUBROUTINE PJAC (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F,
3     +   JAC, RPAR, IPAR)
4C***BEGIN PROLOGUE  PJAC
5C***SUBSIDIARY
6C***PURPOSE  Subsidiary to DEBDF
7C***LIBRARY   SLATEC
8C***TYPE      SINGLE PRECISION (PJAC-S, DPJAC-D)
9C***AUTHOR  Watts, H. A., (SNLA)
10C***DESCRIPTION
11C
12C   PJAC sets up the iteration matrix (involving the Jacobian) for the
13C   integration package DEBDF.
14C
15C***SEE ALSO  DEBDF
16C***ROUTINES CALLED  SGBFA, SGEFA, VNWRMS
17C***COMMON BLOCKS    DEBDF1
18C***REVISION HISTORY  (YYMMDD)
19C   800901  DATE WRITTEN
20C   890531  Changed all specific intrinsics to generic.  (WRB)
21C   891214  Prologue converted to Version 4.0 format.  (BAB)
22C   900328  Added TYPE section.  (WRB)
23C   910722  Updated AUTHOR section.  (ALS)
24C   920422  Changed DIMENSION statement.  (WRB)
25C***END PROLOGUE  PJAC
26C
27CLLL. OPTIMIZE
28      INTEGER NEQ, NYH, IWM, I, I1, I2, IER, II, IOWND, IOWNS, J, J1,
29     1   JJ, JSTART, KFLAG, L, LENP, MAXORD, MBA, MBAND, MEB1, MEBAND,
30     2   METH, MITER, ML, ML3, MU, N, NFE, NJE, NQ, NQU, NST
31      EXTERNAL F, JAC
32      REAL Y, YH, EWT, FTEM, SAVF, WM,
33     1   ROWND, ROWNS, EL0, H, HMIN, HMXI, HU, TN, UROUND,
34     2   CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ, VNWRMS
35      DIMENSION         Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*),
36     1   WM(*), IWM(*), RPAR(*), IPAR(*)
37      COMMON /DEBDF1/ ROWND, ROWNS(210),
38     1   EL0, H, HMIN, HMXI, HU, TN, UROUND, IOWND(14), IOWNS(6),
39     2   IER, JSTART, KFLAG, L, METH, MITER, MAXORD, N, NQ, NST, NFE,
40     3   NJE, NQU
41C-----------------------------------------------------------------------
42C PJAC IS CALLED BY STOD  TO COMPUTE AND PROCESS THE MATRIX
43C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
44C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF
45C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
46C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
47C J IS STORED IN WM AND REPLACED BY P.  IF MITER .NE. 3, P IS THEN
48C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
49C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
50C BY SGEFA IF MITER = 1 OR 2, AND BY SGBFA IF MITER = 4 OR 5.
51C
52C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
53C WITH PJAC USES THE FOLLOWING..
54C Y    = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
55C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STOD ).
56C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
57C WM   = REAL WORK SPACE FOR MATRICES.  ON OUTPUT IT CONTAINS THE
58C        INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
59C        OF P IF MITER IS 1, 2 , 4, OR 5.
60C        STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
61C        WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
62C        WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
63C        WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
64C IWM  = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
65C        IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS THE
66C        BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
67C EL0  = EL(1) (INPUT).
68C IER  = OUTPUT ERROR FLAG,  = 0 IF NO TROUBLE, .NE. 0 IF
69C        P MATRIX FOUND TO BE SINGULAR.
70C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
71C MITER, N, NFE, AND NJE.
72C-----------------------------------------------------------------------
73C***FIRST EXECUTABLE STATEMENT  PJAC
74      NJE = NJE + 1
75      HL0 = H*EL0
76      GO TO (100, 200, 300, 400, 500), MITER
77C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
78 100  LENP = N*N
79      DO 110 I = 1,LENP
80 110    WM(I+2) = 0.0E0
81      CALL JAC (TN, Y, WM(3), N, RPAR, IPAR)
82      CON = -HL0
83      DO 120 I = 1,LENP
84 120    WM(I+2) = WM(I+2)*CON
85      GO TO 240
86C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
87 200  FAC = VNWRMS (N, SAVF, EWT)
88      R0 = 1000.0E0*ABS(H)*UROUND*N*FAC
89      IF (R0 .EQ. 0.0E0) R0 = 1.0E0
90      SRUR = WM(1)
91      J1 = 2
92      DO 230 J = 1,N
93        YJ = Y(J)
94        R = MAX(SRUR*ABS(YJ),R0*EWT(J))
95        Y(J) = Y(J) + R
96        FAC = -HL0/R
97        CALL F (TN, Y, FTEM, RPAR, IPAR)
98        DO 220 I = 1,N
99 220      WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
100        Y(J) = YJ
101        J1 = J1 + N
102 230    CONTINUE
103      NFE = NFE + N
104C ADD IDENTITY MATRIX. -------------------------------------------------
105 240  J = 3
106      DO 250 I = 1,N
107        WM(J) = WM(J) + 1.0E0
108 250    J = J + (N + 1)
109C DO LU DECOMPOSITION ON P. --------------------------------------------
110      CALL SGEFA (WM(3), N, N, IWM(21), IER)
111      RETURN
112C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
113 300  WM(2) = HL0
114      IER = 0
115      R = EL0*0.1E0
116      DO 310 I = 1,N
117 310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
118      CALL F (TN, Y, WM(3), RPAR, IPAR)
119      NFE = NFE + 1
120      DO 320 I = 1,N
121        R0 = H*SAVF(I) - YH(I,2)
122        DI = 0.1E0*R0 - H*(WM(I+2) - SAVF(I))
123        WM(I+2) = 1.0E0
124        IF (ABS(R0) .LT. UROUND*EWT(I)) GO TO 320
125        IF (ABS(DI) .EQ. 0.0E0) GO TO 330
126        WM(I+2) = 0.1E0*R0/DI
127 320    CONTINUE
128      RETURN
129 330  IER = -1
130      RETURN
131C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
132 400  ML = IWM(1)
133      MU = IWM(2)
134      ML3 =  3
135      MBAND = ML + MU + 1
136      MEBAND = MBAND + ML
137      LENP = MEBAND*N
138      DO 410 I = 1,LENP
139 410    WM(I+2) = 0.0E0
140      CALL JAC (TN, Y, WM(ML3), MEBAND, RPAR, IPAR)
141      CON = -HL0
142      DO 420 I = 1,LENP
143 420    WM(I+2) = WM(I+2)*CON
144      GO TO 570
145C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
146 500  ML = IWM(1)
147      MU = IWM(2)
148      MBAND = ML + MU + 1
149      MBA = MIN(MBAND,N)
150      MEBAND = MBAND + ML
151      MEB1 = MEBAND - 1
152      SRUR = WM(1)
153      FAC = VNWRMS (N, SAVF, EWT)
154      R0 = 1000.0E0*ABS(H)*UROUND*N*FAC
155      IF (R0 .EQ. 0.0E0) R0 = 1.0E0
156      DO 560 J = 1,MBA
157        DO 530 I = J,N,MBAND
158          YI = Y(I)
159          R = MAX(SRUR*ABS(YI),R0*EWT(I))
160 530      Y(I) = Y(I) + R
161        CALL F (TN, Y, FTEM, RPAR, IPAR)
162        DO 550 JJ = J,N,MBAND
163          Y(JJ) = YH(JJ,1)
164          YJJ = Y(JJ)
165          R = MAX(SRUR*ABS(YJJ),R0*EWT(JJ))
166          FAC = -HL0/R
167          I1 = MAX(JJ-MU,1)
168          I2 = MIN(JJ+ML,N)
169          II = JJ*MEB1 - ML + 2
170          DO 540 I = I1,I2
171 540        WM(II+I) = (FTEM(I) - SAVF(I))*FAC
172 550      CONTINUE
173 560    CONTINUE
174      NFE = NFE + MBA
175C ADD IDENTITY MATRIX. -------------------------------------------------
176 570  II = MBAND + 2
177      DO 580 I = 1,N
178        WM(II) = WM(II) + 1.0E0
179 580    II = II + MEBAND
180C DO LU DECOMPOSITION OF P. --------------------------------------------
181      CALL SGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
182      RETURN
183C----------------------- END OF SUBROUTINE PJAC -----------------------
184      END
185