subroutine sprja (neq, y, yh, nyh, ewt, ftem, savf, wm, iwm, 1 f, jac) external f, jac integer neq, nyh, iwm integer icf, ierpj, iersl, jcur, jstart, kflag, l, 1 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter, 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu integer jtyp, mused, mxordn, mxords integer i, i1, i2, ier, ii, j, j1, jj, lenp, 1 mba, mband, meb1, meband, ml, ml3, mu, np1 real y, yh, ewt, ftem, savf, wm real ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround real pdnorm real con, fac, hl0, r, r0, srur, yi, yj, yjj, 1 smnorm, sfnorm, sbnorm dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*), 1 wm(*), iwm(*) common /sls001/ ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, 1 icf, ierpj, iersl, jcur, jstart, kflag, l, 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter, 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu common /slsa01/ pdnorm, jtyp, mused, mxordn, mxords C----------------------------------------------------------------------- C SPRJA is called by SSTODA to compute and process the matrix C P = I - H*EL(1)*J , where J is an approximation to the Jacobian. C Here J is computed by the user-supplied routine JAC if C MITER = 1 or 4 or by finite differencing if MITER = 2 or 5. C J, scaled by -H*EL(1), is stored in WM. Then the norm of J (the C matrix norm consistent with the weighted max-norm on vectors given C by SMNORM) is computed, and J is overwritten by P. P is then C subjected to LU decomposition in preparation for later solution C of linear systems with P as coefficient matrix. This is done C by SGEFA if MITER = 1 or 2, and by SGBFA if MITER = 4 or 5. C C In addition to variables described previously, communication C with SPRJA uses the following: C Y = array containing predicted values on entry. C FTEM = work array of length N (ACOR in SSTODA). C SAVF = array containing f evaluated at predicted y. C WM = real work space for matrices. On output it contains the C LU decomposition of P. C Storage of matrix elements starts at WM(3). C WM also contains the following matrix-related data: C WM(1) = SQRT(UROUND), used in numerical Jacobian increments. C IWM = integer work space containing pivot information, starting at C IWM(21). IWM also contains the band parameters C ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5. C EL0 = EL(1) (input). C PDNORM= norm of Jacobian matrix. (Output). C IERPJ = output error flag, = 0 if no trouble, .gt. 0 if C P matrix found to be singular. C JCUR = output flag = 1 to indicate that the Jacobian matrix C (or approximation) is now current. C This routine also uses the Common variables EL0, H, TN, UROUND, C MITER, N, NFE, and NJE. C----------------------------------------------------------------------- nje = nje + 1 ierpj = 0 jcur = 1 hl0 = h*el0 go to (100, 200, 300, 400, 500), miter c if miter = 1, call jac and multiply by scalar. ----------------------- 100 lenp = n*n do 110 i = 1,lenp 110 wm(i+2) = 0.0e0 call jac (neq, tn, y, 0, 0, wm(3), n) con = -hl0 do 120 i = 1,lenp 120 wm(i+2) = wm(i+2)*con go to 240 c if miter = 2, make n calls to f to approximate j. -------------------- 200 fac = smnorm (n, savf, ewt) r0 = 1000.0e0*abs(h)*uround*n*fac if (r0 .eq. 0.0e0) r0 = 1.0e0 srur = wm(1) j1 = 2 do 230 j = 1,n yj = y(j) r = max(srur*abs(yj),r0/ewt(j)) y(j) = y(j) + r fac = -hl0/r call f (neq, tn, y, ftem) do 220 i = 1,n 220 wm(i+j1) = (ftem(i) - savf(i))*fac y(j) = yj j1 = j1 + n 230 continue nfe = nfe + n 240 continue c compute norm of jacobian. -------------------------------------------- pdnorm = sfnorm (n, wm(3), ewt)/abs(hl0) c add identity matrix. ------------------------------------------------- j = 3 np1 = n + 1 do 250 i = 1,n wm(j) = wm(j) + 1.0e0 250 j = j + np1 c do lu decomposition on p. -------------------------------------------- call sgefa (wm(3), n, n, iwm(21), ier) if (ier .ne. 0) ierpj = 1 return c dummy block only, since miter is never 3 in this routine. ------------ 300 return c if miter = 4, call jac and multiply by scalar. ----------------------- 400 ml = iwm(1) mu = iwm(2) ml3 = ml + 3 mband = ml + mu + 1 meband = mband + ml lenp = meband*n do 410 i = 1,lenp 410 wm(i+2) = 0.0e0 call jac (neq, tn, y, ml, mu, wm(ml3), meband) con = -hl0 do 420 i = 1,lenp 420 wm(i+2) = wm(i+2)*con go to 570 c if miter = 5, make mband calls to f to approximate j. ---------------- 500 ml = iwm(1) mu = iwm(2) mband = ml + mu + 1 mba = min(mband,n) meband = mband + ml meb1 = meband - 1 srur = wm(1) fac = smnorm (n, savf, ewt) r0 = 1000.0e0*abs(h)*uround*n*fac if (r0 .eq. 0.0e0) r0 = 1.0e0 do 560 j = 1,mba do 530 i = j,n,mband yi = y(i) r = max(srur*abs(yi),r0/ewt(i)) 530 y(i) = y(i) + r call f (neq, tn, y, ftem) do 550 jj = j,n,mband y(jj) = yh(jj,1) yjj = y(jj) r = max(srur*abs(yjj),r0/ewt(jj)) fac = -hl0/r i1 = max(jj-mu,1) i2 = min(jj+ml,n) ii = jj*meb1 - ml + 2 do 540 i = i1,i2 540 wm(ii+i) = (ftem(i) - savf(i))*fac 550 continue 560 continue nfe = nfe + mba 570 continue c compute norm of jacobian. -------------------------------------------- pdnorm = sbnorm (n, wm(ml+3), meband, ml, mu, ewt)/abs(hl0) c add identity matrix. ------------------------------------------------- ii = mband + 2 do 580 i = 1,n wm(ii) = wm(ii) + 1.0e0 580 ii = ii + meband c do lu decomposition of p. -------------------------------------------- call sgbfa (wm(3), meband, n, ml, mu, iwm(21), ier) if (ier .ne. 0) ierpj = 1 return c----------------------- end of subroutine sprja ----------------------- end