2 DOUBLE PRECISION FUNCTION dumach ()
30 DOUBLE PRECISION u, comp
34 CALL dumsum(1.0d0, u, comp)
35 IF (comp .NE. 1.0d0)
GO TO 10
40 SUBROUTINE dumsum(A,B,C)
42 DOUBLE PRECISION a, b, c
47 SUBROUTINE dcfode (METH, ELCO, TESCO)
90 INTEGER i, ib, nq, nqm1, nqp1
91 DOUBLE PRECISION elco, tesco
92 DOUBLE PRECISION agamq, fnq, fnqm1, pc, pint, ragq,
93 1 rqfac, rq1fac, tsign, xpin
94 dimension elco(13,12), tesco(3,12)
98 GO TO (100, 200), meth
100 100 elco(1,1) = 1.0d0
123 110 pc(i) = pc(i-1) + fnqm1*pc(i)
131 pint = pint + tsign*pc(i)/i
132 120 xpin = xpin + tsign*pc(i)/(i+1)
134 elco(1,nq) = pint*rq1fac
137 130 elco(i+1,nq) = rq1fac*pc(i)/i
141 IF (nq .LT. 12) tesco(1,nqp1) = ragq*rqfac/nqp1
160 210 pc(i) = pc(i-1) + fnq*pc(i)
164 220 elco(i,nq) = pc(i)/pc(2)
167 tesco(2,nq) = nqp1/elco(1,nq)
168 tesco(3,nq) = (nq+2)/elco(1,nq)
175 SUBROUTINE dintdy (T, K, YH, NYH, DKY, IFLAG)
216 INTEGER k, nyh, iflag
217 DOUBLE PRECISION t, yh, dky
218 dimension yh(nyh,*), dky(*)
219 INTEGER iownd, iowns,
220 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
221 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
222 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
223 DOUBLE PRECISION rowns,
224 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
225 COMMON /dls001/ rowns(209),
226 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
227 2 iownd(6), iowns(6),
228 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
229 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
230 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
231 INTEGER i, ic, j, jb, jb2, jj, jj1, jp1
232 DOUBLE PRECISION c, r, s, tp
237 IF (k .LT. 0 .OR. k .GT. nq)
GO TO 80
238 tp = tn - hu - 100.0d0*uround*sign(abs(tn) + abs(hu), hu)
239 IF ((t-tp)*(t-tn) .GT. 0.0d0)
GO TO 90
243 IF (k .EQ. 0)
GO TO 15
249 20 dky(i) = c*yh(i,l)
250 IF (k .EQ. nq)
GO TO 55
256 IF (k .EQ. 0)
GO TO 35
262 40 dky(i) = c*yh(i,jp1) + s*dky(i)
270 80 msg =
'DINTDY- K (=I1) illegal ' 271 CALL xerrwd (msg, 30, 51, 0, 1, k, 0, 0, 0.0d0, 0.0d0)
274 90 msg =
'DINTDY- T (=R1) illegal ' 275 CALL xerrwd (msg, 30, 52, 0, 0, 0, 0, 1, t, 0.0d0)
276 msg=
' T not in interval TCUR - HU (= R1) to TCUR (=R2) ' 277 CALL xerrwd (msg, 60, 52, 0, 0, 0, 0, 2, tp, tn)
283 SUBROUTINE dprepj (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
339 INTEGER neq, nyh, iwm
340 DOUBLE PRECISION y, yh, ewt, ftem, savf, wm
341 dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
343 INTEGER iownd, iowns,
344 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
345 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
346 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
347 DOUBLE PRECISION rowns,
348 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
349 COMMON /dls001/ rowns(209),
350 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
351 2 iownd(6), iowns(6),
352 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
353 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
354 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
355 INTEGER i, i1, i2, ier, ii, j, j1, jj, lenp,
356 1 mba, mband, meb1, meband, ml, ml3, mu, np1
357 DOUBLE PRECISION con, di, fac, hl0, r, r0, srur, yi, yj, yjj,
365 GO TO (100, 200, 300, 400, 500), miter
370 CALL jac (neq, tn, y, 0, 0, wm(3), n)
373 120 wm(i+2) = wm(i+2)*con
376 200 fac = dvnorm(n, savf, ewt)
377 r0 = 1000.0d0*abs(h)*uround*n*fac
378 IF (r0 .EQ. 0.0d0) r0 = 1.0d0
383 r = max(srur*abs(yj),r0/ewt(j))
386 CALL f (neq, tn, y, ftem)
388 220 wm(i+j1) = (ftem(i) - savf(i))*fac
397 wm(j) = wm(j) + 1.0d0
400 CALL dgefa (wm(3), n, n, iwm(21), ier)
401 IF (ier .NE. 0) ierpj = 1
407 310 y(i) = y(i) + r*(h*savf(i) - yh(i,2))
408 CALL f (neq, tn, y, wm(3))
411 r0 = h*savf(i) - yh(i,2)
412 di = 0.1d0*r0 - h*(wm(i+2) - savf(i))
414 IF (abs(r0) .LT. uround/ewt(i))
GO TO 320
415 IF (abs(di) .EQ. 0.0d0)
GO TO 330
416 wm(i+2) = 0.1d0*r0/di
430 CALL jac (neq, tn, y, ml, mu, wm(ml3), meband)
433 420 wm(i+2) = wm(i+2)*con
443 fac = dvnorm(n, savf, ewt)
444 r0 = 1000.0d0*abs(h)*uround*n*fac
445 IF (r0 .EQ. 0.0d0) r0 = 1.0d0
449 r = max(srur*abs(yi),r0/ewt(i))
451 CALL f (neq, tn, y, ftem)
452 DO 550 jj = j,n,mband
455 r = max(srur*abs(yjj),r0/ewt(jj))
459 ii = jj*meb1 - ml + 2
461 540 wm(ii+i) = (ftem(i) - savf(i))*fac
468 wm(ii) = wm(ii) + 1.0d0
471 CALL dgbfa (wm(3), meband, n, ml, mu, iwm(21), ier)
472 IF (ier .NE. 0) ierpj = 1
477 SUBROUTINE dsolsy (WM, IWM, X, TEM)
522 DOUBLE PRECISION wm, x, tem
523 dimension wm(*), iwm(*), x(*), tem(*)
524 INTEGER iownd, iowns,
525 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
526 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
527 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
528 DOUBLE PRECISION rowns,
529 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
530 COMMON /dls001/ rowns(209),
531 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
532 2 iownd(6), iowns(6),
533 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
534 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
535 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
536 INTEGER i, meband, ml, mu
537 DOUBLE PRECISION di, hl0, phl0, r
541 GO TO (100, 100, 300, 400, 400), miter
542 100
CALL dgesl (wm(3), n, n, iwm(21), x, 0)
548 IF (hl0 .EQ. phl0)
GO TO 330
551 di = 1.0d0 - r*(1.0d0 - 1.0d0/wm(i+2))
552 IF (abs(di) .EQ. 0.0d0)
GO TO 390
553 320 wm(i+2) = 1.0d0/di
555 340 x(i) = wm(i+2)*x(i)
562 meband = 2*ml + mu + 1
563 CALL dgbsl (wm(3), meband, n, ml, mu, iwm(21), x, 0)
568 SUBROUTINE dsrcom (RSAV, ISAV, JOB)
605 INTEGER i, lenils, lenrls
606 DOUBLE PRECISION rsav, rls
607 dimension rsav(*), isav(*)
609 COMMON /dls001/ rls(218), ils(37)
610 DATA lenrls/218/, lenils/37/
613 IF (job .EQ. 2)
GO TO 100
630 SUBROUTINE dstode (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
631 1 WM, IWM, F, JAC, PJAC, SLVS)
724 EXTERNAL f, jac, pjac, slvs
725 INTEGER neq, nyh, iwm
726 DOUBLE PRECISION y, yh, yh1, ewt, savf, acor, wm
727 dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
728 1 acor(*), wm(*), iwm(*)
729 INTEGER iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
730 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
731 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
732 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
733 INTEGER i, i1, iredo, iret, j, jb, m, ncf, newq
734 DOUBLE PRECISION conit, crate, el, elco, hold, rmax, tesco,
735 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
736 DOUBLE PRECISION dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup,
737 1 r, rh, rhdn, rhsm, rhup, told, dvnorm
738 COMMON /dls001/ conit, crate, el(13), elco(13,12),
739 1 hold, rmax, tesco(3,12),
740 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
741 3 iownd(6), ialth, ipup, lmax, meo, nqnyh, nslp,
742 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
743 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
744 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
755 IF (jstart .GT. 0)
GO TO 200
756 IF (jstart .EQ. -1)
GO TO 100
757 IF (jstart .EQ. -2)
GO TO 160
795 IF (ialth .EQ. 1) ialth = 2
796 IF (meth .EQ. meo)
GO TO 110
797 CALL dcfode (meth, elco, tesco)
799 IF (nq .GT. maxord)
GO TO 120
803 110
IF (nq .LE. maxord)
GO TO 160
807 125 el(i) = elco(i,nq)
812 ddn = dvnorm(n, savf, ewt)/tesco(1,l)
814 rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
817 IF (h .EQ. hold)
GO TO 170
818 rh = min(rh,abs(h/hold))
826 140
CALL dcfode (meth, elco, tesco)
828 155 el(i) = elco(i,nq)
833 GO TO (160, 170, 200), iret
840 160
IF (h .EQ. hold)
GO TO 200
845 170 rh = max(rh,hmin/abs(h))
846 175 rh = min(rh,rmax)
847 rh = rh/max(1.0d0,abs(h)*hmxi*rh)
852 180 yh(i,j) = yh(i,j)*r
856 IF (iredo .EQ. 0)
GO TO 690
865 200
IF (abs(rc-1.0d0) .GT. ccmax) ipup = miter
866 IF (nst .GE. nslp+msbp) ipup = miter
873 210 yh1(i) = yh1(i) + yh1(i+nyh)
884 CALL f (neq, tn, y, savf)
886 IF (ipup .LE. 0)
GO TO 250
892 CALL pjac (neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac)
897 IF (ierpj .NE. 0)
GO TO 430
900 270
IF (miter .NE. 0)
GO TO 350
906 savf(i) = h*savf(i) - yh(i,2)
907 290 y(i) = savf(i) - acor(i)
908 del = dvnorm(n, y, ewt)
910 y(i) = yh(i,1) + el(1)*savf(i)
911 300 acor(i) = savf(i)
919 360 y(i) = h*savf(i) - (yh(i,2) + acor(i))
920 CALL slvs (wm, iwm, y, savf)
921 IF (iersl .LT. 0)
GO TO 430
922 IF (iersl .GT. 0)
GO TO 410
923 del = dvnorm(n, y, ewt)
925 acor(i) = acor(i) + y(i)
926 380 y(i) = yh(i,1) + el(1)*acor(i)
931 400
IF (m .NE. 0) crate = max(0.2d0*crate,del/delp)
932 dcon = del*min(1.0d0,1.5d0*crate)/(tesco(2,nq)*conit)
933 IF (dcon .LE. 1.0d0)
GO TO 450
935 IF (m .EQ. maxcor)
GO TO 410
936 IF (m .GE. 2 .AND. del .GT. 2.0d0*delp)
GO TO 410
938 CALL f (neq, tn, y, savf)
948 410
IF (miter .EQ. 0 .OR. jcur .EQ. 1)
GO TO 430
961 440 yh1(i) = yh1(i) - yh1(i+nyh)
963 IF (ierpj .LT. 0 .OR. iersl .LT. 0)
GO TO 680
964 IF (abs(h) .LE. hmin*1.00001d0)
GO TO 670
965 IF (ncf .EQ. mxncf)
GO TO 670
977 IF (m .EQ. 0) dsm = del/tesco(2,nq)
978 IF (m .GT. 0) dsm = dvnorm(n, acor, ewt)/tesco(2,nq)
979 IF (dsm .GT. 1.0d0)
GO TO 500
997 470 yh(i,j) = yh(i,j) + el(j)*acor(i)
999 IF (ialth .EQ. 0)
GO TO 520
1000 IF (ialth .GT. 1)
GO TO 700
1001 IF (l .EQ. lmax)
GO TO 700
1003 490 yh(i,lmax) = acor(i)
1012 500 kflag = kflag - 1
1019 510 yh1(i) = yh1(i) - yh1(i+nyh)
1022 IF (abs(h) .LE. hmin*1.00001d0)
GO TO 660
1023 IF (kflag .LE. -3)
GO TO 640
1037 IF (l .EQ. lmax)
GO TO 540
1039 530 savf(i) = acor(i) - yh(i,lmax)
1040 dup = dvnorm(n, savf, ewt)/tesco(3,nq)
1042 rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0)
1044 rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
1046 IF (nq .EQ. 1)
GO TO 560
1047 ddn = dvnorm(n, yh(1,l), ewt)/tesco(1,nq)
1049 rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
1050 560
IF (rhsm .GE. rhup)
GO TO 570
1051 IF (rhup .GT. rhdn)
GO TO 590
1053 570
IF (rhsm .LT. rhdn)
GO TO 580
1059 IF (kflag .LT. 0 .AND. rh .GT. 1.0d0) rh = 1.0d0
1063 IF (rh .LT. 1.1d0)
GO TO 610
1066 600 yh(i,newq+1) = acor(i)*r
1070 620
IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1d0))
GO TO 610
1071 IF (kflag .LE. -2) rh = min(rh,0.2d0)
1077 IF (newq .EQ. nq)
GO TO 170
1091 640
IF (kflag .EQ. -10)
GO TO 660
1093 rh = max(hmin/abs(h),rh)
1097 CALL f (neq, tn, y, savf)
1100 650 yh(i,2) = h*savf(i)
1103 IF (nq .EQ. 1)
GO TO 200
1119 700 r = 1.0d0/tesco(2,nqu)
1121 710 acor(i) = acor(i)*r
1128 SUBROUTINE dewset (N, ITOL, RTOL, ATOL, YCUR, EWT)
1152 DOUBLE PRECISION rtol, atol, ycur, ewt
1153 dimension rtol(*), atol(*), ycur(n), ewt(n)
1156 GO TO (10, 20, 30, 40), itol
1159 15 ewt(i) = rtol(1)*abs(ycur(i)) + atol(1)
1163 25 ewt(i) = rtol(1)*abs(ycur(i)) + atol(i)
1167 35 ewt(i) = rtol(i)*abs(ycur(i)) + atol(1)
1171 45 ewt(i) = rtol(i)*abs(ycur(i)) + atol(i)
1176 DOUBLE PRECISION FUNCTION dvnorm (N, V, W)
1199 DOUBLE PRECISION v, w, sum
1200 dimension v(n), w(n)
1205 10 sum = sum + (v(i)*w(i))**2
1206 dvnorm = sqrt(sum/n)
1211 SUBROUTINE diprep (NEQ, Y, RWORK, IA, JA, IPFLAG, F, JAC)
1213 INTEGER neq, ia, ja, ipflag
1214 DOUBLE PRECISION y, rwork
1215 dimension neq(*), y(*), rwork(*), ia(*), ja(*)
1216 INTEGER iownd, iowns,
1217 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
1218 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1219 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1220 INTEGER iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
1221 1 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
1222 2 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
1223 3 nslj, ngp, nlu, nnz, nsp, nzl, nzu
1224 DOUBLE PRECISION rowns,
1225 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1226 DOUBLE PRECISION rlss
1227 COMMON /dls001/ rowns(209),
1228 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1229 2 iownd(6), iowns(6),
1230 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
1231 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1232 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1233 COMMON /dlss01/ rlss(6),
1234 1 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
1235 2 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
1236 3 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
1237 4 nslj, ngp, nlu, nnz, nsp, nzl, nzu
1238 INTEGER i, imax, lewtn, lyhd, lyhn
1254 CALL dprep (neq, y, rwork(lyh), rwork(lsavf), rwork(lewt),
1255 1 rwork(lacor), ia, ja, rwork(lwm), rwork(lwm), ipflag, f, jac)
1256 lenwk = max(lreq,lwmin)
1257 IF (ipflag .LT. 0)
RETURN 1260 IF (lyhn .GT. lyh)
RETURN 1262 IF (lyhd .EQ. 0)
GO TO 20
1263 imax = lyhn - 1 + lenyhm
1265 10 rwork(i) = rwork(i+lyhd)
1268 20 lsavf = lyh + lenyh
1271 IF (istatc .EQ. 3)
GO TO 40
1273 IF (lewtn .GT. lewt)
RETURN 1275 30 rwork(i+lewtn-1) = rwork(i+lewt-1)
1281 SUBROUTINE dprep (NEQ, Y, YH, SAVF, EWT, FTEM, IA, JA,
1282 1 WK, IWK, IPPER, F, JAC)
1284 INTEGER neq, ia, ja, iwk, ipper
1285 DOUBLE PRECISION y, yh, savf, ewt, ftem, wk
1286 dimension neq(*), y(*), yh(*), savf(*), ewt(*), ftem(*),
1287 1 ia(*), ja(*), wk(*), iwk(*)
1288 INTEGER iownd, iowns,
1289 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
1290 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1291 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1292 INTEGER iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
1293 1 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
1294 2 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
1295 3 nslj, ngp, nlu, nnz, nsp, nzl, nzu
1296 DOUBLE PRECISION rowns,
1297 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1298 DOUBLE PRECISION con0, conmin, ccmxj, psmall, rbig, seth
1299 COMMON /dls001/ rowns(209),
1300 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1301 2 iownd(6), iowns(6),
1302 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
1303 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1304 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1305 COMMON /dlss01/ con0, conmin, ccmxj, psmall, rbig, seth,
1306 1 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
1307 2 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
1308 3 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
1309 4 nslj, ngp, nlu, nnz, nsp, nzl, nzu
1310 INTEGER i, ibr, ier, ipil, ipiu, iptt1, iptt2, j, jfound, k,
1311 1 knew, kmax, kmin, ldif, lenigp, liwk, maxg, np1, nzsut
1312 DOUBLE PRECISION dq, dyj, erwt, fac, yj
1354 IF (ipjan+n-1 .GT. liwk)
GO TO 210
1355 IF (moss .EQ. 0)
GO TO 30
1357 IF (istatc .EQ. 3)
GO TO 20
1361 fac = 1.0d0 + 1.0d0/(i + 1.0d0)
1362 y(i) = y(i) + fac*sign(erwt,y(i))
1364 GO TO (70, 100), moss
1370 GO TO (70, 100), moss
1379 IF (kmin .GT. kmax)
GO TO 45
1382 IF (i .EQ. j) jfound = 1
1383 IF (knew .GT. liwk)
GO TO 210
1387 IF (jfound .EQ. 1)
GO TO 50
1388 45
IF (knew .GT. liwk)
GO TO 210
1391 50 iwk(ipian+j) = knew + 1 - ipjan
1399 CALL f (neq, tn, y, savf)
1403 IF (k .GT. liwk)
GO TO 210
1408 CALL jac (neq, tn, y, j, iwk(ipian), iwk(ipjan), savf)
1410 IF (abs(savf(i)) .LE. seth)
GO TO 80
1411 IF (i .EQ. j)
GO TO 80
1412 IF (k .GT. liwk)
GO TO 210
1416 iwk(ipian+j) = k + 1 - ipjan
1423 CALL f (neq, tn, y, savf)
1425 IF (k .GT. liwk)
GO TO 210
1432 CALL f (neq, tn, y, ftem)
1435 dq = (ftem(i) - savf(i))/dyj
1436 IF (abs(dq) .LE. seth)
GO TO 110
1437 IF (i .EQ. j)
GO TO 110
1438 IF (k .GT. liwk)
GO TO 210
1442 iwk(ipian+j) = k + 1 - ipjan
1446 IF (moss .EQ. 0 .OR. istatc .NE. 1)
GO TO 150
1450 150 nnz = iwk(ipian+n) - 1
1453 IF (miter .NE. 2)
GO TO 160
1462 lreq = iptt2 + n - 1
1463 IF (lreq .GT. liwk)
GO TO 220
1464 CALL jgroup (n, iwk(ipian), iwk(ipjan), maxg, ngp, iwk(ipigp),
1465 1 iwk(ipjgp), iwk(iptt1), iwk(iptt2), ier)
1466 IF (ier .NE. 0)
GO TO 220
1470 160 ipr = ipigp + lenigp
1474 iprsp = (ipisp - 2)/lrat + 2
1475 iesp = lenwk + 1 - iprsp
1476 IF (iesp .LT. 0)
GO TO 230
1480 nsp = liwk + 1 - ipisp
1481 CALL odrv (n, iwk(ipian), iwk(ipjan), wk, iwk(ipr), iwk(ipic),
1482 1 nsp, iwk(ipisp), 1, iys)
1483 IF (iys .EQ. 11*n+1)
GO TO 240
1484 IF (iys .NE. 0)
GO TO 230
1487 ipa = lenwk + 1 - nnz
1489 lreq = max(12*n/lrat, 6*n/lrat+2*n+nnz) + 3
1490 lreq = lreq + iprsp - 1 + nnz
1491 IF (lreq .GT. lenwk)
GO TO 250
1494 180 wk(iba+i) = 0.0d0
1495 ipisp = lrat*(iprsp - 1) + 1
1496 CALL cdrv (n,iwk(ipr),iwk(ipc),iwk(ipic),iwk(ipian),iwk(ipjan),
1497 1 wk(ipa),wk(ipa),wk(ipa),nsp,iwk(ipisp),wk(iprsp),iesp,5,iys)
1499 IF (iys .EQ. 10*n+1)
GO TO 250
1500 IF (iys .NE. 0)
GO TO 260
1502 ipiu = ipil + 2*n + 1
1503 nzu = iwk(ipil+n) - iwk(ipil)
1504 nzl = iwk(ipiu+n) - iwk(ipiu)
1505 IF (lrat .GT. 1)
GO TO 190
1506 CALL adjlr (n, iwk(ipisp), ldif)
1509 IF (lrat .EQ. 2 .AND. nnz .EQ. n) lreq = lreq + 1
1510 nsp = nsp + lreq - lenwk
1511 ipa = lreq + 1 - nnz
1517 lreq = 2 + (2*n + 1)/lrat
1518 lreq = max(lenwk+1,lreq)
1522 lreq = (lreq - 1)/lrat + 1
1526 CALL cntnzu (n, iwk(ipian), iwk(ipjan), nzsut)
1527 lreq = lenwk - iesp + (3*n + 4*nzsut - 1)/lrat + 1
1542 SUBROUTINE jgroup (N,IA,JA,MAXG,NGRP,IGP,JGP,INCL,JDONE,IER)
1543 INTEGER n, ia, ja, maxg, ngrp, igp, jgp, incl, jdone, ier
1544 dimension ia(*), ja(*), igp(*), jgp(*), incl(*), jdone(*)
1565 INTEGER i, j, k, kmin, kmax, ncol, ng
1577 IF (jdone(j) .EQ. 1)
GO TO 50
1583 IF (incl(i) .EQ. 1)
GO TO 50
1594 IF (ncol .EQ. igp(ng))
GO TO 70
1597 IF (ncol .LE. n)
GO TO 80
1606 SUBROUTINE adjlr (N, ISP, LDIF)
1607 INTEGER n, isp, ldif
1616 INTEGER ip, jlmax, jumax, lnfc, lsfc, nzlu
1623 nzlu = isp(n+1) - isp(1) + isp(ip+n+1) - isp(ip+1)
1624 lsfc = 12*n + 3 + 2*max(jlmax,jumax)
1625 lnfc = 9*n + 2 + jlmax + jumax + nzlu
1626 ldif = max(0, lsfc - lnfc)
1631 SUBROUTINE cntnzu (N, IA, JA, NZSUT)
1632 INTEGER n, ia, ja, nzsut
1633 dimension ia(*), ja(*)
1641 INTEGER ii, jj, j, jmin, jmax, k, kmin, kmax, num
1647 IF (jmin .GT. jmax)
GO TO 50
1649 IF (ja(j) - ii) 10, 40, 30
1653 IF (kmin .GT. kmax)
GO TO 30
1655 IF (ja(k) .EQ. ii)
GO TO 40
1665 SUBROUTINE dprjs (NEQ,Y,YH,NYH,EWT,FTEM,SAVF,WK,IWK,F,JAC)
1667 INTEGER neq, nyh, iwk
1668 DOUBLE PRECISION y, yh, ewt, ftem, savf, wk
1669 dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
1671 INTEGER iownd, iowns,
1672 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
1673 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1674 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1675 INTEGER iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
1676 1 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
1677 2 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
1678 3 nslj, ngp, nlu, nnz, nsp, nzl, nzu
1679 DOUBLE PRECISION rowns,
1680 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1681 DOUBLE PRECISION con0, conmin, ccmxj, psmall, rbig, seth
1682 COMMON /dls001/ rowns(209),
1683 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1684 2 iownd(6), iowns(6),
1685 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
1686 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1687 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1688 COMMON /dlss01/ con0, conmin, ccmxj, psmall, rbig, seth,
1689 1 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
1690 2 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
1691 3 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
1692 4 nslj, ngp, nlu, nnz, nsp, nzl, nzu
1693 INTEGER i, imul, j, jj, jok, jmax, jmin, k, kmax, kmin, ng
1694 DOUBLE PRECISION con, di, fac, hl0, pij, r, r0, rcon, rcont,
1737 IF (miter .EQ. 3)
GO TO 300
1740 IF (nst .EQ. 0 .OR. nst .GE. nslj+msbj) jok = 0
1741 IF (icf .EQ. 1 .AND. abs(rc - 1.0d0) .LT. ccmxj) jok = 0
1742 IF (icf .EQ. 2) jok = 0
1743 IF (jok .EQ. 1)
GO TO 250
1751 GO TO (100, 200), miter
1757 kmax = iwk(ipian+j) - 1
1760 CALL jac (neq, tn, y, j, iwk(ipian), iwk(ipjan), ftem)
1761 DO 120 k = kmin, kmax
1763 wk(iba+k) = ftem(i)*con
1764 IF (i .EQ. j) wk(iba+k) = wk(iba+k) + 1.0d0
1772 fac = dvnorm(n, savf, ewt)
1773 r0 = 1000.0d0 * abs(h) * uround * n * fac
1774 IF (r0 .EQ. 0.0d0) r0 = 1.0d0
1778 jmax = iwk(ipigp+ng) - 1
1779 DO 210 j = jmin,jmax
1781 r = max(srur*abs(y(jj)),r0/ewt(jj))
1782 210 y(jj) = y(jj) + r
1783 CALL f (neq, tn, y, ftem)
1784 DO 230 j = jmin,jmax
1787 r = max(srur*abs(y(jj)),r0/ewt(jj))
1790 kmax =iwk(ibian+jj+1) - 1
1791 DO 220 k = kmin,kmax
1793 wk(iba+k) = (ftem(i) - savf(i))*fac
1794 IF (i .EQ. jj) wk(iba+k) = wk(iba+k) + 1.0d0
1805 rcont = abs(con)/conmin
1806 IF (rcont .GT. rbig .AND. iplost .EQ. 1)
GO TO 20
1809 kmax = iwk(ipian+j) - 1
1810 DO 270 k = kmin,kmax
1813 IF (i .NE. j)
GO TO 260
1815 IF (abs(pij) .GE. psmall)
GO TO 260
1817 conmin = min(abs(con0),conmin)
1819 IF (i .EQ. j) pij = pij + 1.0d0
1831 CALL cdrv (n,iwk(ipr),iwk(ipc),iwk(ipic),iwk(ipian),iwk(ipjan),
1832 1 wk(ipa),ftem,ftem,nsp,iwk(ipisp),wk(iprsp),iesp,2,iys)
1833 IF (iys .EQ. 0)
RETURN 1836 IF (imul .EQ. 8) ierpj = 1
1837 IF (imul .EQ. 10) ierpj = -1
1848 310 y(i) = y(i) + r*(h*savf(i) - yh(i,2))
1849 CALL f (neq, tn, y, wk(3))
1852 r0 = h*savf(i) - yh(i,2)
1853 di = 0.1d0*r0 - h*(wk(i+2) - savf(i))
1855 IF (abs(r0) .LT. uround/ewt(i))
GO TO 320
1856 IF (abs(di) .EQ. 0.0d0)
GO TO 330
1857 wk(i+2) = 0.1d0*r0/di
1865 SUBROUTINE dsolss (WK, IWK, X, TEM)
1867 DOUBLE PRECISION wk, x, tem
1868 dimension wk(*), iwk(*), x(*), tem(*)
1869 INTEGER iownd, iowns,
1870 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
1871 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1872 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1873 INTEGER iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
1874 1 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
1875 2 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
1876 3 nslj, ngp, nlu, nnz, nsp, nzl, nzu
1877 DOUBLE PRECISION rowns,
1878 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1879 DOUBLE PRECISION rlss
1880 COMMON /dls001/ rowns(209),
1881 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1882 2 iownd(6), iowns(6),
1883 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
1884 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
1885 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1886 COMMON /dlss01/ rlss(6),
1887 1 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
1888 2 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
1889 3 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
1890 4 nslj, ngp, nlu, nnz, nsp, nzl, nzu
1892 DOUBLE PRECISION di, hl0, phl0, r
1920 GO TO (100, 100, 300), miter
1921 100
CALL cdrv (n,iwk(ipr),iwk(ipc),iwk(ipic),iwk(ipian),iwk(ipjan),
1922 1 wk(ipa),x,x,nsp,iwk(ipisp),wk(iprsp),iesp,4,iersl)
1923 IF (iersl .NE. 0) iersl = -1
1929 IF (hl0 .EQ. phl0)
GO TO 330
1932 di = 1.0d0 - r*(1.0d0 - 1.0d0/wk(i+2))
1933 IF (abs(di) .EQ. 0.0d0)
GO TO 390
1934 320 wk(i+2) = 1.0d0/di
1936 340 x(i) = wk(i+2)*x(i)
1944 SUBROUTINE dsrcms (RSAV, ISAV, JOB)
1959 INTEGER i, lenils, leniss, lenrls, lenrss
1960 DOUBLE PRECISION rsav, rls, rlss
1961 dimension rsav(*), isav(*)
1962 SAVE lenrls, lenils, lenrss, leniss
1963 COMMON /dls001/ rls(218), ils(37)
1964 COMMON /dlss01/ rlss(6), ilss(34)
1965 DATA lenrls/218/, lenils/37/, lenrss/6/, leniss/34/
1967 IF (job .EQ. 2)
GO TO 100
1971 15 rsav(lenrls+i) = rlss(i)
1976 25 isav(lenils+i) = ilss(i)
1982 110 rls(i) = rsav(i)
1984 115 rlss(i) = rsav(lenrls+i)
1987 120 ils(i) = isav(i)
1989 125 ilss(i) = isav(lenils+i)
1996 * (n, ia,ja,a, p,ip, nsp,isp, path, flag)
2128 integer ia(*), ja(*), p(*), ip(*), isp(*), path, flag,
2129 * v, l, head, tmp, q
2131 double precision a(*)
2136 if (path.lt.1 .or. 5.lt.path)
go to 111
2139 if ((path-1) * (path-2) * (path-4) .ne. 0)
go to 1
2145 if (max.lt.n)
go to 110
2148 * (n, ia,ja, max,isp(v),isp(l), isp(head),p,ip, isp(v), flag)
2149 if (flag.ne.0)
go to 100
2152 1
if ((path-2) * (path-3) * (path-4) * (path-5) .ne. 0)
go to 2
2154 q = tmp - (ia(n+1)-1)
2155 if (q.lt.1)
go to 110
2157 dflag = path.eq.4 .or. path.eq.5
2159 * (n, ip, ia, ja, a, isp(tmp), isp(q), dflag)
2173 * (n, ia,ja, max, v,l, head,last,next, mark, flag)
2257 integer ia(*), ja(*), v(*), l(*), head(*), last(*), next(*),
2258 * mark(*), flag, tag, dmin, vk,ek, tail
2264 * (n, ia,ja, max,v,l, head,last,next, mark,tag, flag)
2265 if (flag.ne.0)
return 2271 1
if (k.ge.n)
go to 4
2274 2
if (head(dmin).gt.0)
go to 3
2280 head(dmin) = next(vk)
2281 if (head(dmin).gt.0) last(head(dmin)) = -dmin
2287 tag = tag + last(ek)
2292 * (vk,tail, v,l, last,next, mark)
2296 * (k,ek,tail, v,l, head,last,next, mark)
2300 * (ek,dmin, v,l, head,last,next, mark)
2312 * (n, ia,ja, max,v,l, head,last,next, mark,tag, flag)
2316 integer ia(*), ja(*), v(*), l(*), head(*), last(*), next(*),
2317 * mark(*), tag, flag, sfs, vi,dvi, vj
2331 if (jmin.gt.jmax)
go to 6
2340 if (kmax .eq. 0)
go to 4
2343 if (v(lvk).eq.vj)
go to 5
2346 4
if (sfs.ge.max)
go to 101
2349 mark(vi) = mark(vi) + 1
2356 mark(vj) = mark(vj) + 1
2367 next(vi) = head(dvi)
2371 if (nextvi.gt.0) last(nextvi) = vi
2381 * (vk,tail, v,l, last,next, mark)
2385 integer vk, tail, v(*), l(*), last(*), next(*), mark(*),
2386 * tag, s,ls,vs,es, b,lb,vb, blp,blpmax
2399 if (next(vs).lt.0)
go to 2
2419 if (mark(vb).ge.tag)
go to 3
2436 * (k,ek,tail, v,l, head,last,next, mark)
2440 integer ek, tail, v(*), l(*), head(*), last(*), next(*),
2441 * mark(*), tag, free, li,vi,lvi,evi, s,ls,es, ilp,ilpmax
2449 if (ilpmax.le.0)
go to 12
2456 if (last(vi).eq.0)
go to 3
2457 if (last(vi).gt.0)
go to 1
2458 head(-last(vi)) = next(vi)
2460 1 next(last(vi)) = next(vi)
2461 2
if (next(vi).gt.0) last(next(vi)) = last(vi)
2467 if (ls.eq.0)
go to 6
2469 if (mark(es).lt.tag)
go to 5
2477 if (lvi.ne.0)
go to 7
2483 last(ek) = last(ek) - 1
2488 7
if (l(lvi).ne.0)
go to 9
2490 if (next(evi).ge.0)
go to 9
2491 if (mark(evi).lt.0)
go to 8
2507 mark(evi) = mark(evi) - 1
2525 * (ek,dmin, v,l, head,last,next, mark)
2529 integer ek, dmin, v(*), l(*), head(*), last(*), next(*),
2530 * mark(*), tag, vi,evi,dvi, s,vs,es, b,vb, ilp,ilpmax,
2535 tag = mark(ek) - last(ek)
2540 if (ilpmax.le.0)
go to 11
2544 if (last(vi)) 1, 10, 8
2556 if (next(vs).lt.0)
go to 3
2565 3
if (mark(es).lt.0)
go to 6
2575 if (mark(vb).ge.tag)
go to 4
2585 mark(es) = mark(es) - 1
2587 if (s.eq.0)
go to 10
2589 if (mark(es).lt.0) mark(es) = mark(es) - 1
2595 dvi = last(ek) + last(evi) + mark(evi)
2599 9 next(vi) = head(dvi)
2602 if (next(vi).gt.0) last(next(vi)) = vi
2603 if (dvi.lt.dmin) dmin = dvi
2610 * (n, ip, ia,ja,a, q, r, dflag)
2641 integer ip(*), ia(*), ja(*), q(*), r(*)
2643 double precision a(*), ak
2656 if (jmin.gt.jmax)
go to 3
2661 if (ip(k).lt.ip(i)) ja(j) = i
2662 if (ip(k).ge.ip(i)) k = i
2673 ia(i+1) = ia(i) + q(i)
2682 do 6 jdummy=jmin,jmax
2684 if (.not.dflag .or. ja(j).ne.i .or. i.eq.ilast)
go to 5
2700 7
if (r(j).eq.j)
go to 8
2717 * (n, r,c,ic, ia,ja,a, b, z, nsp,isp,rsp,esp, path, flag)
2879 integer r(*), c(*), ic(*), ia(*), ja(*), isp(*), esp, path,
2880 * flag, d, u, q, row, tmp, ar, umax
2882 double precision a(*), b(*), z(*), rsp(*)
2890 if (path.lt.1 .or. 5.lt.path)
go to 111
2901 if ((path-1) * (path-5) .ne. 0)
go to 5
2902 max = (lratio*nsp + 1 - jl) - (n+1) - 5*n
2911 jumax = lratio*nsp + 1 - jutmp
2913 if (jlmax.le.0 .or. jumax.le.0)
go to 110
2916 if (c(i).ne.i)
go to 2
2921 * (n, ic, ia,ja,a, isp(il), rsp(ar), isp(iu), flag)
2922 if (flag.ne.0)
go to 100
2926 * jlmax, isp(il), isp(jl), isp(ijl),
2927 * jumax, isp(iu), isp(jutmp), isp(iju),
2928 * isp(q), isp(ira), isp(jra), isp(irac),
2929 * isp(irl), isp(jrl), isp(iru), isp(jru), flag)
2930 if(flag .ne. 0)
go to 100
2932 jlmax = isp(ijl+n-1)
2934 jumax = isp(iju+n-1)
2935 if (jumax.le.0)
go to 5
2937 4 isp(ju+j-1) = isp(jutmp+j-1)
2940 5 jlmax = isp(ijl+n-1)
2942 jumax = isp(iju+n-1)
2943 l = (ju + jumax - 2 + lratio) / lratio + 1
2944 lmax = isp(il+n) - 1
2950 esp = umax - (isp(iu+n) - 1)
2952 if ((path-1) * (path-2) .ne. 0)
go to 6
2953 if (umax.lt.0)
go to 110
2955 * (n, r, c, ic, ia, ja, a, z, b,
2956 * lmax, isp(il), isp(jl), isp(ijl), rsp(l), rsp(d),
2957 * umax, isp(iu), isp(ju), isp(iju), rsp(u),
2958 * rsp(row), rsp(tmp), isp(irl), isp(jrl), flag)
2959 if(flag .ne. 0)
go to 100
2961 6
if ((path-3) .ne. 0)
go to 7
2963 * (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l),
2964 * rsp(d), isp(iu), isp(ju), isp(iju), rsp(u),
2967 7
if ((path-4) .ne. 0)
go to 8
2969 * (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l),
2970 * rsp(d), isp(iu), isp(ju), isp(iju), rsp(u),
2983 subroutine nroc (n, ic, ia, ja, a, jar, ar, p, flag)
3180 integer ic(*), ia(*), ja(*), jar(*), p(*), flag
3182 double precision a(*), ar(*)
3188 if(jmin .gt. jmax)
go to 5
3194 1
if(p(i) .ge. newj)
go to 2
3197 2
if(p(i) .eq. newj)
go to 102
3218 * (n, r, ic, ia,ja, jlmax,il,jl,ijl, jumax,iu,ju,iju,
3219 * q, ira,jra, irac, irl,jrl, iru,jru, flag)
3272 integer cend, qm, rend, rk, vj
3273 integer ia(*), ja(*), ira(*), jra(*), il(*), jl(*), ijl(*)
3274 integer iu(*), ju(*), iju(*), irl(*), jrl(*), iru(*), jru(*)
3275 integer r(*), ic(*), q(*), irac(*), flag
3294 if (iak .ge. ia(rk+1))
go to 101
3296 if (jaiak .gt. k)
go to 105
3297 jra(k) = irac(jaiak)
3309 if (vj .eq. 0)
go to 5
3313 if (qm .lt. vj)
go to 4
3314 if (qm .eq. vj)
go to 102
3319 if (vj .ne. 0)
go to 3
3326 if (i .eq. 0)
go to 10
3329 jmax = ijl(i) + il(i+1) - il(i) - 1
3331 if (long .lt. 0)
go to 6
3333 if (jtmp .ne. k) long = long + 1
3334 if (jtmp .eq. k) r(i) = -r(i)
3335 if (lastid .ge. long)
go to 7
3343 if (qm .lt. vj)
go to 8
3344 if (qm .eq. vj)
go to 9
3354 if (qm .ne. k)
go to 105
3355 if (luk .eq. 0)
go to 17
3356 if (lastid .ne. luk)
go to 11
3360 if (jl(irll) .ne. k) ijl(k) = ijl(k) - 1
3363 11
if (jlmin .gt. jlptr)
go to 15
3366 if (jl(j) - qm) 12, 13, 15
3371 if (jl(i) .ne. qm)
go to 15
3373 if (qm .gt. n)
go to 17
3377 15 jlmin = jlptr + 1
3379 if (luk .eq. 0)
go to 17
3381 if (jlptr .gt. jlmax)
go to 103
3387 il(k+1) = il(k) + luk
3396 if (jmin .gt. jmax)
go to 20
3402 if (qm .lt. vj)
go to 18
3403 if (qm .eq. vj)
go to 102
3415 if (i .eq. 0)
go to 26
3419 jmax = iju(i) + iu(i+1) - iu(i) - 1
3421 if (long .lt. 0)
go to 21
3423 if (jtmp .eq. k)
go to 22
3426 cend = ijl(i) + il(i+1) - il(i)
3428 if (irl(i) .ge. cend)
go to 22
3432 22
if (lastid .ge. long)
go to 23
3436 23
do 25 j=jmin,jmax
3440 if (qm .lt. vj)
go to 24
3441 if (qm .eq. vj)
go to 25
3449 26
if (il(k+1) .le. il(k))
go to 27
3456 if (qm .ne. k)
go to 105
3457 if (luk .eq. 0)
go to 34
3458 if (lastid .ne. luk)
go to 28
3462 if (ju(irul) .ne. k) iju(k) = iju(k) - 1
3465 28
if (jumin .gt. juptr)
go to 32
3468 if (ju(j) - qm) 29, 30, 32
3473 if (ju(i) .ne. qm)
go to 32
3475 if (qm .gt. n)
go to 34
3479 32 jumin = juptr + 1
3481 if (luk .eq. 0)
go to 34
3483 if (juptr .gt. jumax)
go to 106
3489 iu(k+1) = iu(k) + luk
3494 if (r(i) .lt. 0)
go to 36
3495 rend = iju(i) + iu(i+1) - iu(i)
3496 if (iru(i) .ge. rend)
go to 37
3503 if (i .eq. 0)
go to 38
3509 if (i .eq. 0)
go to 41
3512 if (ira(i) .ge. ia(r(i)+1))
go to 40
3514 jairai = ic(ja(irai))
3515 if (jairai .gt. i)
go to 40
3516 jra(i) = irac(jairai)
3519 if (i .ne. 0)
go to 39
3544 * (n, r,c,ic, ia,ja,a, z, b,
3545 * lmax,il,jl,ijl,l, d, umax,iu,ju,iju,u,
3546 * row, tmp, irl,jrl, flag)
3578 integer r(*), c(*), ic(*), ia(*), ja(*), il(*), jl(*), ijl(*)
3579 integer iu(*), ju(*), iju(*), irl(*), jrl(*), flag
3582 double precision a(*), l(*), d(*), u(*), z(*), b(*), row(*)
3583 double precision tmp(*), lki, sum, dk
3586 if(il(n+1)-1 .gt. lmax)
go to 104
3587 if(iu(n+1)-1 .gt. umax)
go to 107
3598 if (jrl(k) .eq. 0)
go to 3
3605 if (i .ne. 0)
go to 2
3608 jmax = jmin + iu(k+1) - iu(k) - 1
3609 if (jmin .gt. jmax)
go to 5
3617 row(ic(ja(j))) = a(j)
3622 if (i .eq. 0)
go to 10
3627 sum = sum + lki * tmp(i)
3630 if (jmin .gt. jmax)
go to 9
3633 8 row(ju(mu+j)) = row(ju(mu+j)) + lki * u(j)
3635 if (i .ne. 0)
go to 7
3638 10
if (row(k) .eq. 0.0d0)
go to 108
3642 if (k .eq. n)
go to 19
3645 if (jmin .gt. jmax)
go to 12
3648 11 u(j) = row(ju(mu+j)) * dk
3653 if (i .eq. 0)
go to 18
3654 14 irl(i) = irl(i) + 1
3656 if (irl(i) .ge. il(i+1))
go to 17
3657 ijlb = irl(i) - il(i) + ijl(i)
3659 15
if (i .gt. jrl(j))
go to 16
3665 if (i .ne. 0)
go to 14
3666 18
if (irl(k) .ge. il(k+1))
go to 19
3678 if (jmin .gt. jmax)
go to 21
3681 20 sum = sum - u(j) * tmp(ju(mu+j))
3699 * (n, r, c, il, jl, ijl, l, d, iu, ju, iju, u, z, b, tmp)
3716 integer r(*), c(*), il(*), jl(*), ijl(*), iu(*), ju(*), iju(*)
3718 double precision l(*), d(*), u(*), b(*), z(*), tmp(*), tmpk,sum
3727 tmpk = -d(k) * tmp(k)
3729 if (jmin .gt. jmax)
go to 3
3732 2 tmp(jl(ml+j)) = tmp(jl(ml+j)) + tmpk * l(j)
3740 if (jmin .gt. jmax)
go to 5
3743 4 sum = sum + u(j) * tmp(ju(mu+j))
3751 * (n, r, c, il, jl, ijl, l, d, iu, ju, iju, u, z, b, tmp)
3769 integer r(*), c(*), il(*), jl(*), ijl(*), iu(*), ju(*), iju(*)
3771 double precision l(*), d(*), u(*), b(*), z(*), tmp(*), tmpk,sum
3781 if (jmin .gt. jmax)
go to 3
3784 2 tmp(ju(mu+j)) = tmp(ju(mu+j)) + tmpk * u(j)
3792 if (jmin .gt. jmax)
go to 5
3795 4 sum = sum + l(j) * tmp(jl(ml+j))
3796 5 tmp(k) = -sum * d(k)
3803 SUBROUTINE dstoda (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
3804 1 WM, IWM, F, JAC, PJAC, SLVS)
3805 EXTERNAL f, jac, pjac, slvs
3806 INTEGER neq, nyh, iwm
3807 DOUBLE PRECISION y, yh, yh1, ewt, savf, acor, wm
3808 dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
3809 1 acor(*), wm(*), iwm(*)
3810 INTEGER iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
3811 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
3812 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
3813 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
3814 INTEGER iownd2, icount, irflag, jtyp, mused, mxordn, mxords
3815 DOUBLE PRECISION conit, crate, el, elco, hold, rmax, tesco,
3816 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
3817 DOUBLE PRECISION rownd2, cm1, cm2, pdest, pdlast, ratio,
3819 COMMON /dls001/ conit, crate, el(13), elco(13,12),
3820 1 hold, rmax, tesco(3,12),
3821 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
3822 3 iownd(6), ialth, ipup, lmax, meo, nqnyh, nslp,
3823 4 icf, ierpj, iersl, jcur, jstart, kflag, l,
3824 5 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
3825 6 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
3826 COMMON /dlsa01/ rownd2, cm1(12), cm2(5), pdest, pdlast, ratio,
3828 2 iownd2(3), icount, irflag, jtyp, mused, mxordn, mxords
3829 INTEGER i, i1, iredo, iret, j, jb, m, ncf, newq
3830 INTEGER lm1, lm1p1, lm2, lm2p1, nqm1, nqm2
3831 DOUBLE PRECISION dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup,
3832 1 r, rh, rhdn, rhsm, rhup, told, dmnorm
3833 DOUBLE PRECISION alpha, dm1,dm2, exm1,exm2,
3834 1 pdh, pnorm, rate, rh1, rh1it, rh2, rm, sm1(12)
3836 DATA sm1/0.5d0, 0.575d0, 0.55d0, 0.45d0, 0.35d0, 0.25d0,
3837 1 0.20d0, 0.15d0, 0.10d0, 0.075d0, 0.050d0, 0.025d0/
3926 IF (jstart .GT. 0)
GO TO 200
3927 IF (jstart .EQ. -1)
GO TO 100
3928 IF (jstart .EQ. -2)
GO TO 160
3956 CALL dcfode (2, elco, tesco)
3958 10 cm2(i) = tesco(2,i)*elco(i+1,i)
3959 CALL dcfode (1, elco, tesco)
3961 20 cm1(i) = tesco(2,i)*elco(i+1,i)
3976 IF (ialth .EQ. 1) ialth = 2
3977 IF (meth .EQ. mused)
GO TO 160
3978 CALL dcfode (meth, elco, tesco)
3986 155 el(i) = elco(i,nq)
3990 conit = 0.5d0/(nq+2)
3991 GO TO (160, 170, 200), iret
3998 160
IF (h .EQ. hold)
GO TO 200
4003 170 rh = max(rh,hmin/abs(h))
4004 175 rh = min(rh,rmax)
4005 rh = rh/max(1.0d0,abs(h)*hmxi*rh)
4011 IF (meth .EQ. 2)
GO TO 178
4013 pdh = max(abs(h)*pdlast,0.000001d0)
4014 IF (rh*pdh*1.00001d0 .LT. sm1(nq))
GO TO 178
4022 180 yh(i,j) = yh(i,j)*r
4026 IF (iredo .EQ. 0)
GO TO 690
4035 200
IF (abs(rc-1.0d0) .GT. ccmax) ipup = miter
4036 IF (nst .GE. nslp+msbp) ipup = miter
4043 210 yh1(i) = yh1(i) + yh1(i+nyh)
4045 pnorm = dmnorm(n, yh1, ewt)
4057 CALL f (neq, tn, y, savf)
4059 IF (ipup .LE. 0)
GO TO 250
4065 CALL pjac (neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac)
4070 IF (ierpj .NE. 0)
GO TO 430
4073 270
IF (miter .NE. 0)
GO TO 350
4079 savf(i) = h*savf(i) - yh(i,2)
4080 290 y(i) = savf(i) - acor(i)
4081 del = dmnorm(n, y, ewt)
4083 y(i) = yh(i,1) + el(1)*savf(i)
4084 300 acor(i) = savf(i)
4092 360 y(i) = h*savf(i) - (yh(i,2) + acor(i))
4093 CALL slvs (wm, iwm, y, savf)
4094 IF (iersl .LT. 0)
GO TO 430
4095 IF (iersl .GT. 0)
GO TO 410
4096 del = dmnorm(n, y, ewt)
4098 acor(i) = acor(i) + y(i)
4099 380 y(i) = yh(i,1) + el(1)*acor(i)
4113 IF (del .LE. 100.0d0*pnorm*uround)
GO TO 450
4114 IF (m .EQ. 0 .AND. meth .EQ. 1)
GO TO 405
4115 IF (m .EQ. 0)
GO TO 402
4117 IF (del .LE. 1024.0d0*delp) rm = del/delp
4119 crate = max(0.2d0*crate,rm)
4120 402 dcon = del*min(1.0d0,1.5d0*crate)/(tesco(2,nq)*conit)
4121 IF (dcon .GT. 1.0d0)
GO TO 405
4122 pdest = max(pdest,rate/abs(h*el(1)))
4123 IF (pdest .NE. 0.0d0) pdlast = pdest
4127 IF (m .EQ. maxcor)
GO TO 410
4128 IF (m .GE. 2 .AND. del .GT. 2.0d0*delp)
GO TO 410
4130 CALL f (neq, tn, y, savf)
4140 410
IF (miter .EQ. 0 .OR. jcur .EQ. 1)
GO TO 430
4153 440 yh1(i) = yh1(i) - yh1(i+nyh)
4155 IF (ierpj .LT. 0 .OR. iersl .LT. 0)
GO TO 680
4156 IF (abs(h) .LE. hmin*1.00001d0)
GO TO 670
4157 IF (ncf .EQ. mxncf)
GO TO 670
4169 IF (m .EQ. 0) dsm = del/tesco(2,nq)
4170 IF (m .GT. 0) dsm = dmnorm(n, acor, ewt)/tesco(2,nq)
4171 IF (dsm .GT. 1.0d0)
GO TO 500
4193 460 yh(i,j) = yh(i,j) + el(j)*acor(i)
4195 IF (icount .GE. 0)
GO TO 488
4196 IF (meth .EQ. 2)
GO TO 480
4215 IF (nq .GT. 5)
GO TO 488
4216 IF (dsm .GT. 100.0d0*pnorm*uround .AND. pdest .NE. 0.0d0)
4218 IF (irflag .EQ. 0)
GO TO 488
4220 nqm2 = min(nq,mxords)
4224 rh1 = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
4227 IF (pdh*rh1 .GT. 0.00001d0) rh1it = sm1(nq)/pdh
4228 rh1 = min(rh1,rh1it)
4229 IF (nq .LE. mxords)
GO TO 474
4234 dm2 = dmnorm(n, yh(1,lm2p1), ewt)/cm2(mxords)
4235 rh2 = 1.0d0/(1.2d0*dm2**exm2 + 0.0000012d0)
4237 474 dm2 = dsm*(cm1(nq)/cm2(nq))
4238 rh2 = 1.0d0/(1.2d0*dm2**exsm + 0.0000012d0)
4241 IF (rh2 .LT. ratio*rh1)
GO TO 488
4263 IF (mxordn .GE. nq)
GO TO 484
4268 dm1 = dmnorm(n, yh(1,lm1p1), ewt)/cm1(mxordn)
4269 rh1 = 1.0d0/(1.2d0*dm1**exm1 + 0.0000012d0)
4271 484 dm1 = dsm*(cm2(nq)/cm1(nq))
4272 rh1 = 1.0d0/(1.2d0*dm1**exsm + 0.0000012d0)
4275 486 rh1it = 2.0d0*rh1
4277 IF (pdh*rh1 .GT. 0.00001d0) rh1it = sm1(nqm1)/pdh
4278 rh1 = min(rh1,rh1it)
4279 rh2 = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
4280 IF (rh1*ratio .LT. 5.0d0*rh2)
GO TO 488
4281 alpha = max(0.001d0,rh1)
4282 dm1 = (alpha**exm1)*dm1
4283 IF (dm1 .LE. 1000.0d0*uround*pnorm)
GO TO 488
4297 IF (ialth .EQ. 0)
GO TO 520
4298 IF (ialth .GT. 1)
GO TO 700
4299 IF (l .EQ. lmax)
GO TO 700
4301 490 yh(i,lmax) = acor(i)
4310 500 kflag = kflag - 1
4317 510 yh1(i) = yh1(i) - yh1(i+nyh)
4320 IF (abs(h) .LE. hmin*1.00001d0)
GO TO 660
4321 IF (kflag .LE. -3)
GO TO 640
4335 IF (l .EQ. lmax)
GO TO 540
4337 530 savf(i) = acor(i) - yh(i,lmax)
4338 dup = dmnorm(n, savf, ewt)/tesco(3,nq)
4340 rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0)
4342 rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
4344 IF (nq .EQ. 1)
GO TO 550
4345 ddn = dmnorm(n, yh(1,l), ewt)/tesco(1,nq)
4347 rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
4349 550
IF (meth .EQ. 2)
GO TO 560
4350 pdh = max(abs(h)*pdlast,0.000001d0)
4351 IF (l .LT. lmax) rhup = min(rhup,sm1(l)/pdh)
4352 rhsm = min(rhsm,sm1(nq)/pdh)
4353 IF (nq .GT. 1) rhdn = min(rhdn,sm1(nq-1)/pdh)
4355 560
IF (rhsm .GE. rhup)
GO TO 570
4356 IF (rhup .GT. rhdn)
GO TO 590
4358 570
IF (rhsm .LT. rhdn)
GO TO 580
4364 IF (kflag .LT. 0 .AND. rh .GT. 1.0d0) rh = 1.0d0
4368 IF (rh .LT. 1.1d0)
GO TO 610
4371 600 yh(i,newq+1) = acor(i)*r
4376 620
IF (meth .EQ. 2)
GO TO 622
4377 IF (rh*pdh*1.00001d0 .GE. sm1(newq))
GO TO 625
4378 622
IF (kflag .EQ. 0 .AND. rh .LT. 1.1d0)
GO TO 610
4379 625
IF (kflag .LE. -2) rh = min(rh,0.2d0)
4385 IF (newq .EQ. nq)
GO TO 170
4399 640
IF (kflag .EQ. -10)
GO TO 660
4401 rh = max(hmin/abs(h),rh)
4405 CALL f (neq, tn, y, savf)
4408 650 yh(i,2) = h*savf(i)
4411 IF (nq .EQ. 1)
GO TO 200
4427 700 r = 1.0d0/tesco(2,nqu)
4429 710 acor(i) = acor(i)*r
4436 SUBROUTINE dprja (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
4439 INTEGER neq, nyh, iwm
4440 DOUBLE PRECISION y, yh, ewt, ftem, savf, wm
4441 dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
4443 INTEGER iownd, iowns,
4444 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
4445 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
4446 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
4447 INTEGER iownd2, iowns2, jtyp, mused, mxordn, mxords
4448 DOUBLE PRECISION rowns,
4449 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
4450 DOUBLE PRECISION rownd2, rowns2, pdnorm
4451 COMMON /dls001/ rowns(209),
4452 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
4453 2 iownd(6), iowns(6),
4454 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
4455 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
4456 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
4457 COMMON /dlsa01/ rownd2, rowns2(20), pdnorm,
4458 1 iownd2(3), iowns2(2), jtyp, mused, mxordn, mxords
4459 INTEGER i, i1, i2, ier, ii, j, j1, jj, lenp,
4460 1 mba, mband, meb1, meband, ml, ml3, mu, np1
4461 DOUBLE PRECISION con, fac, hl0, r, r0, srur, yi, yj, yjj,
4462 1 dmnorm, dfnorm, dbnorm
4501 GO TO (100, 200, 300, 400, 500), miter
4506 CALL jac (neq, tn, y, 0, 0, wm(3), n)
4509 120 wm(i+2) = wm(i+2)*con
4512 200 fac = dmnorm(n, savf, ewt)
4513 r0 = 1000.0d0*abs(h)*uround*n*fac
4514 IF (r0 .EQ. 0.0d0) r0 = 1.0d0
4519 r = max(srur*abs(yj),r0/ewt(j))
4522 CALL f (neq, tn, y, ftem)
4524 220 wm(i+j1) = (ftem(i) - savf(i))*fac
4531 pdnorm = dfnorm(n, wm(3), ewt)/abs(hl0)
4536 wm(j) = wm(j) + 1.0d0
4539 CALL dgefa (wm(3), n, n, iwm(21), ier)
4540 IF (ier .NE. 0) ierpj = 1
4553 CALL jac (neq, tn, y, ml, mu, wm(ml3), meband)
4556 420 wm(i+2) = wm(i+2)*con
4566 fac = dmnorm(n, savf, ewt)
4567 r0 = 1000.0d0*abs(h)*uround*n*fac
4568 IF (r0 .EQ. 0.0d0) r0 = 1.0d0
4570 DO 530 i = j,n,mband
4572 r = max(srur*abs(yi),r0/ewt(i))
4574 CALL f (neq, tn, y, ftem)
4575 DO 550 jj = j,n,mband
4578 r = max(srur*abs(yjj),r0/ewt(jj))
4582 ii = jj*meb1 - ml + 2
4584 540 wm(ii+i) = (ftem(i) - savf(i))*fac
4590 pdnorm = dbnorm(n, wm(ml+3), meband, ml, mu, ewt)/abs(hl0)
4594 wm(ii) = wm(ii) + 1.0d0
4595 580 ii = ii + meband
4597 CALL dgbfa (wm(3), meband, n, ml, mu, iwm(21), ier)
4598 IF (ier .NE. 0) ierpj = 1
4603 DOUBLE PRECISION FUNCTION dmnorm (N, V, W)
4611 DOUBLE PRECISION v, w, vm
4612 dimension v(n), w(n)
4615 10 vm = max(vm,abs(v(i))*w(i))
4621 DOUBLE PRECISION FUNCTION dfnorm (N, A, W)
4629 DOUBLE PRECISION a, w, an, sum
4630 dimension a(n,n), w(n)
4635 10 sum = sum + abs(a(i,j))/w(j)
4636 an = max(an,sum*w(i))
4643 DOUBLE PRECISION FUNCTION dbnorm (N, A, NRA, ML, MU, W)
4653 INTEGER n, nra, ml, mu
4654 INTEGER i, i1, jlo, jhi, j
4655 DOUBLE PRECISION a, w
4656 DOUBLE PRECISION an, sum
4657 dimension a(nra,n), w(n)
4665 10 sum = sum + abs(a(i1-j,j))/w(j)
4666 an = max(an,sum*w(i))
4673 SUBROUTINE dsrcma (RSAV, ISAV, JOB)
4688 INTEGER i, lenrls, lenils, lenrla, lenila
4689 DOUBLE PRECISION rsav
4690 DOUBLE PRECISION rls, rlsa
4691 dimension rsav(*), isav(*)
4692 SAVE lenrls, lenils, lenrla, lenila
4693 COMMON /dls001/ rls(218), ils(37)
4694 COMMON /dlsa01/ rlsa(22), ilsa(9)
4695 DATA lenrls/218/, lenils/37/, lenrla/22/, lenila/9/
4697 IF (job .EQ. 2)
GO TO 100
4701 15 rsav(lenrls+i) = rlsa(i)
4706 25 isav(lenils+i) = ilsa(i)
4712 110 rls(i) = rsav(i)
4714 115 rlsa(i) = rsav(lenrls+i)
4717 120 ils(i) = isav(i)
4719 125 ilsa(i) = isav(lenils+i)
4725 SUBROUTINE drchek (JOB, G, NEQ, Y, YH,NYH, G0, G1, GX, JROOT, IRT)
4727 INTEGER job, neq, nyh, jroot, irt
4728 DOUBLE PRECISION y, yh, g0, g1, gx
4729 dimension neq(*), y(*), yh(nyh,*), g0(*), g1(*), gx(*), jroot(*)
4730 INTEGER iownd, iowns,
4731 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
4732 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
4733 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
4734 INTEGER iownd3, iownr3, irfnd, itaskc, ngc, nge
4735 DOUBLE PRECISION rowns,
4736 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
4737 DOUBLE PRECISION rownr3, t0, tlast, toutc
4738 COMMON /dls001/ rowns(209),
4739 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
4740 2 iownd(6), iowns(6),
4741 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
4742 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
4743 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
4744 COMMON /dlsr01/ rownr3(2), t0, tlast, toutc,
4745 1 iownd3(3), iownr3(2), irfnd, itaskc, ngc, nge
4746 INTEGER i, iflag, jflag
4747 DOUBLE PRECISION hming, t1, temp1, temp2, x
4787 hming = (abs(tn) + abs(h))*uround*100.0d0
4789 GO TO (100, 200, 300), job
4794 CALL g (neq, t0, y, ngc, g0)
4798 110
IF (abs(g0(i)) .LE. 0.0d0) zroot = .true.
4799 IF (.NOT. zroot)
GO TO 190
4801 temp2 = max(hming/abs(h), 0.1d0)
4805 120 y(i) = y(i) + temp2*yh(i,2)
4806 CALL g (neq, t0, y, ngc, g0)
4810 130
IF (abs(g0(i)) .LE. 0.0d0) zroot = .true.
4811 IF (.NOT. zroot)
GO TO 190
4821 IF (irfnd .EQ. 0)
GO TO 260
4823 CALL dintdy (t0, 0, yh, nyh, y, iflag)
4824 CALL g (neq, t0, y, ngc, g0)
4828 210
IF (abs(g0(i)) .LE. 0.0d0) zroot = .true.
4829 IF (.NOT. zroot)
GO TO 260
4831 temp1 = sign(hming,h)
4833 IF ((t0 - tn)*h .LT. 0.0d0)
GO TO 230
4836 220 y(i) = y(i) + temp2*yh(i,2)
4838 230
CALL dintdy (t0, 0, yh, nyh, y, iflag)
4839 240
CALL g (neq, t0, y, ngc, g0)
4843 IF (abs(g0(i)) .GT. 0.0d0)
GO TO 250
4847 IF (.NOT. zroot)
GO TO 260
4852 260
IF (tn .EQ. tlast)
GO TO 390
4856 IF (itaskc.EQ.2 .OR. itaskc.EQ.3 .OR. itaskc.EQ.5)
GO TO 310
4857 IF ((toutc - tn)*h .GE. 0.0d0)
GO TO 310
4859 IF ((t1 - t0)*h .LE. 0.0d0)
GO TO 390
4860 CALL dintdy (t1, 0, yh, nyh, y, iflag)
4865 330
CALL g (neq, t1, y, ngc, g1)
4870 CALL droots (ngc, hming, jflag, t0, t1, g0, g1, gx, x, jroot)
4871 IF (jflag .GT. 1)
GO TO 360
4872 CALL dintdy (x, 0, yh, nyh, y, iflag)
4873 CALL g (neq, x, y, ngc, gx)
4877 CALL dcopy (ngc, gx, 1, g0, 1)
4878 IF (jflag .EQ. 4)
GO TO 390
4880 CALL dintdy (x, 0, yh, nyh, y, iflag)
4889 SUBROUTINE droots (NG, HMIN, JFLAG, X0, X1, G0, G1, GX, X, JROOT)
4890 INTEGER ng, jflag, jroot
4891 DOUBLE PRECISION hmin, x0, x1, g0, g1, gx, x
4892 dimension g0(ng), g1(ng), gx(ng), jroot(ng)
4893 INTEGER iownd3, imax, last, idum3
4894 DOUBLE PRECISION alpha, x2, rdum3
4895 COMMON /dlsr01/ alpha, x2, rdum3(3),
4896 1 iownd3(3), imax, last, idum3(4)
4970 INTEGER i, imxold, nxlast
4971 DOUBLE PRECISION t2, tmax, fracint, fracsub, zero,half,tenth,five
4972 LOGICAL zroot, sgnchg, xroot
4973 SAVE zero, half, tenth, five
4974 DATA zero/0.0d0/, half/0.5d0/, tenth/0.1d0/, five/5.0d0/
4976 IF (jflag .EQ. 1)
GO TO 200
4982 IF (abs(g1(i)) .GT. zero)
GO TO 110
4986 110
IF (sign(1.0d0,g0(i)) .EQ. sign(1.0d0,g1(i)))
GO TO 120
4987 t2 = abs(g1(i)/(g1(i)-g0(i)))
4988 IF (t2 .LE. tmax)
GO TO 120
4992 IF (imax .GT. 0)
GO TO 130
4996 140
IF (.NOT. sgnchg)
GO TO 400
5004 IF (xroot)
GO TO 300
5005 IF (nxlast .EQ. last)
GO TO 160
5008 160
IF (last .EQ. 0)
GO TO 170
5011 170 alpha = 2.0d0*alpha
5012 180 x2 = x1 - (x1 - x0)*g1(imax) / (g1(imax) - alpha*g0(imax))
5015 IF (abs(x2 - x0) < half*hmin)
THEN 5016 fracint = abs(x1 - x0)/hmin
5018 IF (fracint .LE. five) fracsub = half/fracint
5019 x2 = x0 + fracsub*(x1 - x0)
5021 IF (abs(x1 - x2) < half*hmin)
THEN 5022 fracint = abs(x1 - x0)/hmin
5024 IF (fracint .LE. five) fracsub = half/fracint
5025 x2 = x1 - fracsub*(x1 - x0)
5037 IF (abs(gx(i)) .GT. zero)
GO TO 210
5041 210
IF (sign(1.0d0,g0(i)) .EQ. sign(1.0d0,gx(i)))
GO TO 220
5042 t2 = abs(gx(i)/(gx(i) - g0(i)))
5043 IF (t2 .LE. tmax)
GO TO 220
5047 IF (imax .GT. 0)
GO TO 230
5053 IF (.NOT. sgnchg)
GO TO 250
5056 CALL dcopy (ng, gx, 1, g1, 1)
5060 250
IF (.NOT. zroot)
GO TO 260
5063 CALL dcopy (ng, gx, 1, g1, 1)
5068 CALL dcopy (ng, gx, 1, g0, 1)
5072 270
IF (abs(x1-x0) .LE. hmin) xroot = .true.
5078 CALL dcopy (ng, g1, 1, gx, 1)
5081 IF (abs(g1(i)) .GT. zero)
GO TO 310
5084 310
IF (sign(1.0d0,g0(i)) .NE. sign(1.0d0,g1(i))) jroot(i) = 1
5089 400
IF (.NOT. zroot)
GO TO 420
5093 CALL dcopy (ng, g1, 1, gx, 1)
5096 IF (abs(g1(i)) .LE. zero) jroot(i) = 1
5102 420
CALL dcopy (ng, g1, 1, gx, 1)
5109 SUBROUTINE dsrcar (RSAV, ISAV, JOB)
5123 INTEGER ils, ilsa, ilsr
5124 INTEGER i, ioff, lenrls, lenils, lenrla, lenila, lenrlr, lenilr
5125 DOUBLE PRECISION rsav
5126 DOUBLE PRECISION rls, rlsa, rlsr
5127 dimension rsav(*), isav(*)
5128 SAVE lenrls, lenils, lenrla, lenila, lenrlr, lenilr
5129 COMMON /dls001/ rls(218), ils(37)
5130 COMMON /dlsa01/ rlsa(22), ilsa(9)
5131 COMMON /dlsr01/ rlsr(5), ilsr(9)
5132 DATA lenrls/218/, lenils/37/, lenrla/22/, lenila/9/
5133 DATA lenrlr/5/, lenilr/9/
5135 IF (job .EQ. 2)
GO TO 100
5139 15 rsav(lenrls+i) = rlsa(i)
5140 ioff = lenrls + lenrla
5142 20 rsav(ioff+i) = rlsr(i)
5147 35 isav(lenils+i) = ilsa(i)
5148 ioff = lenils + lenila
5150 40 isav(ioff+i) = ilsr(i)
5156 110 rls(i) = rsav(i)
5158 115 rlsa(i) = rsav(lenrls+i)
5159 ioff = lenrls + lenrla
5161 120 rlsr(i) = rsav(ioff+i)
5164 130 ils(i) = isav(i)
5166 135 ilsa(i) = isav(lenils+i)
5167 ioff = lenils + lenila
5169 140 ilsr(i) = isav(ioff+i)
5175 SUBROUTINE dstodpk (NEQ, Y, YH, NYH, YH1, EWT, SAVF, SAVX, ACOR,
5176 1 WM, IWM, F, JAC, PSOL)
5177 EXTERNAL f, jac, psol
5178 INTEGER neq, nyh, iwm
5179 DOUBLE PRECISION y, yh, yh1, ewt, savf, savx, acor, wm
5180 dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
5181 1 savx(*), acor(*), wm(*), iwm(*)
5182 INTEGER iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
5183 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
5184 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
5185 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5186 INTEGER jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
5187 1 nni, nli, nps, ncfn, ncfl
5188 DOUBLE PRECISION conit, crate, el, elco, hold, rmax, tesco,
5189 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
5190 DOUBLE PRECISION delt, epcon, sqrtn, rsqrtn
5191 COMMON /dls001/ conit, crate, el(13), elco(13,12),
5192 1 hold, rmax, tesco(3,12),
5193 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
5194 3 iownd(6), ialth, ipup, lmax, meo, nqnyh, nslp,
5195 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
5196 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
5197 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5198 COMMON /dlpk01/ delt, epcon, sqrtn, rsqrtn,
5199 1 jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
5200 2 nni, nli, nps, ncfn, ncfl
5287 INTEGER i, i1, iredo, iret, j, jb, m, ncf, newq
5288 DOUBLE PRECISION dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup,
5289 1 r, rh, rhdn, rhsm, rhup, told, dvnorm
5299 IF (jstart .GT. 0)
GO TO 200
5300 IF (jstart .EQ. -1)
GO TO 100
5301 IF (jstart .EQ. -2)
GO TO 160
5339 IF (ialth .EQ. 1) ialth = 2
5340 IF (meth .EQ. meo)
GO TO 110
5341 CALL dcfode (meth, elco, tesco)
5343 IF (nq .GT. maxord)
GO TO 120
5347 110
IF (nq .LE. maxord)
GO TO 160
5351 125 el(i) = elco(i,nq)
5355 conit = 0.5d0/(nq+2)
5356 epcon = conit*tesco(2,nq)
5357 ddn = dvnorm(n, savf, ewt)/tesco(1,l)
5359 rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
5360 rh = min(rhdn,1.0d0)
5362 IF (h .EQ. hold)
GO TO 170
5363 rh = min(rh,abs(h/hold))
5371 140
CALL dcfode (meth, elco, tesco)
5373 155 el(i) = elco(i,nq)
5377 conit = 0.5d0/(nq+2)
5378 epcon = conit*tesco(2,nq)
5379 GO TO (160, 170, 200), iret
5386 160
IF (h .EQ. hold)
GO TO 200
5391 170 rh = max(rh,hmin/abs(h))
5392 175 rh = min(rh,rmax)
5393 rh = rh/max(1.0d0,abs(h)*hmxi*rh)
5398 180 yh(i,j) = yh(i,j)*r
5402 IF (iredo .EQ. 0)
GO TO 690
5412 200
IF (jacflg .NE. 0)
GO TO 202
5416 202
IF (abs(rc-1.0d0) .GT. ccmax) ipup = miter
5417 IF (nst .GE. nslp+msbp) ipup = miter
5424 210 yh1(i) = yh1(i) + yh1(i+nyh)
5436 CALL f (neq, tn, y, savf)
5438 IF (ipup .LE. 0)
GO TO 250
5444 CALL dpkset (neq, y, yh1, ewt, acor, savf, wm, iwm, f, jac)
5449 IF (ierpj .NE. 0)
GO TO 430
5452 270
IF (miter .NE. 0)
GO TO 350
5458 savf(i) = h*savf(i) - yh(i,2)
5459 290 y(i) = savf(i) - acor(i)
5460 del = dvnorm(n, y, ewt)
5462 y(i) = yh(i,1) + el(1)*savf(i)
5463 300 acor(i) = savf(i)
5471 360 savx(i) = h*savf(i) - (yh(i,2) + acor(i))
5472 CALL dsolpk (neq, y, savf, savx, ewt, wm, iwm, f, psol)
5473 IF (iersl .LT. 0)
GO TO 430
5474 IF (iersl .GT. 0)
GO TO 410
5475 del = dvnorm(n, savx, ewt)
5477 acor(i) = acor(i) + savx(i)
5478 380 y(i) = yh(i,1) + el(1)*acor(i)
5483 400
IF (m .NE. 0) crate = max(0.2d0*crate,del/delp)
5484 dcon = del*min(1.0d0,1.5d0*crate)/epcon
5485 IF (dcon .LE. 1.0d0)
GO TO 450
5487 IF (m .EQ. maxcor)
GO TO 410
5488 IF (m .GE. 2 .AND. del .GT. 2.0d0*delp)
GO TO 410
5491 CALL f (neq, tn, y, savf)
5501 410
IF (miter.EQ.0 .OR. jcur.EQ.1 .OR. jacflg.EQ.0)
GO TO 430
5515 440 yh1(i) = yh1(i) - yh1(i+nyh)
5517 IF (ierpj .LT. 0 .OR. iersl .LT. 0)
GO TO 680
5518 IF (abs(h) .LE. hmin*1.00001d0)
GO TO 670
5519 IF (ncf .EQ. mxncf)
GO TO 670
5531 IF (m .EQ. 0) dsm = del/tesco(2,nq)
5532 IF (m .GT. 0) dsm = dvnorm(n, acor, ewt)/tesco(2,nq)
5533 IF (dsm .GT. 1.0d0)
GO TO 500
5551 470 yh(i,j) = yh(i,j) + el(j)*acor(i)
5553 IF (ialth .EQ. 0)
GO TO 520
5554 IF (ialth .GT. 1)
GO TO 700
5555 IF (l .EQ. lmax)
GO TO 700
5557 490 yh(i,lmax) = acor(i)
5566 500 kflag = kflag - 1
5573 510 yh1(i) = yh1(i) - yh1(i+nyh)
5576 IF (abs(h) .LE. hmin*1.00001d0)
GO TO 660
5577 IF (kflag .LE. -3)
GO TO 640
5591 IF (l .EQ. lmax)
GO TO 540
5593 530 savf(i) = acor(i) - yh(i,lmax)
5594 dup = dvnorm(n, savf, ewt)/tesco(3,nq)
5596 rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0)
5598 rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
5600 IF (nq .EQ. 1)
GO TO 560
5601 ddn = dvnorm(n, yh(1,l), ewt)/tesco(1,nq)
5603 rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
5604 560
IF (rhsm .GE. rhup)
GO TO 570
5605 IF (rhup .GT. rhdn)
GO TO 590
5607 570
IF (rhsm .LT. rhdn)
GO TO 580
5613 IF (kflag .LT. 0 .AND. rh .GT. 1.0d0) rh = 1.0d0
5617 IF (rh .LT. 1.1d0)
GO TO 610
5620 600 yh(i,newq+1) = acor(i)*r
5624 620
IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1d0))
GO TO 610
5625 IF (kflag .LE. -2) rh = min(rh,0.2d0)
5631 IF (newq .EQ. nq)
GO TO 170
5645 640
IF (kflag .EQ. -10)
GO TO 660
5647 rh = max(hmin/abs(h),rh)
5651 CALL f (neq, tn, y, savf)
5654 650 yh(i,2) = h*savf(i)
5657 IF (nq .EQ. 1)
GO TO 200
5673 700 r = 1.0d0/tesco(2,nqu)
5675 710 acor(i) = acor(i)*r
5682 SUBROUTINE dpkset (NEQ, Y, YSV, EWT, FTEM, SAVF, WM, IWM, F, JAC)
5685 DOUBLE PRECISION y, ysv, ewt, ftem, savf, wm
5686 dimension neq(*), y(*), ysv(*), ewt(*), ftem(*), savf(*),
5688 INTEGER iownd, iowns,
5689 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
5690 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
5691 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5692 INTEGER jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
5693 1 nni, nli, nps, ncfn, ncfl
5694 DOUBLE PRECISION rowns,
5695 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
5696 DOUBLE PRECISION delt, epcon, sqrtn, rsqrtn
5697 COMMON /dls001/ rowns(209),
5698 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
5699 2 iownd(6), iowns(6),
5700 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
5701 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
5702 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5703 COMMON /dlpk01/ delt, epcon, sqrtn, rsqrtn,
5704 1 jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
5705 2 nni, nli, nps, ncfn, ncfl
5729 DOUBLE PRECISION hl0
5734 CALL jac (f, neq, tn, y, ysv, ewt, savf, ftem, hl0,
5735 1 wm(locwp), iwm(lociwp), ier)
5737 IF (ier .EQ. 0)
RETURN 5743 SUBROUTINE dsolpk (NEQ, Y, SAVF, X, EWT, WM, IWM, F, PSOL)
5746 DOUBLE PRECISION y, savf, x, ewt, wm
5747 dimension neq(*), y(*), savf(*), x(*), ewt(*), wm(*), iwm(*)
5748 INTEGER iownd, iowns,
5749 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
5750 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
5751 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5752 INTEGER jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
5753 1 nni, nli, nps, ncfn, ncfl
5754 DOUBLE PRECISION rowns,
5755 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
5756 DOUBLE PRECISION delt, epcon, sqrtn, rsqrtn
5757 COMMON /dls001/ rowns(209),
5758 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
5759 2 iownd(6), iowns(6),
5760 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
5761 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
5762 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5763 COMMON /dlpk01/ delt, epcon, sqrtn, rsqrtn,
5764 1 jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
5765 2 nni, nli, nps, ncfn, ncfl
5791 INTEGER iflag, lb, ldl, lhes, liom, lgmr, lpcg, lp, lq, lr,
5792 1 lv, lw, lwk, lz, maxlp1, npsl
5793 DOUBLE PRECISION delta, hl0
5798 GO TO (100, 200, 300, 400, 900, 900, 900, 900, 900), miter
5806 lwk = lhes + maxl*maxl
5807 CALL dcopy (n, x, 1, wm(lb), 1)
5808 CALL dscal (n, rsqrtn, ewt, 1)
5809 CALL dspiom (neq, tn, y, savf, wm(lb), ewt, n, maxl, kmp, delta,
5810 1 hl0, jpre, mnewt, f, psol, npsl, x, wm(lv), wm(lhes), iwm,
5811 2 liom, wm(locwp), iwm(lociwp), wm(lwk), iflag)
5815 CALL dscal (n, sqrtn, ewt, 1)
5816 IF (iflag .NE. 0) ncfl = ncfl + 1
5817 IF (iflag .GE. 2) iersl = 1
5818 IF (iflag .LT. 0) iersl = -1
5828 lq = lhes + maxl*maxlp1
5830 ldl = lwk + min(1,maxl-kmp)*n
5831 CALL dcopy (n, x, 1, wm(lb), 1)
5832 CALL dscal (n, rsqrtn, ewt, 1)
5833 CALL dspigmr (neq, tn, y, savf, wm(lb), ewt, n, maxl, maxlp1, kmp,
5834 1 delta, hl0, jpre, mnewt, f, psol, npsl, x, wm(lv), wm(lhes),
5835 2 wm(lq), lgmr, wm(locwp), iwm(lociwp), wm(lwk), wm(ldl), iflag)
5839 CALL dscal (n, sqrtn, ewt, 1)
5840 IF (iflag .NE. 0) ncfl = ncfl + 1
5841 IF (iflag .GE. 2) iersl = 1
5842 IF (iflag .LT. 0) iersl = -1
5853 CALL dcopy (n, x, 1, wm(lr), 1)
5854 CALL dpcg (neq, tn, y, savf, wm(lr), ewt, n, maxl, delta, hl0,
5855 1 jpre, mnewt, f, psol, npsl, x, wm(lp), wm(lw), wm(lz),
5856 2 lpcg, wm(locwp), iwm(lociwp), wm(lwk), iflag)
5860 IF (iflag .NE. 0) ncfl = ncfl + 1
5861 IF (iflag .GE. 2) iersl = 1
5862 IF (iflag .LT. 0) iersl = -1
5873 CALL dcopy (n, x, 1, wm(lr), 1)
5874 CALL dpcgs (neq, tn, y, savf, wm(lr), ewt, n, maxl, delta, hl0,
5875 1 jpre, mnewt, f, psol, npsl, x, wm(lp), wm(lw), wm(lz),
5876 2 lpcg, wm(locwp), iwm(lociwp), wm(lwk), iflag)
5880 IF (iflag .NE. 0) ncfl = ncfl + 1
5881 IF (iflag .GE. 2) iersl = 1
5882 IF (iflag .LT. 0) iersl = -1
5891 CALL dcopy (n, x, 1, wm(lb), 1)
5892 CALL dusol (neq, tn, y, savf, wm(lb), ewt, n, delta, hl0, mnewt,
5893 1 psol, npsl, x, wm(locwp), iwm(lociwp), wm(lwk), iflag)
5896 IF (iflag .NE. 0) ncfl = ncfl + 1
5897 IF (iflag .EQ. 3) iersl = 1
5898 IF (iflag .LT. 0) iersl = -1
5903 SUBROUTINE dspiom (NEQ, TN, Y, SAVF, B, WGHT, N, MAXL, KMP, DELTA,
5904 1 HL0, JPRE, MNEWT, F, PSOL, NPSL, X, V, HES, IPVT,
5905 2 LIOM, WP, IWP, WK, IFLAG)
5907 INTEGER neq,n,maxl,kmp,jpre,mnewt,npsl,ipvt,liom,iwp,iflag
5908 DOUBLE PRECISION tn,y,savf,b,wght,delta,hl0,x,v,hes,wp,wk
5909 dimension neq(*), y(*), savf(*), b(*), wght(*), x(*), v(n,*),
5910 1 hes(maxl,maxl), ipvt(*), wp(*), iwp(*), wk(*)
5990 INTEGER i, ier, info, j, k, ll, lm1
5991 DOUBLE PRECISION bnrm, bnrm0, prod, rho, snormw, dnrm2, tem
6001 10 v(i,1) = b(i)*wght(i)
6002 bnrm0 = dnrm2(n, v, 1)
6004 IF (bnrm0 .GT. delta)
GO TO 30
6005 IF (mnewt .GT. 0)
GO TO 20
6006 CALL dcopy (n, b, 1, x, 1)
6014 IF (jpre .EQ. 0 .OR. jpre .EQ. 2)
GO TO 55
6015 CALL psol (neq, tn, y, savf, wk, hl0, wp, iwp, b, 1, ier)
6017 IF (ier .NE. 0)
GO TO 300
6020 50 v(i,1) = b(i)*wght(i)
6021 bnrm = dnrm2(n, v, 1)
6022 delta = delta*(bnrm/bnrm0)
6024 CALL dscal (n, tem, v(1,1), 1)
6043 CALL datv (neq, y, savf, v(1,ll), wght, x, f, psol, v(1,ll+1),
6044 1 wk, wp, iwp, hl0, jpre, ier, npsl)
6045 IF (ier .NE. 0)
GO TO 300
6046 CALL dorthog (v(1,ll+1), v, hes, n, ll, maxl, kmp, snormw)
6047 CALL dhefa (hes, maxl, ll, ipvt, info, ll)
6049 IF (ll .GT. 1 .AND. ipvt(lm1) .EQ. lm1) prod = prod*hes(ll,lm1)
6050 IF (info .NE. ll)
GO TO 70
6056 IF (snormw .EQ. 0.0d0)
GO TO 120
6057 IF (ll .EQ. maxl)
GO TO 120
6065 rho = bnrm*snormw*abs(prod/hes(ll,ll))
6066 IF (rho .LE. delta)
GO TO 200
6067 IF (ll .EQ. maxl)
GO TO 100
6070 hes(ll+1,ll) = snormw
6072 CALL dscal (n, tem, v(1,ll+1), 1)
6080 IF (rho .LE. 1.0d0)
GO TO 150
6081 IF (rho .LE. bnrm .AND. mnewt .EQ. 0)
GO TO 150
6096 CALL dhesl (hes, maxl, ll, ipvt, b)
6100 CALL daxpy (n, b(i), v(1,i), 1, x, 1)
6103 240 x(i) = x(i)/wght(i)
6104 IF (jpre .LE. 1)
RETURN 6105 CALL psol (neq, tn, y, savf, wk, hl0, wp, iwp, x, 2, ier)
6107 IF (ier .NE. 0)
GO TO 300
6113 IF (ier .LT. 0) iflag = -1
6114 IF (ier .GT. 0) iflag = 3
6119 SUBROUTINE datv (NEQ, Y, SAVF, V, WGHT, FTEM, F, PSOL, Z, VTEM,
6120 1 WP, IWP, HL0, JPRE, IER, NPSL)
6122 INTEGER neq, iwp, jpre, ier, npsl
6123 DOUBLE PRECISION y, savf, v, wght, ftem, z, vtem, wp, hl0
6124 dimension neq(*), y(*), savf(*), v(*), wght(*), ftem(*), z(*),
6125 1 vtem(*), wp(*), iwp(*)
6126 INTEGER iownd, iowns,
6127 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
6128 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
6129 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
6130 DOUBLE PRECISION rowns,
6131 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
6132 COMMON /dls001/ rowns(209),
6133 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
6134 2 iownd(6), iowns(6),
6135 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
6136 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
6137 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
6189 DOUBLE PRECISION fac, rnorm, dnrm2, tempn
6193 10 vtem(i) = v(i)/wght(i)
6195 IF (jpre .GE. 2)
GO TO 30
6198 CALL dcopy (n, y, 1, z, 1)
6200 20 y(i) = z(i) + vtem(i)
6206 CALL psol (neq, tn, y, savf, ftem, hl0, wp, iwp, vtem, 2, ier)
6208 IF (ier .NE. 0)
RETURN 6211 40 z(i) = vtem(i)*wght(i)
6212 tempn = dnrm2(n, z, 1)
6215 CALL dcopy (n, y, 1, z, 1)
6217 50 y(i) = z(i) + vtem(i)*rnorm
6222 CALL f (neq, tn, y, ftem)
6224 CALL dcopy (n, z, 1, y, 1)
6227 70 z(i) = ftem(i) - savf(i)
6229 80 z(i) = vtem(i) - fac*z(i)
6231 IF (jpre .EQ. 0 .OR. jpre .EQ. 2)
GO TO 85
6232 CALL psol (neq, tn, y, savf, ftem, hl0, wp, iwp, z, 1, ier)
6234 IF (ier .NE. 0)
RETURN 6238 90 z(i) = z(i)*wght(i)
6243 SUBROUTINE dorthog (VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
6244 INTEGER n, ll, ldhes, kmp
6245 DOUBLE PRECISION vnew, v, hes, snormw
6246 dimension vnew(*), v(n,*), hes(ldhes,*)
6288 DOUBLE PRECISION arg, ddot, dnrm2, sumdsq, tem, vnrm
6291 vnrm = dnrm2(n, vnew, 1)
6297 i0 = max(1,ll-kmp+1)
6299 hes(i,ll) = ddot(n, v(1,i), 1, vnew, 1)
6301 CALL daxpy (n, tem, v(1,i), 1, vnew, 1)
6310 snormw = dnrm2(n, vnew, 1)
6311 IF (vnrm + 0.001d0*snormw .NE. vnrm)
RETURN 6314 tem = -ddot(n, v(1,i), 1, vnew, 1)
6315 IF (hes(i,ll) + 0.001d0*tem .EQ. hes(i,ll))
GO TO 30
6316 hes(i,ll) = hes(i,ll) - tem
6317 CALL daxpy (n, tem, v(1,i), 1, vnew, 1)
6318 sumdsq = sumdsq + tem**2
6320 IF (sumdsq .EQ. 0.0d0)
RETURN 6321 arg = max(0.0d0,snormw**2 - sumdsq)
6328 SUBROUTINE dspigmr (NEQ, TN, Y, SAVF, B, WGHT, N, MAXL, MAXLP1,
6329 1 KMP, DELTA, HL0, JPRE, MNEWT, F, PSOL, NPSL, X, V, HES, Q,
6330 2 LGMR, WP, IWP, WK, DL, IFLAG)
6332 INTEGER neq,n,maxl,maxlp1,kmp,jpre,mnewt,npsl,lgmr,iwp,iflag
6333 DOUBLE PRECISION tn,y,savf,b,wght,delta,hl0,x,v,hes,q,wp,wk,dl
6334 dimension neq(*), y(*), savf(*), b(*), wght(*), x(*), v(n,*),
6335 1 hes(maxlp1,*), q(*), wp(*), iwp(*), wk(*), dl(*)
6423 INTEGER i, ier, info, ip1, i2, j, k, ll, llp1
6424 DOUBLE PRECISION bnrm,bnrm0,c,dlnrm,prod,rho,s,snormw,dnrm2,tem
6434 10 v(i,1) = b(i)*wght(i)
6435 bnrm0 = dnrm2(n, v, 1)
6437 IF (bnrm0 .GT. delta)
GO TO 30
6438 IF (mnewt .GT. 0)
GO TO 20
6439 CALL dcopy (n, b, 1, x, 1)
6447 IF (jpre .EQ. 0 .OR. jpre .EQ. 2)
GO TO 55
6448 CALL psol (neq, tn, y, savf, wk, hl0, wp, iwp, b, 1, ier)
6450 IF (ier .NE. 0)
GO TO 300
6453 50 v(i,1) = b(i)*wght(i)
6454 bnrm = dnrm2(n, v, 1)
6455 delta = delta*(bnrm/bnrm0)
6457 CALL dscal (n, tem, v(1,1), 1)
6476 CALL datv (neq, y, savf, v(1,ll), wght, x, f, psol, v(1,ll+1),
6477 1 wk, wp, iwp, hl0, jpre, ier, npsl)
6478 IF (ier .NE. 0)
GO TO 300
6479 CALL dorthog (v(1,ll+1), v, hes, n, ll, maxlp1, kmp, snormw)
6480 hes(ll+1,ll) = snormw
6481 CALL dheqr (hes, maxlp1, ll, q, info, ll)
6482 IF (info .EQ. ll)
GO TO 120
6490 rho = abs(prod*bnrm)
6491 IF ((ll.GT.kmp) .AND. (kmp.LT.maxl))
THEN 6492 IF (ll .EQ. kmp+1)
THEN 6493 CALL dcopy (n, v(1,1), 1, dl, 1)
6500 70 dl(k) = s*dl(k) + c*v(k,ip1)
6504 c = q(2*ll-1)/snormw
6507 80 dl(k) = s*dl(k) + c*v(k,llp1)
6508 dlnrm = dnrm2(n, dl, 1)
6515 IF (rho .LE. delta)
GO TO 200
6516 IF (ll .EQ. maxl)
GO TO 100
6521 CALL dscal (n, tem, v(1,ll+1), 1)
6524 IF (rho .LE. 1.0d0)
GO TO 150
6525 IF (rho .LE. bnrm .AND. mnewt .EQ. 0)
GO TO 150
6541 CALL dhels (hes, maxlp1, ll, q, b)
6545 CALL daxpy (n, b(i), v(1,i), 1, x, 1)
6548 240 x(i) = x(i)/wght(i)
6549 IF (jpre .LE. 1)
RETURN 6550 CALL psol (neq, tn, y, savf, wk, hl0, wp, iwp, x, 2, ier)
6552 IF (ier .NE. 0)
GO TO 300
6558 IF (ier .LT. 0) iflag = -1
6559 IF (ier .GT. 0) iflag = 3
6565 SUBROUTINE dpcg (NEQ, TN, Y, SAVF, R, WGHT, N, MAXL, DELTA, HL0,
6566 1 JPRE, MNEWT, F, PSOL, NPSL, X, P, W, Z, LPCG, WP, IWP, WK, IFLAG)
6568 INTEGER neq, n, maxl, jpre, mnewt, npsl, lpcg, iwp, iflag
6569 DOUBLE PRECISION tn,y,savf,r,wght,delta,hl0,x,p,w,z,wp,wk
6570 dimension neq(*), y(*), savf(*), r(*), wght(*), x(*), p(*), w(*),
6571 1 z(*), wp(*), iwp(*), wk(*)
6640 DOUBLE PRECISION alpha,beta,bnrm,ptw,rnrm,ddot,dvnorm,ztr,ztr0
6647 bnrm = dvnorm(n, r, wght)
6649 IF (bnrm .GT. delta)
GO TO 20
6650 IF (mnewt .GT. 0)
RETURN 6651 CALL dcopy (n, r, 1, x, 1)
6658 CALL dcopy (n, r, 1, z, 1)
6660 IF (jpre .EQ. 0)
GO TO 40
6661 CALL psol (neq, tn, y, savf, wk, hl0, wp, iwp, z, 3, ier)
6663 IF (ier .NE. 0)
GO TO 100
6666 ztr = ddot(n, z, 1, r, 1)
6667 IF (lpcg .NE. 1)
GO TO 50
6668 CALL dcopy (n, z, 1, p, 1)
6671 IF (ztr0 .EQ. 0.0d0)
GO TO 200
6674 60 p(i) = z(i) + beta*p(i)
6679 CALL datp (neq, y, savf, p, wght, hl0, wk, f, w)
6681 ptw = ddot(n, p, 1, w, 1)
6682 IF (ptw .EQ. 0.0d0)
GO TO 200
6684 CALL daxpy (n, alpha, p, 1, x, 1)
6686 CALL daxpy (n, alpha, w, 1, r, 1)
6687 rnrm = dvnorm(n, r, wght)
6688 IF (rnrm .LE. delta)
RETURN 6689 IF (lpcg .LT. maxl)
GO TO 30
6691 IF (rnrm .LE. 1.0d0) iflag = 1
6692 IF (rnrm .LE. bnrm .AND. mnewt .EQ. 0) iflag = 1
6698 IF (ier .LT. 0) iflag = -1
6699 IF (ier .GT. 0) iflag = 3
6710 SUBROUTINE dpcgs (NEQ, TN, Y, SAVF, R, WGHT, N, MAXL, DELTA, HL0,
6711 1 JPRE, MNEWT, F, PSOL, NPSL, X, P, W, Z, LPCG, WP, IWP, WK, IFLAG)
6713 INTEGER neq, n, maxl, jpre, mnewt, npsl, lpcg, iwp, iflag
6714 DOUBLE PRECISION tn,y,savf,r,wght,delta,hl0,x,p,w,z,wp,wk
6715 dimension neq(*), y(*), savf(*), r(*), wght(*), x(*), p(*), w(*),
6716 1 z(*), wp(*), iwp(*), wk(*)
6786 DOUBLE PRECISION alpha, beta, bnrm, ptw, rnrm, dvnorm, ztr, ztr0
6793 bnrm = dvnorm(n, r, wght)
6795 IF (bnrm .GT. delta)
GO TO 20
6796 IF (mnewt .GT. 0)
RETURN 6797 CALL dcopy (n, r, 1, x, 1)
6804 CALL dcopy (n, r, 1, z, 1)
6806 IF (jpre .EQ. 0)
GO TO 40
6807 CALL psol (neq, tn, y, savf, wk, hl0, wp, iwp, z, 3, ier)
6809 IF (ier .NE. 0)
GO TO 100
6814 45 ztr = ztr + z(i)*r(i)*wght(i)**2
6815 IF (lpcg .NE. 1)
GO TO 50
6816 CALL dcopy (n, z, 1, p, 1)
6819 IF (ztr0 .EQ. 0.0d0)
GO TO 200
6822 60 p(i) = z(i) + beta*p(i)
6827 CALL datp (neq, y, savf, p, wght, hl0, wk, f, w)
6831 80 ptw = ptw + p(i)*w(i)*wght(i)**2
6832 IF (ptw .EQ. 0.0d0)
GO TO 200
6834 CALL daxpy (n, alpha, p, 1, x, 1)
6836 CALL daxpy (n, alpha, w, 1, r, 1)
6837 rnrm = dvnorm(n, r, wght)
6838 IF (rnrm .LE. delta)
RETURN 6839 IF (lpcg .LT. maxl)
GO TO 30
6841 IF (rnrm .LE. 1.0d0) iflag = 1
6842 IF (rnrm .LE. bnrm .AND. mnewt .EQ. 0) iflag = 1
6848 IF (ier .LT. 0) iflag = -1
6849 IF (ier .GT. 0) iflag = 3
6860 SUBROUTINE datp (NEQ, Y, SAVF, P, WGHT, HL0, WK, F, W)
6863 DOUBLE PRECISION y, savf, p, wght, hl0, wk, w
6864 dimension neq(*), y(*), savf(*), p(*), wght(*), wk(*), w(*)
6865 INTEGER iownd, iowns,
6866 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
6867 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
6868 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
6869 DOUBLE PRECISION rowns,
6870 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
6871 COMMON /dls001/ rowns(209),
6872 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
6873 2 iownd(6), iowns(6),
6874 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
6875 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
6876 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
6909 DOUBLE PRECISION fac, pnrm, rpnrm, dvnorm
6911 pnrm = dvnorm(n, p, wght)
6913 CALL dcopy (n, y, 1, w, 1)
6915 20 y(i) = w(i) + p(i)*rpnrm
6916 CALL f (neq, tn, y, wk)
6918 CALL dcopy (n, w, 1, y, 1)
6921 40 w(i) = p(i) - fac*(wk(i) - savf(i))
6926 SUBROUTINE dusol (NEQ, TN, Y, SAVF, B, WGHT, N, DELTA, HL0, MNEWT,
6927 1 PSOL, NPSL, X, WP, IWP, WK, IFLAG)
6929 INTEGER neq, n, mnewt, npsl, iwp, iflag
6930 DOUBLE PRECISION tn, y, savf, b, wght, delta, hl0, x, wp, wk
6931 dimension neq(*), y(*), savf(*), b(*), wght(*), x(*),
6932 1 wp(*), iwp(*), wk(*)
6986 DOUBLE PRECISION bnrm, dvnorm
6993 bnrm = dvnorm(n, b, wght)
6994 IF (bnrm .GT. delta)
GO TO 30
6995 IF (mnewt .GT. 0)
GO TO 10
6996 CALL dcopy (n, b, 1, x, 1)
7003 CALL psol (neq, tn, y, savf, wk, hl0, wp, iwp, b, 0, ier)
7005 IF (ier .NE. 0)
GO TO 100
7006 CALL dcopy (n, b, 1, x, 1)
7012 IF (ier .LT. 0) iflag = -1
7013 IF (ier .GT. 0) iflag = 3
7018 SUBROUTINE dsrcpk (RSAV, ISAV, JOB)
7033 INTEGER i, lenilp, lenrlp, lenils, lenrls
7034 DOUBLE PRECISION rsav, rls, rlsp
7035 dimension rsav(*), isav(*)
7036 SAVE lenrls, lenils, lenrlp, lenilp
7037 COMMON /dls001/ rls(218), ils(37)
7038 COMMON /dlpk01/ rlsp(4), ilsp(13)
7039 DATA lenrls/218/, lenils/37/, lenrlp/4/, lenilp/13/
7041 IF (job .EQ. 2)
GO TO 100
7042 CALL dcopy (lenrls, rls, 1, rsav, 1)
7043 CALL dcopy (lenrlp, rlsp, 1, rsav(lenrls+1), 1)
7047 40 isav(lenils+i) = ilsp(i)
7051 CALL dcopy (lenrls, rsav, 1, rls, 1)
7052 CALL dcopy (lenrlp, rsav(lenrls+1), 1, rlsp, 1)
7054 120 ils(i) = isav(i)
7056 140 ilsp(i) = isav(lenils+i)
7061 SUBROUTINE dhefa (A, LDA, N, IPVT, INFO, JOB)
7062 INTEGER lda, n, ipvt(*), info, job
7063 DOUBLE PRECISION a(lda,*)
7115 INTEGER idamax, j, k, km1, kp1, l, nm1
7118 IF (job .GT. 1)
GO TO 80
7128 IF (nm1 .LT. 1)
GO TO 70
7134 l = idamax(2, a(k,k), 1) + k - 1
7139 IF (a(l,k) .EQ. 0.0d0)
GO TO 40
7143 IF (l .EQ. k)
GO TO 10
7152 a(k+1,k) = a(k+1,k)*t
7158 IF (l .EQ. k)
GO TO 20
7162 CALL daxpy (n-k, t, a(k+1,k), 1, a(k+1,j), 1)
7171 IF (a(n,n) .EQ. 0.0d0) info = n
7184 IF (nm1 .LE. 1)
GO TO 105
7189 IF (l .EQ. km1)
GO TO 90
7193 a(k,n) = a(k,n) + a(k,km1)*t
7203 l = idamax(2, a(nm1,nm1), 1) + nm1 - 1
7208 IF (a(l,nm1) .EQ. 0.0d0)
GO TO 140
7212 IF (l .EQ. nm1)
GO TO 110
7214 a(l,nm1) = a(nm1,nm1)
7220 t = -1.0d0/a(nm1,nm1)
7221 a(n,nm1) = a(n,nm1)*t
7226 IF (l .EQ. nm1)
GO TO 120
7230 a(n,n) = a(n,n) + t*a(n,nm1)
7236 IF (a(n,n) .EQ. 0.0d0) info = n
7241 SUBROUTINE dhesl (A, LDA, N, IPVT, B)
7242 INTEGER lda, n, ipvt(*)
7243 DOUBLE PRECISION a(lda,*), b(*)
7277 INTEGER k, kb, l, nm1
7285 IF (nm1 .LT. 1)
GO TO 30
7289 IF (l .EQ. k)
GO TO 10
7293 b(k+1) = b(k+1) + t*a(k+1,k)
7303 CALL daxpy (k-1, t, a(1,k), 1, b(1), 1)
7309 SUBROUTINE dheqr (A, LDA, N, Q, INFO, IJOB)
7310 INTEGER lda, n, info, ijob
7311 DOUBLE PRECISION a(lda,*), q(*)
7361 INTEGER i, iq, j, k, km1, kp1, nm1
7362 DOUBLE PRECISION c, s, t, t1, t2
7364 IF (ijob .GT. 1)
GO TO 70
7379 IF (km1 .LT. 1)
GO TO 20
7386 a(j,k) = c*t1 - s*t2
7387 a(j+1,k) = s*t1 + c*t2
7396 IF (t2 .NE. 0.0d0)
GO TO 30
7401 IF (abs(t2) .LT. abs(t1))
GO TO 40
7403 s = -1.0d0/sqrt(1.0d0+t*t)
7408 c = 1.0d0/sqrt(1.0d0+t*t)
7413 a(k,k) = c*t1 - s*t2
7414 IF (a(k,k) .EQ. 0.0d0) info = k
7433 a(k,n) = c*t1 - s*t2
7434 a(k+1,n) = s*t1 + c*t2
7443 IF (t2 .NE. 0.0d0)
GO TO 110
7448 IF (abs(t2) .LT. abs(t1))
GO TO 120
7450 s = -1.0d0/sqrt(1.0d0+t*t)
7455 c = 1.0d0/sqrt(1.0d0+t*t)
7461 a(n,n) = c*t1 - s*t2
7462 IF (a(n,n) .EQ. 0.0d0) info = n
7467 SUBROUTINE dhels (A, LDA, N, Q, B)
7469 DOUBLE PRECISION a(lda,*), b(*), q(*)
7508 INTEGER iq, k, kb, kp1
7509 DOUBLE PRECISION c, s, t, t1, t2
7522 b(kp1) = s*t1 + c*t2
7531 CALL daxpy (k-1, t, a(1,k), 1, b(1), 1)
7537 SUBROUTINE dlhin (NEQ, N, T0, Y0, YDOT, F, TOUT, UROUND,
7538 1 EWT, ITOL, ATOL, Y, TEMP, H0, NITER, IER)
7540 DOUBLE PRECISION t0, y0, ydot, tout, uround, ewt, atol, y,
7542 INTEGER neq, n, itol, niter, ier
7543 dimension neq(*), y0(*), ydot(*), ewt(*), atol(*), y(*), temp(*)
7585 DOUBLE PRECISION afi, atoli, delyi, half, hg, hlb, hnew, hrat,
7586 1 hub, hun, pt1, t1, tdist, tround, two, dvnorm, yddnrm
7592 SAVE half, hun, pt1, two
7593 DATA half /0.5d0/, hun /100.0d0/, pt1 /0.1d0/, two /2.0d0/
7596 tdist = abs(tout - t0)
7597 tround = uround*max(abs(t0),abs(tout))
7598 IF (tdist .LT. two*tround)
GO TO 100
7606 IF (itol .EQ. 2 .OR. itol .EQ. 4) atoli = atol(i)
7607 delyi = pt1*abs(y0(i)) + atoli
7609 IF (afi*hub .GT. delyi) hub = delyi/afi
7616 IF (hub .LT. hlb)
THEN 7626 60 y(i) = y0(i) + hg*ydot(i)
7627 CALL f (neq, t1, y, temp)
7629 70 temp(i) = (temp(i) - ydot(i))/hg
7630 yddnrm = dvnorm(n, temp, ewt)
7632 IF (yddnrm*hub*hub .GT. two)
THEN 7633 hnew = sqrt(two/yddnrm)
7645 IF (iter .GE. 4)
GO TO 80
7647 IF ( (hrat .GT. half) .AND. (hrat .LT. two) )
GO TO 80
7648 IF ( (iter .GE. 2) .AND. (hnew .GT. two*hg) )
THEN 7657 IF (h0 .LT. hlb) h0 = hlb
7658 IF (h0 .GT. hub) h0 = hub
7659 90 h0 = sign(h0, tout - t0)
7661 CALL dcopy (n, y0, 1, y, 1)
7671 SUBROUTINE dstoka (NEQ, Y, YH, NYH, YH1, EWT, SAVF, SAVX, ACOR,
7672 1 WM, IWM, F, JAC, PSOL)
7673 EXTERNAL f, jac, psol
7674 INTEGER neq, nyh, iwm
7675 DOUBLE PRECISION y, yh, yh1, ewt, savf, savx, acor, wm
7676 dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
7677 1 savx(*), acor(*), wm(*), iwm(*)
7678 INTEGER iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
7679 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
7680 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
7681 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
7682 INTEGER newt, nsfi, nslj, njev
7683 INTEGER jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
7684 1 nni, nli, nps, ncfn, ncfl
7685 DOUBLE PRECISION conit, crate, el, elco, hold, rmax, tesco,
7686 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
7687 DOUBLE PRECISION stifr
7688 DOUBLE PRECISION delt, epcon, sqrtn, rsqrtn
7689 COMMON /dls001/ conit, crate, el(13), elco(13,12),
7690 1 hold, rmax, tesco(3,12),
7691 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
7692 3 iownd(6), ialth, ipup, lmax, meo, nqnyh, nslp,
7693 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
7694 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
7695 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
7696 COMMON /dls002/ stifr, newt, nsfi, nslj, njev
7697 COMMON /dlpk01/ delt, epcon, sqrtn, rsqrtn,
7698 1 jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
7699 2 nni, nli, nps, ncfn, ncfl
7774 INTEGER i, i1, iredo, iret, j, jb, jok, m, ncf, newq, nslow
7775 DOUBLE PRECISION dcon, ddn, del, delp, drc, dsm, dup, exdn, exsm,
7776 1 exup, dfnorm, r, rh, rhdn, rhsm, rhup, roc, stiff, told, dvnorm
7786 IF (jstart .GT. 0)
GO TO 200
7787 IF (jstart .EQ. -1)
GO TO 100
7788 IF (jstart .EQ. -2)
GO TO 160
7829 IF (ialth .EQ. 1) ialth = 2
7830 IF (meth .EQ. meo)
GO TO 110
7831 CALL dcfode (meth, elco, tesco)
7833 IF (nq .GT. maxord)
GO TO 120
7837 110
IF (nq .LE. maxord)
GO TO 160
7841 125 el(i) = elco(i,nq)
7845 conit = 0.5d0/(nq+2)
7846 epcon = conit*tesco(2,nq)
7847 ddn = dvnorm(n, savf, ewt)/tesco(1,l)
7849 rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
7850 rh = min(rhdn,1.0d0)
7852 IF (h .EQ. hold)
GO TO 170
7853 rh = min(rh,abs(h/hold))
7861 140
CALL dcfode (meth, elco, tesco)
7863 155 el(i) = elco(i,nq)
7867 conit = 0.5d0/(nq+2)
7868 epcon = conit*tesco(2,nq)
7869 GO TO (160, 170, 200), iret
7876 160
IF (h .EQ. hold)
GO TO 200
7881 170 rh = max(rh,hmin/abs(h))
7882 175 rh = min(rh,rmax)
7883 rh = rh/max(1.0d0,abs(h)*hmxi*rh)
7888 180 yh(i,j) = yh(i,j)*r
7892 IF (iredo .EQ. 0)
GO TO 690
7902 200
IF (newt .EQ. 0 .OR. jacflg .EQ. 0)
THEN 7907 drc = abs(rc - 1.0d0)
7908 IF (drc .GT. ccmax) ipup = miter
7909 IF (nst .GE. nslp+msbp) ipup = miter
7917 210 yh1(i) = yh1(i) + yh1(i+nyh)
7935 CALL f (neq, tn, y, savf)
7937 IF (newt .EQ. 0 .OR. ipup .LE. 0)
GO TO 250
7945 IF (nst .EQ. 0 .OR. nst .GT. nslj+50) jok = -1
7946 IF (icf .EQ. 1 .AND. drc .LT. 0.2d0) jok = -1
7947 IF (icf .EQ. 2) jok = -1
7948 IF (jok .EQ. -1)
THEN 7952 CALL dsetpk (neq, y, yh1, ewt, acor, savf, jok, wm, iwm, f, jac)
7958 IF (ierpj .NE. 0)
GO TO 430
7961 270
IF (newt .NE. 0)
GO TO 350
7967 savf(i) = h*savf(i) - yh(i,2)
7968 290 y(i) = savf(i) - acor(i)
7969 del = dvnorm(n, y, ewt)
7971 y(i) = yh(i,1) + el(1)*savf(i)
7972 300 acor(i) = savf(i)
7982 360 savx(i) = h*savf(i) - (yh(i,2) + acor(i))
7983 dfnorm = dvnorm(n, savx, ewt)
7984 CALL dsolpk (neq, y, savf, savx, ewt, wm, iwm, f, psol)
7985 IF (iersl .LT. 0)
GO TO 430
7986 IF (iersl .GT. 0)
GO TO 410
7987 del = dvnorm(n, savx, ewt)
7988 IF (del .GT. 1.0d-8) stiff = max(stiff, dfnorm/del)
7990 acor(i) = acor(i) + savx(i)
7991 380 y(i) = yh(i,1) + el(1)*acor(i)
7998 400
IF (m .NE. 0)
THEN 7999 roc = max(0.05d0, del/delp)
8000 crate = max(0.2d0*crate,roc)
8002 dcon = del*min(1.0d0,1.5d0*crate)/epcon
8003 IF (dcon .LE. 1.0d0)
GO TO 450
8005 IF (m .EQ. maxcor)
GO TO 410
8006 IF (m .GE. 2 .AND. del .GT. 2.0d0*delp)
GO TO 410
8007 IF (roc .GT. 10.0d0)
GO TO 410
8008 IF (roc .GT. 0.8d0) nslow = nslow + 1
8009 IF (nslow .GE. 2)
GO TO 410
8012 CALL f (neq, tn, y, savf)
8029 IF (newt .EQ. 0)
THEN 8030 IF (nst .EQ. 0)
GO TO 430
8031 IF (miter .EQ. 0)
GO TO 430
8037 IF (jcur.EQ.1 .OR. jacflg.EQ.0)
GO TO 430
8050 440 yh1(i) = yh1(i) - yh1(i+nyh)
8052 IF (ierpj .LT. 0 .OR. iersl .LT. 0)
GO TO 680
8053 IF (abs(h) .LE. hmin*1.00001d0)
GO TO 670
8054 IF (ncf .EQ. mxncf)
GO TO 670
8067 IF (newt .GT. 0) stifr = 0.5d0*(stifr + stiff)
8068 IF (m .EQ. 0) dsm = del/tesco(2,nq)
8069 IF (m .GT. 0) dsm = dvnorm(n, acor, ewt)/tesco(2,nq)
8070 IF (dsm .GT. 1.0d0)
GO TO 500
8086 IF (newt .EQ. 0) nsfi = nsfi + 1
8087 IF (newt .GT. 0 .AND. stifr .LT. 1.5d0) newt = 0
8092 470 yh(i,j) = yh(i,j) + el(j)*acor(i)
8094 IF (ialth .EQ. 0)
GO TO 520
8095 IF (ialth .GT. 1)
GO TO 700
8096 IF (l .EQ. lmax)
GO TO 700
8098 490 yh(i,lmax) = acor(i)
8107 500 kflag = kflag - 1
8114 510 yh1(i) = yh1(i) - yh1(i+nyh)
8117 IF (abs(h) .LE. hmin*1.00001d0)
GO TO 660
8118 IF (kflag .LE. -3)
GO TO 640
8132 IF (l .EQ. lmax)
GO TO 540
8134 530 savf(i) = acor(i) - yh(i,lmax)
8135 dup = dvnorm(n, savf, ewt)/tesco(3,nq)
8137 rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0)
8139 rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
8141 IF (nq .EQ. 1)
GO TO 560
8142 ddn = dvnorm(n, yh(1,l), ewt)/tesco(1,nq)
8144 rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
8145 560
IF (rhsm .GE. rhup)
GO TO 570
8146 IF (rhup .GT. rhdn)
GO TO 590
8148 570
IF (rhsm .LT. rhdn)
GO TO 580
8154 IF (kflag .LT. 0 .AND. rh .GT. 1.0d0) rh = 1.0d0
8158 IF (rh .LT. 1.1d0)
GO TO 610
8161 600 yh(i,newq+1) = acor(i)*r
8165 620
IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1d0))
GO TO 610
8166 IF (kflag .LE. -2) rh = min(rh,0.2d0)
8172 IF (newq .EQ. nq)
GO TO 170
8186 640
IF (kflag .EQ. -10)
GO TO 660
8188 rh = max(hmin/abs(h),rh)
8192 CALL f (neq, tn, y, savf)
8195 650 yh(i,2) = h*savf(i)
8198 IF (nq .EQ. 1)
GO TO 200
8214 700 r = 1.0d0/tesco(2,nqu)
8216 710 acor(i) = acor(i)*r
8223 SUBROUTINE dsetpk (NEQ, Y, YSV, EWT, FTEM, SAVF, JOK, WM, IWM,
8226 INTEGER neq, jok, iwm
8227 DOUBLE PRECISION y, ysv, ewt, ftem, savf, wm
8228 dimension neq(*), y(*), ysv(*), ewt(*), ftem(*), savf(*),
8230 INTEGER iownd, iowns,
8231 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
8232 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
8233 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
8234 INTEGER jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
8235 1 nni, nli, nps, ncfn, ncfl
8236 DOUBLE PRECISION rowns,
8237 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
8238 DOUBLE PRECISION delt, epcon, sqrtn, rsqrtn
8239 COMMON /dls001/ rowns(209),
8240 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
8241 2 iownd(6), iowns(6),
8242 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
8243 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
8244 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
8245 COMMON /dlpk01/ delt, epcon, sqrtn, rsqrtn,
8246 1 jpre, jacflg, locwp, lociwp, lsavx, kmp, maxl, mnewt,
8247 2 nni, nli, nps, ncfn, ncfl
8274 DOUBLE PRECISION hl0
8278 IF (jok .EQ. -1) jcur = 1
8280 CALL jac (f, neq, tn, y, ysv, ewt, savf, ftem, hl0, jok,
8281 1 wm(locwp), iwm(lociwp), ier)
8283 IF (ier .EQ. 0)
RETURN 8289 SUBROUTINE dsrckr (RSAV, ISAV, JOB)
8303 INTEGER ils, ils2, ilsr, ilsp
8304 INTEGER i, ioff, lenilp, lenrlp, lenils, lenrls, lenilr, lenrlr
8305 DOUBLE PRECISION rsav, rls, rls2, rlsr, rlsp
8306 dimension rsav(*), isav(*)
8307 SAVE lenrls, lenils, lenrlp, lenilp, lenrlr, lenilr
8308 COMMON /dls001/ rls(218), ils(37)
8309 COMMON /dls002/ rls2, ils2(4)
8310 COMMON /dlsr01/ rlsr(5), ilsr(9)
8311 COMMON /dlpk01/ rlsp(4), ilsp(13)
8312 DATA lenrls/218/, lenils/37/, lenrlp/4/, lenilp/13/
8313 DATA lenrlr/5/, lenilr/9/
8315 IF (job .EQ. 2)
GO TO 100
8316 CALL dcopy (lenrls, rls, 1, rsav, 1)
8317 rsav(lenrls+1) = rls2
8318 CALL dcopy (lenrlr, rlsr, 1, rsav(lenrls+2), 1)
8319 CALL dcopy (lenrlp, rlsp, 1, rsav(lenrls+lenrlr+2), 1)
8322 isav(lenils+1) = ils2(1)
8323 isav(lenils+2) = ils2(2)
8324 isav(lenils+3) = ils2(3)
8325 isav(lenils+4) = ils2(4)
8328 30 isav(ioff+i) = ilsr(i)
8329 ioff = ioff + lenilr
8331 40 isav(ioff+i) = ilsp(i)
8335 CALL dcopy (lenrls, rsav, 1, rls, 1)
8336 rls2 = rsav(lenrls+1)
8337 CALL dcopy (lenrlr, rsav(lenrls+2), 1, rlsr, 1)
8338 CALL dcopy (lenrlp, rsav(lenrls+lenrlr+2), 1, rlsp, 1)
8340 120 ils(i) = isav(i)
8341 ils2(1) = isav(lenils+1)
8342 ils2(2) = isav(lenils+2)
8343 ils2(3) = isav(lenils+3)
8344 ils2(4) = isav(lenils+4)
8347 130 ilsr(i) = isav(ioff+i)
8348 ioff = ioff + lenilr
8350 140 ilsp(i) = isav(ioff+i)
8355 SUBROUTINE dainvg (RES, ADDA, NEQ, T, Y, YDOT, MITER,
8356 1 ML, MU, PW, IPVT, IER )
8358 INTEGER neq, miter, ml, mu, ipvt, ier
8359 INTEGER i, lenpw, mlp1, nrowpw
8360 DOUBLE PRECISION t, y, ydot, pw
8361 dimension y(*), ydot(*), pw(*), ipvt(*)
8374 IF (miter .GE. 4)
GO TO 100
8383 CALL res ( neq, t, y, pw, ydot, ier )
8384 IF (ier .GT. 1)
RETURN 8386 CALL adda ( neq, t, y, 0, 0, pw, neq )
8387 CALL dgefa ( pw, neq, neq, ipvt, ier )
8388 IF (ier .EQ. 0)
GO TO 20
8391 20
CALL dgesl ( pw, neq, neq, ipvt, ydot, 0 )
8397 nrowpw = 2*ml + mu + 1
8398 lenpw = neq * nrowpw
8403 CALL res ( neq, t, y, pw, ydot, ier )
8404 IF (ier .GT. 1)
RETURN 8407 CALL adda ( neq, t, y, ml, mu, pw(mlp1), nrowpw )
8408 CALL dgbfa ( pw, nrowpw, neq, ml, mu, ipvt, ier )
8409 IF (ier .EQ. 0)
GO TO 120
8412 120
CALL dgbsl ( pw, nrowpw, neq, ml, mu, ipvt, ydot, 0 )
8417 SUBROUTINE dstodi (NEQ, Y, YH, NYH, YH1, EWT, SAVF, SAVR,
8418 1 ACOR, WM, IWM, RES, ADDA, JAC, PJAC, SLVS )
8419 EXTERNAL res, adda, jac, pjac, slvs
8420 INTEGER neq, nyh, iwm
8421 DOUBLE PRECISION y, yh, yh1, ewt, savf, savr, acor, wm
8422 dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
8423 1 savr(*), acor(*), wm(*), iwm(*)
8424 INTEGER iownd, ialth, ipup, lmax, meo, nqnyh, nslp,
8425 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
8426 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
8427 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
8428 DOUBLE PRECISION conit, crate, el, elco, hold, rmax, tesco,
8429 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
8430 COMMON /dls001/ conit, crate, el(13), elco(13,12),
8431 1 hold, rmax, tesco(3,12),
8432 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
8433 3 iownd(6), ialth, ipup, lmax, meo, nqnyh, nslp,
8434 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
8435 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
8436 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
8437 INTEGER i, i1, iredo, ires, iret, j, jb, kgo, m, ncf, newq
8438 DOUBLE PRECISION dcon, ddn, del, delp, dsm, dup,
8439 1 eljh, el1h, exdn, exsm, exup,
8440 2 r, rh, rhdn, rhsm, rhup, told, dvnorm
8527 IF (jstart .GT. 0)
GO TO 200
8528 IF (jstart .EQ. -1)
GO TO 100
8529 IF (jstart .EQ. -2)
GO TO 160
8567 IF (ialth .EQ. 1) ialth = 2
8568 IF (meth .EQ. meo)
GO TO 110
8569 CALL dcfode (meth, elco, tesco)
8571 IF (nq .GT. maxord)
GO TO 120
8575 110
IF (nq .LE. maxord)
GO TO 160
8579 125 el(i) = elco(i,nq)
8583 conit = 0.5d0/(nq+2)
8584 ddn = dvnorm(n, savf, ewt)/tesco(1,l)
8586 rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
8587 rh = min(rhdn,1.0d0)
8589 IF (h .EQ. hold)
GO TO 170
8590 rh = min(rh,abs(h/hold))
8598 140
CALL dcfode (meth, elco, tesco)
8600 155 el(i) = elco(i,nq)
8604 conit = 0.5d0/(nq+2)
8605 GO TO (160, 170, 200), iret
8612 160
IF (h .EQ. hold)
GO TO 200
8617 170 rh = max(rh,hmin/abs(h))
8618 175 rh = min(rh,rmax)
8619 rh = rh/max(1.0d0,abs(h)*hmxi*rh)
8624 180 yh(i,j) = yh(i,j)*r
8628 IF (iredo .EQ. 0)
GO TO 690
8637 200
IF (abs(rc-1.0d0) .GT. ccmax) ipup = miter
8638 IF (nst .GE. nslp+msbp) ipup = miter
8645 210 yh1(i) = yh1(i) + yh1(i+nyh)
8655 savf(i) = yh(i,2) / h
8657 IF (ipup .LE. 0)
GO TO 240
8663 CALL pjac (neq, y, yh, nyh, ewt, acor, savr, savf, wm, iwm,
8669 IF (ierpj .EQ. 0)
GO TO 250
8670 IF (ierpj .LT. 0)
GO TO 435
8672 GO TO (430, 435, 430), ires
8675 CALL res ( neq, tn, y, savf, savr, ires )
8678 GO TO ( 250, 435, 430 ) , kgo
8686 CALL slvs (wm, iwm, savr, savf)
8687 IF (iersl .LT. 0)
GO TO 430
8688 IF (iersl .GT. 0)
GO TO 410
8690 del = dvnorm(n, savr, ewt) * abs(h)
8692 acor(i) = acor(i) + savr(i)
8693 savf(i) = acor(i) + yh(i,2)/h
8694 380 y(i) = yh(i,1) + el1h*acor(i)
8699 IF (m .NE. 0) crate = max(0.2d0*crate,del/delp)
8700 dcon = del*min(1.0d0,1.5d0*crate)/(tesco(2,nq)*conit)
8701 IF (dcon .LE. 1.0d0)
GO TO 460
8703 IF (m .EQ. maxcor)
GO TO 410
8704 IF (m .GE. 2 .AND. del .GT. 2.0d0*delp)
GO TO 410
8707 CALL res ( neq, tn, y, savf, savr, ires )
8710 GO TO ( 270, 435, 410 ) , kgo
8720 IF (jcur .EQ. 1)
GO TO 430
8732 440 yh1(i) = yh1(i) - yh1(i+nyh)
8734 IF (ires .EQ. 2)
GO TO 680
8735 IF (ierpj .LT. 0 .OR. iersl .LT. 0)
GO TO 685
8736 IF (abs(h) .LE. hmin*1.00001d0)
GO TO 450
8737 IF (ncf .EQ. mxncf)
GO TO 450
8742 450
IF (ires .EQ. 3)
GO TO 680
8751 IF (m .EQ. 0) dsm = del/tesco(2,nq)
8752 IF (m .GT. 0) dsm = abs(h) * dvnorm(n, acor, ewt)/tesco(2,nq)
8753 IF (dsm .GT. 1.0d0)
GO TO 500
8772 470 yh(i,j) = yh(i,j) + eljh*acor(i)
8774 IF (ialth .EQ. 0)
GO TO 520
8775 IF (ialth .GT. 1)
GO TO 700
8776 IF (l .EQ. lmax)
GO TO 700
8778 490 yh(i,lmax) = acor(i)
8787 500 kflag = kflag - 1
8794 510 yh1(i) = yh1(i) - yh1(i+nyh)
8797 IF (abs(h) .LE. hmin*1.00001d0)
GO TO 660
8798 IF (kflag .LE. -7)
GO TO 660
8812 IF (l .EQ. lmax)
GO TO 540
8814 530 savf(i) = acor(i) - yh(i,lmax)
8815 dup = abs(h) * dvnorm(n, savf, ewt)/tesco(3,nq)
8817 rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0)
8819 rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0)
8821 IF (nq .EQ. 1)
GO TO 560
8822 ddn = dvnorm(n, yh(1,l), ewt)/tesco(1,nq)
8824 rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0)
8825 560
IF (rhsm .GE. rhup)
GO TO 570
8826 IF (rhup .GT. rhdn)
GO TO 590
8828 570
IF (rhsm .LT. rhdn)
GO TO 580
8834 IF (kflag .LT. 0 .AND. rh .GT. 1.0d0) rh = 1.0d0
8838 IF (rh .LT. 1.1d0)
GO TO 610
8841 600 yh(i,newq+1) = acor(i)*r
8845 620
IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1d0))
GO TO 610
8846 IF (kflag .LE. -2) rh = min(rh,0.1d0)
8852 IF (newq .EQ. nq)
GO TO 170
8865 680 kflag = -1 - ires
8870 700 r = h/tesco(2,nqu)
8872 710 acor(i) = acor(i)*r
8879 SUBROUTINE dprepji (NEQ, Y, YH, NYH, EWT, RTEM, SAVR, S, WM, IWM,
8881 EXTERNAL res, jac, adda
8882 INTEGER neq, nyh, iwm
8883 DOUBLE PRECISION y, yh, ewt, rtem, savr, s, wm
8884 dimension neq(*), y(*), yh(nyh,*), ewt(*), rtem(*),
8885 1 s(*), savr(*), wm(*), iwm(*)
8886 INTEGER iownd, iowns,
8887 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
8888 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
8889 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
8890 DOUBLE PRECISION rowns,
8891 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
8892 COMMON /dls001/ rowns(209),
8893 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
8894 2 iownd(6), iowns(6),
8895 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
8896 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
8897 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
8898 INTEGER i, i1, i2, ier, ii, ires, j, j1, jj, lenp,
8899 1 mba, mband, meb1, meband, ml, ml3, mu
8900 DOUBLE PRECISION con, fac, hl0, r, srur, yi, yj, yjj
8941 GO TO (100, 200, 300, 400, 500), miter
8944 CALL res (neq, tn, y, s, savr, ires)
8946 IF (ires .GT. 1)
GO TO 600
8950 CALL jac ( neq, tn, y, s, 0, 0, wm(3), n )
8953 120 wm(i+2) = wm(i+2)*con
8958 CALL res (neq, tn, y, s, savr, ires)
8960 IF (ires .GT. 1)
GO TO 600
8965 r = max(srur*abs(yj),0.01d0/ewt(j))
8968 CALL res ( neq, tn, y, s, rtem, ires )
8970 IF (ires .GT. 1)
GO TO 600
8972 220 wm(i+j1) = (rtem(i) - savr(i))*fac
8977 CALL res (neq, tn, y, s, savr, ires)
8979 IF (ires .GT. 1)
GO TO 600
8982 CALL adda(neq, tn, y, 0, 0, wm(3), n)
8984 CALL dgefa (wm(3), n, n, iwm(21), ier)
8985 IF (ier .NE. 0) ierpj = 1
8991 CALL res (neq, tn, y, s, savr, ires)
8993 IF (ires .GT. 1)
GO TO 600
9002 CALL jac ( neq, tn, y, s, ml, mu, wm(ml3), meband)
9005 420 wm(i+2) = wm(i+2)*con
9010 CALL res (neq, tn, y, s, savr, ires)
9012 IF (ires .GT. 1)
GO TO 600
9022 DO 530 i = j,n,mband
9024 r = max(srur*abs(yi),0.01d0/ewt(i))
9026 CALL res ( neq, tn, y, s, rtem, ires)
9028 IF (ires .GT. 1)
GO TO 600
9029 DO 550 jj = j,n,mband
9032 r = max(srur*abs(yjj),0.01d0/ewt(jj))
9036 ii = jj*meb1 - ml + 2
9038 540 wm(ii+i) = (rtem(i) - savr(i))*fac
9042 CALL res (neq, tn, y, s, savr, ires)
9044 IF (ires .GT. 1)
GO TO 600
9047 CALL adda(neq, tn, y, ml, mu, wm(ml3), meband)
9049 CALL dgbfa (wm(3), meband, n, ml, mu, iwm(21), ier)
9050 IF (ier .NE. 0) ierpj = 1
9058 SUBROUTINE daigbt (RES, ADDA, NEQ, T, Y, YDOT,
9059 1 MB, NB, PW, IPVT, IER )
9061 INTEGER neq, mb, nb, ipvt, ier
9062 INTEGER i, lenpw, lblox, lpb, lpc
9063 DOUBLE PRECISION t, y, ydot, pw
9064 dimension y(*), ydot(*), pw(*), ipvt(*), neq(*)
9084 CALL res (neq, t, y, pw, ydot, ier)
9085 IF (ier .GT. 1)
RETURN 9086 CALL adda (neq, t, y, mb, nb, pw(1), pw(lpb), pw(lpc) )
9087 CALL ddecbt (mb, nb, pw, pw(lpb), pw(lpc), ipvt, ier)
9088 IF (ier .EQ. 0)
GO TO 20
9091 20
CALL dsolbt (mb, nb, pw, pw(lpb), pw(lpc), ydot, ipvt)
9096 SUBROUTINE dpjibt (NEQ, Y, YH, NYH, EWT, RTEM, SAVR, S, WM, IWM,
9098 EXTERNAL res, jac, adda
9099 INTEGER neq, nyh, iwm
9100 DOUBLE PRECISION y, yh, ewt, rtem, savr, s, wm
9101 dimension neq(*), y(*), yh(nyh,*), ewt(*), rtem(*),
9102 1 s(*), savr(*), wm(*), iwm(*)
9103 INTEGER iownd, iowns,
9104 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
9105 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
9106 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
9107 DOUBLE PRECISION rowns,
9108 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
9109 COMMON /dls001/ rowns(209),
9110 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
9111 2 iownd(6), iowns(6),
9112 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
9113 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
9114 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
9115 INTEGER i, ier, iia, iib, iic, ipa, ipb, ipc, ires, j, j1, j2,
9116 1 k, k1, lenp, lblox, lpb, lpc, mb, mbsq, mwid, nb
9117 DOUBLE PRECISION con, fac, hl0, r, srur
9163 GO TO (100, 200), miter
9166 CALL res (neq, tn, y, s, savr, ires)
9168 IF (ires .GT. 1)
GO TO 600
9171 CALL jac (neq, tn, y, s, mb, nb, wm(3), wm(lpb), wm(lpc))
9174 120 wm(i+2) = wm(i+2)*con
9180 CALL res (neq, tn, y, s, savr, ires)
9182 IF (ires .GT. 1)
GO TO 600
9191 DO 210 i = j1,n,mwid
9192 r = max(srur*abs(y(i)),0.01d0/ewt(i))
9195 CALL res (neq, tn, y, s, rtem, ires)
9197 IF (ires .GT. 1)
GO TO 600
9199 215 rtem(i) = rtem(i) - savr(i)
9201 DO 230 i = j1,n,mwid
9204 r = max(srur*abs(y(i)),0.01d0/ewt(i))
9208 ipa = 2 + (j-1)*mb + (k1-1)*mbsq
9210 221 wm(ipa+j2) = rtem(iia+j2)*fac
9211 IF (k1 .LE. 1)
GO TO 223
9214 ipb = ipa + lblox - mbsq
9216 222 wm(ipb+j2) = rtem(iib+j2)*fac
9218 IF (k1 .GE. nb)
GO TO 225
9221 ipc = ipa + 2*lblox + mbsq
9223 224 wm(ipc+j2) = rtem(iic+j2)*fac
9225 IF (k1 .NE. 3)
GO TO 227
9227 ipc = ipa - 2*mbsq + 2*lblox
9229 226 wm(ipc+j2) = rtem(j2)*fac
9231 IF (k1 .NE. nb-2)
GO TO 229
9234 ipb = ipa + 2*mbsq + lblox
9236 228 wm(ipb+j2) = rtem(iib+j2)*fac
9243 CALL res (neq, tn, y, s, savr, ires)
9245 IF (ires .GT. 1)
GO TO 600
9248 CALL adda (neq, tn, y, mb, nb, wm(3), wm(lpb), wm(lpc))
9250 CALL ddecbt (mb, nb, wm(3), wm(lpb), wm(lpc), iwm(21), ier)
9251 IF (ier .NE. 0) ierpj = 1
9259 SUBROUTINE dslsbt (WM, IWM, X, TEM)
9261 INTEGER lblox, lpb, lpc, mb, nb
9262 DOUBLE PRECISION wm, x, tem
9263 dimension wm(*), iwm(*), x(*), tem(*)
9283 CALL dsolbt (mb, nb, wm(3), wm(lpb), wm(lpc), x, iwm(21))
9288 SUBROUTINE ddecbt (M, N, A, B, C, IP, IER)
9289 INTEGER m, n, ip(m,n), ier
9290 DOUBLE PRECISION a(m,m,n), b(m,m,n), c(m,m,n)
9336 INTEGER nm1, nm2, km1, i, j, k
9337 DOUBLE PRECISION dp, ddot
9338 IF (m .LT. 1 .OR. n .LT. 4)
GO TO 210
9342 CALL dgefa (a, m, m, ip, ier)
9344 IF (ier .NE. 0)
GO TO 200
9346 CALL dgesl (a, m, m, ip, b(1,j,1), 0)
9347 CALL dgesl (a, m, m, ip, c(1,j,1), 0)
9352 dp = ddot(m, c(i,1,2), m, c(1,j,1), 1)
9353 b(i,j,2) = b(i,j,2) - dp
9361 dp = ddot(m, c(i,1,k), m, b(1,j,km1), 1)
9362 a(i,j,k) = a(i,j,k) - dp
9365 CALL dgefa (a(1,1,k), m, m, ip(1,k), ier)
9366 IF (ier .NE. 0)
GO TO 200
9368 80
CALL dgesl (a(1,1,k), m, m, ip(1,k), b(1,j,k), 0)
9373 dp = ddot(m, b(i,1,n), m, b(1,j,nm2), 1)
9374 c(i,j,n) = c(i,j,n) - dp
9379 dp = ddot(m, c(i,1,n), m, b(1,j,nm1), 1)
9380 a(i,j,n) = a(i,j,n) - dp
9383 CALL dgefa (a(1,1,n), m, m, ip(1,n), ier)
9385 IF (ier .NE. 0)
GO TO 200
9395 SUBROUTINE dsolbt (M, N, A, B, C, Y, IP)
9396 INTEGER m, n, ip(m,n)
9397 DOUBLE PRECISION a(m,m,n), b(m,m,n), c(m,m,n), y(m,n)
9417 INTEGER nm1, nm2, i, k, kb, km1, kp1
9418 DOUBLE PRECISION dp, ddot
9422 CALL dgesl (a, m, m, ip, y, 0)
9426 dp = ddot(m, c(i,1,k), m, y(1,km1), 1)
9427 y(i,k) = y(i,k) - dp
9429 CALL dgesl (a(1,1,k), m, m, ip(1,k), y(1,k), 0)
9432 dp = ddot(m, c(i,1,n), m, y(1,nm1), 1)
9433 1 + ddot(m, b(i,1,n), m, y(1,nm2), 1)
9434 y(i,n) = y(i,n) - dp
9436 CALL dgesl (a(1,1,n), m, m, ip(1,n), y(1,n), 0)
9442 dp = ddot(m, b(i,1,k), m, y(1,kp1), 1)
9443 y(i,k) = y(i,k) - dp
9447 dp = ddot(m, c(i,1,1), m, y(1,3), 1)
9448 y(i,1) = y(i,1) - dp
9454 SUBROUTINE diprepi (NEQ, Y, S, RWORK, IA, JA, IC, JC, IPFLAG,
9456 EXTERNAL res, jac, adda
9457 INTEGER neq, ia, ja, ic, jc, ipflag
9458 DOUBLE PRECISION y, s, rwork
9459 dimension neq(*), y(*), s(*), rwork(*), ia(*), ja(*), ic(*), jc(*)
9460 INTEGER iownd, iowns,
9461 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
9462 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
9463 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
9464 INTEGER iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
9465 1 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
9466 2 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
9467 3 nslj, ngp, nlu, nnz, nsp, nzl, nzu
9468 DOUBLE PRECISION rowns,
9469 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
9470 DOUBLE PRECISION rlss
9471 COMMON /dls001/ rowns(209),
9472 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
9473 2 iownd(6), iowns(6),
9474 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
9475 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
9476 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
9477 COMMON /dlss01/ rlss(6),
9478 1 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
9479 2 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
9480 3 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
9481 4 nslj, ngp, nlu, nnz, nsp, nzl, nzu
9482 INTEGER i, imax, lewtn, lyhd, lyhn
9497 CALL dprepi (neq, y, s, rwork(lyh), rwork(lsavf), rwork(lewt),
9498 1 rwork(lacor), ia, ja, ic, jc, rwork(lwm), rwork(lwm), ipflag,
9500 lenwk = max(lreq,lwmin)
9501 IF (ipflag .LT. 0)
RETURN 9504 IF (lyhn .GT. lyh)
RETURN 9506 IF (lyhd .EQ. 0)
GO TO 20
9507 imax = lyhn - 1 + lenyhm
9509 10 rwork(i) = rwork(i+lyhd)
9512 20 lsavf = lyh + lenyh
9515 IF (istatc .EQ. 3)
GO TO 40
9517 IF (lewtn .GT. lewt)
RETURN 9519 30 rwork(i+lewtn-1) = rwork(i+lewt-1)
9525 SUBROUTINE dprepi (NEQ, Y, S, YH, SAVR, EWT, RTEM, IA, JA, IC, JC,
9526 1 WK, IWK, IPPER, RES, JAC, ADDA)
9527 EXTERNAL res, jac, adda
9528 INTEGER neq, ia, ja, ic, jc, iwk, ipper
9529 DOUBLE PRECISION y, s, yh, savr, ewt, rtem, wk
9530 dimension neq(*), y(*), s(*), yh(*), savr(*), ewt(*), rtem(*),
9531 1 ia(*), ja(*), ic(*), jc(*), wk(*), iwk(*)
9532 INTEGER iownd, iowns,
9533 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
9534 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
9535 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
9536 INTEGER iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
9537 1 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
9538 2 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
9539 3 nslj, ngp, nlu, nnz, nsp, nzl, nzu
9540 DOUBLE PRECISION rowns,
9541 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
9542 DOUBLE PRECISION rlss
9543 COMMON /dls001/ rowns(209),
9544 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
9545 2 iownd(6), iowns(6),
9546 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
9547 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
9548 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
9549 COMMON /dlss01/ rlss(6),
9550 1 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
9551 2 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
9552 3 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
9553 4 nslj, ngp, nlu, nnz, nsp, nzl, nzu
9554 INTEGER i, ibr, ier, ipil, ipiu, iptt1, iptt2, j, k, knew, kamax,
9555 1 kamin, kcmax, kcmin, ldif, lenigp, lenwk1, liwk, ljfo, maxg,
9557 DOUBLE PRECISION erwt, fac, yj
9605 IF (moss .EQ. 0) liwk = liwk - n
9606 IF (moss .EQ. 1 .OR. moss .EQ. 2) liwk = lenwk1*lrat
9607 IF (ipjan+n-1 .GT. liwk)
GO TO 310
9608 IF (moss .EQ. 0)
GO TO 30
9610 IF (istatc .EQ. 3)
GO TO 20
9615 fac = 1.0d0 + 1.0d0/(i + 1.0d0)
9616 y(i) = y(i) + fac*sign(erwt,y(i))
9617 s(i) = 1.0d0 + fac*erwt
9619 GO TO (70, 100, 150, 200), moss
9626 GO TO (70, 100, 150, 200), moss
9637 IF (kamin .GT. kamax)
GO TO 45
9638 DO 40 k = kamin,kamax
9641 IF (knew .GT. liwk)
GO TO 310
9645 45 kamin = kamax + 1
9647 IF (kcmin .GT. kcmax)
GO TO 55
9648 DO 50 k = kcmin,kcmax
9650 IF (iwk(liwk+i) .NE. 0)
GO TO 50
9651 IF (knew .GT. liwk)
GO TO 310
9655 55 iwk(ipian+j) = knew + 1 - ipjan
9664 CALL res (neq, tn, y, s, savr, ier)
9665 IF (ier .GT. 1)
GO TO 370
9668 75 wk(lenwk1+i) = 0.0d0
9672 CALL adda (neq, tn, y, j, iwk(ipian), iwk(ipjan), wk(lenwk1+1))
9673 CALL jac (neq, tn, y, s, j, iwk(ipian), iwk(ipjan), savr)
9676 IF (wk(ljfo) .EQ. 0.0d0)
GO TO 80
9680 80
IF (savr(i) .EQ. 0.0d0)
GO TO 90
9682 85
IF (k .GT. liwk)
GO TO 310
9686 iwk(ipian+j) = k + 1 - ipjan
9692 105 wk(lenwk1+i) = 0.0d0
9696 IF (miter .EQ. 1) ier = 1
9697 CALL res (neq, tn, y, s, savr, ier)
9698 IF (ier .GT. 1)
GO TO 370
9700 CALL adda (neq, tn, y, j, iwk(ipian), iwk(ipjan), wk(lenwk1+1))
9703 y(j) = yj + sign(erwt,yj)
9704 CALL res (neq, tn, y, s, rtem, ier)
9705 IF (ier .GT. 1)
RETURN 9709 IF (wk(ljfo) .EQ. 0.0d0)
GO TO 110
9712 110
IF (rtem(i) .EQ. savr(i))
GO TO 120
9713 115
IF (k .GT. liwk)
GO TO 310
9717 iwk(ipian+j) = k + 1 - ipjan
9725 CALL res (neq, tn, y, s, savr, ier)
9726 IF (ier .GT. 1)
GO TO 370
9733 CALL jac (neq, tn, y, s, j, iwk(ipian), iwk(ipjan), savr)
9735 IF (kamin .GT. kamax)
GO TO 170
9736 DO 160 k = kamin,kamax
9739 IF (knew .GT. liwk)
GO TO 310
9743 170 kamin = kamax + 1
9745 IF (savr(i) .EQ. 0.0d0)
GO TO 180
9747 IF (knew .GT. liwk)
GO TO 310
9751 iwk(ipian+j) = knew + 1 - ipjan
9760 IF (miter .EQ. 1) ier = 1
9761 CALL res (neq, tn, y, s, savr, ier)
9762 IF (ier .GT. 1)
GO TO 370
9766 y(j) = yj + sign(erwt,yj)
9767 CALL res (neq, tn, y, s, rtem, ier)
9768 IF (ier .GT. 1)
RETURN 9771 IF (kamin .GT. kamax)
GO TO 225
9772 DO 220 k = kamin,kamax
9775 IF (knew .GT. liwk)
GO TO 310
9779 225 kamin = kamax + 1
9781 IF (rtem(i) .EQ. savr(i))
GO TO 230
9782 IF (knew .GT. liwk)
GO TO 310
9786 iwk(ipian+j) = knew + 1 - ipjan
9790 IF (moss .EQ. 0 .OR. istatc .EQ. 3)
GO TO 250
9794 250 nnz = iwk(ipian+n) - 1
9799 IF (miter .NE. 2)
GO TO 260
9809 lreq = iptt2 + n - 1
9810 IF (lreq .GT. liwk)
GO TO 320
9811 CALL jgroup (n, iwk(ipian), iwk(ipjan), maxg, ngp, iwk(ipigp),
9812 1 iwk(ipjgp), iwk(iptt1), iwk(iptt2), ier)
9813 IF (ier .NE. 0)
GO TO 320
9817 260 ipr = ipigp + lenigp
9821 iprsp = (ipisp-2)/lrat + 2
9822 iesp = lenwk + 1 - iprsp
9823 IF (iesp .LT. 0)
GO TO 330
9827 nsp = liwk + 1 - ipisp
9828 CALL odrv(n, iwk(ipian), iwk(ipjan), wk, iwk(ipr), iwk(ipic), nsp,
9829 1 iwk(ipisp), 1, iys)
9830 IF (iys .EQ. 11*n+1)
GO TO 340
9831 IF (iys .NE. 0)
GO TO 330
9834 ipa = lenwk + 1 - nnz
9836 lreq = max(12*n/lrat, 6*n/lrat+2*n+nnz) + 3
9837 lreq = lreq + iprsp - 1 + nnz
9838 IF (lreq .GT. lenwk)
GO TO 350
9841 280 wk(iba+i) = 0.0d0
9842 ipisp = lrat*(iprsp - 1) + 1
9843 CALL cdrv(n,iwk(ipr),iwk(ipc),iwk(ipic),iwk(ipian),iwk(ipjan),
9844 1 wk(ipa),wk(ipa),wk(ipa),nsp,iwk(ipisp),wk(iprsp),iesp,5,iys)
9846 IF (iys .EQ. 10*n+1)
GO TO 350
9847 IF (iys .NE. 0)
GO TO 360
9849 ipiu = ipil + 2*n + 1
9850 nzu = iwk(ipil+n) - iwk(ipil)
9851 nzl = iwk(ipiu+n) - iwk(ipiu)
9852 IF (lrat .GT. 1)
GO TO 290
9853 CALL adjlr (n, iwk(ipisp), ldif)
9856 IF (lrat .EQ. 2 .AND. nnz .EQ. n) lreq = lreq + 1
9857 nsp = nsp + lreq - lenwk
9858 ipa = lreq + 1 - nnz
9864 lreq = 2 + (2*n + 1)/lrat
9865 lreq = max(lenwk+1,lreq)
9869 lreq = (lreq - 1)/lrat + 1
9873 CALL cntnzu (n, iwk(ipian), iwk(ipjan), nzsut)
9874 lreq = lenwk - iesp + (3*n + 4*nzsut - 1)/lrat + 1
9887 370 ipper = -ier - 5
9888 lreq = 2 + (2*n + 1)/lrat
9893 SUBROUTINE dainvgs (NEQ, T, Y, WK, IWK, TEM, YDOT, IER, RES, ADDA)
9895 INTEGER neq, iwk, ier
9896 INTEGER iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
9897 1 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
9898 2 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
9899 3 nslj, ngp, nlu, nnz, nsp, nzl, nzu
9900 INTEGER i, imul, j, k, kmin, kmax
9901 DOUBLE PRECISION t, y, wk, tem, ydot
9902 DOUBLE PRECISION rlss
9903 dimension y(*), wk(*), iwk(*), tem(*), ydot(*)
9904 COMMON /dlss01/ rlss(6),
9905 1 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
9906 2 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
9907 3 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
9908 4 nslj, ngp, nlu, nnz, nsp, nzl, nzu
9942 10 wk(iba+i) = 0.0d0
9945 CALL res (neq, t, y, wk(ipa), ydot, ier)
9946 IF (ier .GT. 1)
RETURN 9950 kmax = iwk(ipian+j) - 1
9954 CALL adda (neq, t, y, j, iwk(ipian), iwk(ipjan), tem)
9957 20 wk(iba+k) = tem(i)
9966 CALL cdrv (neq,iwk(ipr),iwk(ipc),iwk(ipic),iwk(ipian),iwk(ipjan),
9967 1 wk(ipa),tem,tem,nsp,iwk(ipisp),wk(iprsp),iesp,2,iys)
9968 IF (iys .EQ. 0)
GO TO 50
9969 imul = (iys - 1)/neq
9971 IF (imul .EQ. 8) ier = 1
9972 IF (imul .EQ. 10) ier = 4
9976 50
CALL cdrv (neq,iwk(ipr),iwk(ipc),iwk(ipic),iwk(ipian),iwk(ipjan),
9977 1 wk(ipa),ydot,ydot,nsp,iwk(ipisp),wk(iprsp),iesp,4,iys)
9978 IF (iys .NE. 0) ier = 5
9983 SUBROUTINE dprjis (NEQ, Y, YH, NYH, EWT, RTEM, SAVR, S, WK, IWK,
9985 EXTERNAL res, jac, adda
9986 INTEGER neq, nyh, iwk
9987 DOUBLE PRECISION y, yh, ewt, rtem, savr, s, wk
9988 dimension neq(*), y(*), yh(nyh,*), ewt(*), rtem(*),
9989 1 s(*), savr(*), wk(*), iwk(*)
9990 INTEGER iownd, iowns,
9991 1 icf, ierpj, iersl, jcur, jstart, kflag, l,
9992 2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
9993 3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
9994 INTEGER iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
9995 1 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
9996 2 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
9997 3 nslj, ngp, nlu, nnz, nsp, nzl, nzu
9998 DOUBLE PRECISION rowns,
9999 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
10000 DOUBLE PRECISION rlss
10001 COMMON /dls001/ rowns(209),
10002 1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
10003 2 iownd(6), iowns(6),
10004 3 icf, ierpj, iersl, jcur, jstart, kflag, l,
10005 4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
10006 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
10007 COMMON /dlss01/ rlss(6),
10008 1 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp,
10009 2 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa,
10010 3 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj,
10011 4 nslj, ngp, nlu, nnz, nsp, nzl, nzu
10012 INTEGER i, imul, ires, j, jj, jmax, jmin, k, kmax, kmin, ng
10013 DOUBLE PRECISION con, fac, hl0, r, srur
10053 GO TO (100, 200), miter
10057 CALL res (neq, tn, y, s, savr, ires)
10059 IF (ires .GT. 1)
GO TO 600
10062 kmax = iwk(ipian+j)-1
10064 110 rtem(i) = 0.0d0
10065 CALL jac (neq, tn, y, s, j, iwk(ipian), iwk(ipjan), rtem)
10067 120 rtem(i) = rtem(i)*con
10068 CALL adda (neq, tn, y, j, iwk(ipian), iwk(ipjan), rtem)
10069 DO 125 k = kmin,kmax
10071 wk(iba+k) = rtem(i)
10080 CALL res (neq, tn, y, s, savr, ires)
10082 IF (ires .GT. 1)
GO TO 600
10086 jmax = iwk(ipigp+ng) - 1
10087 DO 210 j = jmin,jmax
10089 r = max(srur*abs(y(jj)),0.01d0/ewt(jj))
10090 210 y(jj) = y(jj) + r
10091 CALL res (neq,tn,y,s,rtem,ires)
10093 IF (ires .GT. 1)
GO TO 600
10094 DO 230 j = jmin,jmax
10097 r = max(srur*abs(y(jj)),0.01d0/ewt(jj))
10099 kmin = iwk(ibian+jj)
10100 kmax = iwk(ibian+jj+1) - 1
10101 DO 220 k = kmin,kmax
10103 rtem(i) = (rtem(i) - savr(i))*fac
10105 CALL adda (neq, tn, y, jj, iwk(ipian), iwk(ipjan), rtem)
10106 DO 225 k = kmin,kmax
10108 wk(iba+k) = rtem(i)
10114 CALL res (neq, tn, y, s, savr, ires)
10116 IF (ires .GT. 1)
GO TO 600
10122 295 rtem(i) = 0.0d0
10123 CALL cdrv (n,iwk(ipr),iwk(ipc),iwk(ipic),iwk(ipian),iwk(ipjan),
10124 1 wk(ipa),rtem,rtem,nsp,iwk(ipisp),wk(iprsp),iesp,2,iys)
10125 IF (iys .EQ. 0)
RETURN 10128 IF (imul .EQ. 8) ierpj = 1
10129 IF (imul .EQ. 10) ierpj = -1