1C/MEMBR ADD NAME=PREPJI,SSI=0 2 subroutine prepji (neq, y, yh, nyh, ewt, rtem, savr, s, wm, iwm, 3 1 res, jac, adda) 4clll. optimize 5 6 external res, jac, adda 7 integer neq, nyh, iwm 8 integer iownd, iowns, 9 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 10 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu 11 integer i, i1, i2, ier, ii, ires, j, j1, jj, lenp, 12 1 mba, mband, meb1, meband, ml, ml3, mu 13 double precision y, yh, ewt, rtem, savr, s, wm 14 double precision rownd, rowns, 15 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround 16 double precision con, fac, hl0, r, srur, yi, yj, yjj 17 dimension neq(*), y(*), yh(nyh,*), ewt(*), rtem(*), 18 1 s(*), savr(*), wm(*), iwm(*) 19cDEC$ ATTRIBUTES DLLIMPORT:: /ls0001/ 20 common /ls0001/ rownd, rowns(209), 21 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, 22 3 iownd(14), iowns(6), 23 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 24 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu 25c----------------------------------------------------------------------- 26c%purpose 27c prepji is called by stodi to compute and process the matrix 28c p = a - h*el(1)*j , where j is an approximation to the jacobian dr/dy, 29c where r = g(t,y) - a(t,y)*s. here j is computed by the user-supplied 30c routine jac if miter = 1 or 4, or by finite differencing if miter = 31c 2 or 5. j is stored in wm, rescaled, and adda is called to generate 32c p. p is then subjected to lu decomposition in preparation 33c for later solution of linear systems with p as coefficient 34c matrix. this is done by dgefa if miter = 1 or 2, and by 35c dgbfa if miter = 4 or 5. 36c 37c%additional parameters 38c in addition to variables described previously, communication 39c with prepji uses the following.. 40c y = array containing predicted values on entry. 41c rtem = work array of length n (acor in stodi). 42c savr = array used for output only. on output it contains the 43c residual evaluated at current values of t and y. 44c s = array containing predicted values of dy/dt (savf in stodi). 45c wm = real work space for matrices. on output it contains the 46c lu decomposition of p. 47c storage of matrix elements starts at wm(3). 48c wm also contains the following matrix-related data.. 49c wm(1) = sqrt(uround), used in numerical jacobian increments. 50c iwm = integer work space containing pivot information, starting at 51c iwm(21). iwm also contains the band parameters 52c ml = iwm(1) and mu = iwm(2) if miter is 4 or 5. 53c el0 = el(1) (input). 54c ierpj = output error flag. 55c = 0 if no trouble occurred, 56c = 1 if the p matrix was found to be singular, 57c = ires (= 2 or 3) if res returned ires = 2 or 3. 58c jcur = output flag = 1 to indicate that the jacobian matrix 59c (or approximation) is now current. 60c this routine also uses the common variables el0, h, tn, uround, 61c miter, n, nre, and nje. 62c! 63c----------------------------------------------------------------------- 64 nje = nje + 1 65 hl0 = h*el0 66 ierpj = 0 67 jcur = 1 68 go to (100, 200, 300, 400, 500), miter 69c if miter = 1, call res, then jac, and multiply by scalar. ------------ 70 100 ires = 1 71 call res (neq, tn, y, s, savr, ires) 72 if(ierror.gt.0) return 73 nre = nre + 1 74 if (ires .gt. 1) go to 600 75 lenp = n*n 76 do 110 i = 1,lenp 77 110 wm(i+2) = 0.0d+0 78 call jac ( neq, tn, y, s, 0, 0, wm(3), n ) 79 if(ierror.gt.0) return 80 con = -hl0 81 do 120 i = 1,lenp 82 120 wm(i+2) = wm(i+2)*con 83 go to 240 84c if miter = 2, make n + 1 calls to res to approximate j. -------------- 85 200 continue 86 ires = -1 87 call res (neq, tn, y, s, savr, ires) 88 if(ierror.gt.0) return 89 nre = nre + 1 90 if (ires .gt. 1) go to 600 91 srur = wm(1) 92 j1 = 2 93 do 230 j = 1,n 94 yj = y(j) 95 r = max(srur*abs(yj),0.010d+0/ewt(j)) 96 y(j) = y(j) + r 97 fac = -hl0/r 98 call res ( neq, tn, y, s, rtem, ires ) 99 if(ierror.gt.0) return 100 nre = nre + 1 101 if (ires .gt. 1) go to 600 102 do 220 i = 1,n 103 220 wm(i+j1) = (rtem(i) - savr(i))*fac 104 y(j) = yj 105 j1 = j1 + n 106 230 continue 107 ires = 1 108 call res (neq, tn, y, s, savr, ires) 109 if(ierror.gt.0) return 110 nre = nre + 1 111 if (ires .gt. 1) go to 600 112c add matrix a. -------------------------------------------------------- 113 240 continue 114 call adda(neq, tn, y, 0, 0, wm(3), n) 115 if(ierror.gt.0) return 116c do lu decomposition on p. -------------------------------------------- 117 call dgefa (wm(3), n, n, iwm(21), ier) 118 if (ier .ne. 0) ierpj = 1 119 return 120c dummy section for miter = 3 121 300 return 122c if miter = 4, call res, then jac, and multiply by scalar. ------------ 123 400 ires = 1 124 call res (neq, tn, y, s, savr, ires) 125 if(ierror.gt.0) return 126 nre = nre + 1 127 if (ires .gt. 1) go to 600 128 ml = iwm(1) 129 mu = iwm(2) 130cc mod 06-01-89 131cc ml3 = ml + 3 132 ml3 = 3 133cc fin 134 mband = ml + mu + 1 135 meband = mband + ml 136 lenp = meband*n 137 do 410 i = 1,lenp 138 410 wm(i+2) = 0.0d+0 139 call jac ( neq, tn, y, s, ml, mu, wm(ml3), meband) 140 if(ierror.gt.0) return 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 ml + mu + 2 calls to res to approximate j. -------- 146 500 continue 147 ires = -1 148 call res (neq, tn, y, s, savr, ires) 149 if(ierror.gt.0) return 150 nre = nre + 1 151 if (ires .gt. 1) go to 600 152 ml = iwm(1) 153 mu = iwm(2) 154 ml3 = ml + 3 155 mband = ml + mu + 1 156 mba = min(mband,n) 157 meband = mband + ml 158 meb1 = meband - 1 159 srur = wm(1) 160 do 560 j = 1,mba 161 do 530 i = j,n,mband 162 yi = y(i) 163 r = max(srur*abs(yi),0.010d+0/ewt(i)) 164 530 y(i) = y(i) + r 165 call res ( neq, tn, y, s, rtem, ires) 166 if(ierror.gt.0) return 167 nre = nre + 1 168 if (ires .gt. 1) go to 600 169 do 550 jj = j,n,mband 170 y(jj) = yh(jj,1) 171 yjj = y(jj) 172 r = max(srur*abs(yjj),0.010d+0/ewt(jj)) 173 fac = -hl0/r 174 i1 = max(jj-mu,1) 175 i2 = min(jj+ml,n) 176 ii = jj*meb1 - ml + 2 177 do 540 i = i1,i2 178 540 wm(ii+i) = (rtem(i) - savr(i))*fac 179 550 continue 180 560 continue 181 ires = 1 182 call res (neq, tn, y, s, savr, ires) 183 if(ierror.gt.0) return 184 nre = nre + 1 185 if (ires .gt. 1) go to 600 186c add matrix a. -------------------------------------------------------- 187 570 continue 188 call adda(neq, tn, y, ml, mu, wm(ml3), meband) 189 if(ierror.gt.0) return 190c do lu decomposition of p. -------------------------------------------- 191 call dgbfa (wm(3), meband, n, ml, mu, iwm(21), ier) 192 if (ier .ne. 0) ierpj = 1 193 return 194c error return for ires = 2 or ires = 3 return from res. --------------- 195 600 ierpj = ires 196 return 197c----------------------- end of subroutine prepji ---------------------- 198 end 199