1 SUBROUTINE ddastp (X, Y, YPRIME, NEQ, RES, JAC, H, WT, JSTART,
2 + IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, ALPHA, BETA, GAMMA,
3 + PSI, SIGMA, CJ, CJOLD, HOLD, S, HMIN, UROUND, IPHASE, JCALC,
4 + K, KOLD, NS, NONNEG, NTEMP)
94 INTEGER NEQ, JSTART, IDID, IPAR(*), IWM(*), IPHASE, JCALC, K,
95 * KOLD, NS, NONNEG, NTEMP
97 * x, y(*), yprime(*), h, wt(*), rpar(*), phi(neq,*), delta(*),
98 * e(*), wm(*), alpha(*), beta(*), gamma(*), psi(*), sigma(*), cj,
99 * cjold, hold, s, hmin, uround
102 EXTERNAL ddajac, ddanrm, ddaslv, ddatrp
103 DOUBLE PRECISION DDANRM
105 INTEGER I, IER, IRES, J, J1, KDIFF, KM1, KNEW, KP1, KP2, LCTF,
106 * letf, lmxord, lnje, lnre, lnst, m, maxit, ncf, nef, nsf, nsp1
108 * alpha0, alphas, cjlast, ck, delnrm, enorm, erk, erkm1,
109 * erkm2, erkp1, err, est, hnew, oldnrm, pnorm, r, rate, temp1,
110 * temp2, terk, terkm1, terkm2, terkp1, xold, xrate
141 IF(jstart .NE. 0)
GO TO 120
175 IF(h.NE.hold.OR.k .NE. kold) ns = 0
178 IF(kp1 .LT. ns)
GO TO 230
188 beta(i)=beta(i-1)*psi(i-1)/temp2
191 sigma(i)=(i-1)*sigma(i-1)*alpha(i)
192 gamma(i)=gamma(i-1)+alpha(i-1)/h
201 alphas = alphas - 1.0d0/i
202 alpha0 = alpha0 - alpha(i)
210 ck =
abs(alpha(kp1) + alphas - alpha0)
211 ck =
max(ck,alpha(kp1))
214 temp1 = (1.0d0 - xrate)/(1.0d0 + xrate)
216 IF (cj/cjold .LT. temp1 .OR. cj/cjold .GT. temp2) jcalc = -1
217 IF (cj .NE. cjlast) s = 100.d0
220 IF(kp1 .LT. nsp1)
GO TO 280
223 260 phi(i,j)=beta(j)*phi(i,j)
248 320 yprime(i)=yprime(i)+gamma(j)*phi(i,j)
250 pnorm = ddanrm(neq,y,wt,rpar,ipar)
258 iwm(lnre)=iwm(lnre)+1
260 CALL res(x,y,yprime,delta,ires,rpar,ipar)
261 IF (ires .LT. 0)
GO TO 380
269 IF(jcalc .NE. -1)
GO TO 340
270 iwm(lnje)=iwm(lnje)+1
272 CALL ddajac(neq,x,y,yprime,delta,cj,h,
273 * ier,wt,e,wm,iwm,res,ires,uround,jac,rpar,
277 IF (ires .LT. 0)
GO TO 380
278 IF(ier .NE. 0)
GO TO 380
292 temp1 = 2.0d0/(1.0d0 + cj/cjold)
294 355 delta(i) = delta(i) * temp1
298 CALL ddaslv(neq,delta,wm,iwm)
304 360 yprime(i)=yprime(i)-cj*delta(i)
307 delnrm=ddanrm(neq,delta,wt,rpar,ipar)
308 IF (delnrm .LE. 100.d0*uround*pnorm)
GO TO 375
309 IF (m .GT. 0)
GO TO 365
312 365 rate = (delnrm/oldnrm)**(1.0d0/m)
313 IF (rate .GT. 0.90d0)
GO TO 370
314 s = rate/(1.0d0 - rate)
315 367
IF (s*delnrm .LE. 0.33d0)
GO TO 375
322 IF(m.GE.maxit)
GO TO 370
326 iwm(lnre)=iwm(lnre)+1
328 CALL res(x,y,yprime,delta,ires,
330 IF (ires .LT. 0)
GO TO 380
339 IF(jcalc.EQ.0)
GO TO 380
348 375
IF(nonneg .EQ. 0)
GO TO 390
350 377 delta(i) =
min(y(i),0.0d0)
351 delnrm = ddanrm(neq,delta,wt,rpar,ipar)
352 IF(delnrm .GT. 0.33d0)
GO TO 380
354 378 e(i) = e(i) - delta(i)
363 IF(.NOT.convgd)
GO TO 600
378 enorm = ddanrm(neq,e,wt,rpar,ipar)
379 erk = sigma(k+1)*enorm
383 IF(k .EQ. 1)
GO TO 430
385 405 delta(i) = phi(i,kp1) + e(i)
386 erkm1=sigma(k)*ddanrm(neq,delta,wt,rpar,ipar)
388 IF(k .GT. 2)
GO TO 410
389 IF(terkm1 .LE. 0.5d0*terk)
GO TO 420
393 415 delta(i) = phi(i,k) + delta(i)
394 erkm2=sigma(k-1)*ddanrm(neq,delta,wt,rpar,ipar)
396 IF(
max(terkm1,terkm2).GT.terk)
GO TO 430
407 IF(err .GT. 1.0d0)
GO TO 600
421 iwm(lnst)=iwm(lnst)+1
432 IF(knew.EQ.km1.OR.k.EQ.iwm(lmxord))iphase=1
433 IF(iphase .EQ. 0)
GO TO 545
434 IF(knew.EQ.km1)
GO TO 540
435 IF(k.EQ.iwm(lmxord))
GO TO 550
436 IF(kp1.GE.ns.OR.kdiff.EQ.1)
GO TO 550
438 510 delta(i)=e(i)-phi(i,kp2)
439 erkp1 = (1.0d0/(k+2))*ddanrm(neq,delta,wt,rpar,ipar)
442 IF(terkp1.GE.0.5d0*terk)
GO TO 550
444 520
IF(terkm1.LE.
min(terk,terkp1))
GO TO 540
445 IF(terkp1.GE.terk.OR.k.EQ.iwm(lmxord))
GO TO 550
469 r=(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
470 IF(r .LT. 2.0d0)
GO TO 555
473 555
IF(r .GT. 1.0d0)
GO TO 560
474 r =
max(0.5d0,
min(0.9d0,r))
481 IF(kold.EQ.iwm(lmxord))
GO TO 585
486 590 phi(i,kp1)=phi(i,kp1)+e(i)
490 595 phi(i,j)=phi(i,j)+phi(i,j+1)
509 IF(kp1.LT.nsp1)
GO TO 630
513 610 phi(i,j)=temp1*phi(i,j)
517 640 psi(i-1)=psi(i)-h
523 iwm(lctf)=iwm(lctf)+1
529 IF(ier.EQ.0)
GO TO 650
538 IF (nsf .LT. 3 .AND.
abs(h) .GE. hmin)
GO TO 690
548 IF (ires .GT. -2)
GO TO 655
554 IF (ncf .LT. 10 .AND.
abs(h) .GE. hmin)
GO TO 690
556 IF (ires .LT. 0) idid = -10
557 IF (nef .GE. 3) idid = -9
565 iwm(letf)=iwm(letf)+1
566 IF (nef .GT. 1)
GO TO 665
573 r = 0.90d0*(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
574 r =
max(0.25d0,
min(0.9d0,r))
576 IF (
abs(h) .GE. hmin)
GO TO 690
583 665
IF (nef .GT. 2)
GO TO 670
586 IF (
abs(h) .GE. hmin)
GO TO 690
594 IF (
abs(h) .GE. hmin)
GO TO 690
604 CALL ddatrp(x,x,y,yprime,neq,k,phi,psi)
charNDArray max(char d, const charNDArray &m)
charNDArray min(char d, const charNDArray &m)