GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
ddastp.f
Go to the documentation of this file.
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)
5 C***BEGIN PROLOGUE DDASTP
6 C***SUBSIDIARY
7 C***PURPOSE Perform one step of the DDASSL integration.
8 C***LIBRARY SLATEC (DASSL)
9 C***TYPE DOUBLE PRECISION (SDASTP-S, DDASTP-D)
10 C***AUTHOR PETZOLD, LINDA R., (LLNL)
11 C***DESCRIPTION
12 C-----------------------------------------------------------------------
13 C DDASTP SOLVES A SYSTEM OF DIFFERENTIAL/
14 C ALGEBRAIC EQUATIONS OF THE FORM
15 C G(X,Y,YPRIME) = 0, FOR ONE STEP (NORMALLY
16 C FROM X TO X+H).
17 C
18 C THE METHODS USED ARE MODIFIED DIVIDED
19 C DIFFERENCE,FIXED LEADING COEFFICIENT
20 C FORMS OF BACKWARD DIFFERENTIATION
21 C FORMULAS. THE CODE ADJUSTS THE STEPSIZE
22 C AND ORDER TO CONTROL THE LOCAL ERROR PER
23 C STEP.
24 C
25 C
26 C THE PARAMETERS REPRESENT
27 C X -- INDEPENDENT VARIABLE
28 C Y -- SOLUTION VECTOR AT X
29 C YPRIME -- DERIVATIVE OF SOLUTION VECTOR
30 C AFTER SUCCESSFUL STEP
31 C NEQ -- NUMBER OF EQUATIONS TO BE INTEGRATED
32 C RES -- EXTERNAL USER-SUPPLIED SUBROUTINE
33 C TO EVALUATE THE RESIDUAL. THE CALL IS
34 C CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
35 C X,Y,YPRIME ARE INPUT. DELTA IS OUTPUT.
36 C ON INPUT, IRES=0. RES SHOULD ALTER IRES ONLY
37 C IF IT ENCOUNTERS AN ILLEGAL VALUE OF Y OR A
38 C STOP CONDITION. SET IRES=-1 IF AN INPUT VALUE
39 C OF Y IS ILLEGAL, AND DDASTP WILL TRY TO SOLVE
40 C THE PROBLEM WITHOUT GETTING IRES = -1. IF
41 C IRES=-2, DDASTP RETURNS CONTROL TO THE CALLING
42 C PROGRAM WITH IDID = -11.
43 C JAC -- EXTERNAL USER-SUPPLIED ROUTINE TO EVALUATE
44 C THE ITERATION MATRIX (THIS IS OPTIONAL)
45 C THE CALL IS OF THE FORM
46 C CALL JAC(X,Y,YPRIME,PD,CJ,RPAR,IPAR)
47 C PD IS THE MATRIX OF PARTIAL DERIVATIVES,
48 C PD=DG/DY+CJ*DG/DYPRIME
49 C H -- APPROPRIATE STEP SIZE FOR NEXT STEP.
50 C NORMALLY DETERMINED BY THE CODE
51 C WT -- VECTOR OF WEIGHTS FOR ERROR CRITERION.
52 C JSTART -- INTEGER VARIABLE SET 0 FOR
53 C FIRST STEP, 1 OTHERWISE.
54 C IDID -- COMPLETION CODE WITH THE FOLLOWING MEANINGS:
55 C IDID= 1 -- THE STEP WAS COMPLETED SUCCESSFULLY
56 C IDID=-6 -- THE ERROR TEST FAILED REPEATEDLY
57 C IDID=-7 -- THE CORRECTOR COULD NOT CONVERGE
58 C IDID=-8 -- THE ITERATION MATRIX IS SINGULAR
59 C IDID=-9 -- THE CORRECTOR COULD NOT CONVERGE.
60 C THERE WERE REPEATED ERROR TEST
61 C FAILURES ON THIS STEP.
62 C IDID=-10-- THE CORRECTOR COULD NOT CONVERGE
63 C BECAUSE IRES WAS EQUAL TO MINUS ONE
64 C IDID=-11-- IRES EQUAL TO -2 WAS ENCOUNTERED,
65 C AND CONTROL IS BEING RETURNED TO
66 C THE CALLING PROGRAM
67 C RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS THAT
68 C ARE USED FOR COMMUNICATION BETWEEN THE
69 C CALLING PROGRAM AND EXTERNAL USER ROUTINES
70 C THEY ARE NOT ALTERED BY DDASTP
71 C PHI -- ARRAY OF DIVIDED DIFFERENCES USED BY
72 C DDASTP. THE LENGTH IS NEQ*(K+1),WHERE
73 C K IS THE MAXIMUM ORDER
74 C DELTA,E -- WORK VECTORS FOR DDASTP OF LENGTH NEQ
75 C WM,IWM -- REAL AND INTEGER ARRAYS STORING
76 C MATRIX INFORMATION SUCH AS THE MATRIX
77 C OF PARTIAL DERIVATIVES,PERMUTATION
78 C VECTOR,AND VARIOUS OTHER INFORMATION.
79 C
80 C THE OTHER PARAMETERS ARE INFORMATION
81 C WHICH IS NEEDED INTERNALLY BY DDASTP TO
82 C CONTINUE FROM STEP TO STEP.
83 C
84 C-----------------------------------------------------------------------
85 C***ROUTINES CALLED DDAJAC, DDANRM, DDASLV, DDATRP
86 C***REVISION HISTORY (YYMMDD)
87 C 830315 DATE WRITTEN
88 C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
89 C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format.
90 C 901026 Added explicit declarations for all variables and minor
91 C cosmetic changes to prologue. (FNF)
92 C***END PROLOGUE DDASTP
93 C
94  INTEGER NEQ, JSTART, IDID, IPAR(*), IWM(*), IPHASE, JCALC, K,
95  * KOLD, NS, NONNEG, NTEMP
96  DOUBLE PRECISION
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
100  EXTERNAL res, jac
101 C
102  EXTERNAL ddajac, ddanrm, ddaslv, ddatrp
103  DOUBLE PRECISION DDANRM
104 C
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
107  double precision
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
111  LOGICAL CONVGD
112 C
113  PARAMETER (LMXORD=3)
114  parameter(lnst=11)
115  parameter(lnre=12)
116  parameter(lnje=13)
117  parameter(letf=14)
118  parameter(lctf=15)
119 C
120  DATA maxit/4/
121  DATA xrate/0.25d0/
122 C
123 C
124 C
125 C
126 C
127 C-----------------------------------------------------------------------
128 C BLOCK 1.
129 C INITIALIZE. ON THE FIRST CALL,SET
130 C THE ORDER TO 1 AND INITIALIZE
131 C OTHER VARIABLES.
132 C-----------------------------------------------------------------------
133 C
134 C INITIALIZATIONS FOR ALL CALLS
135 C***FIRST EXECUTABLE STATEMENT DDASTP
136  idid=1
137  xold=x
138  ncf=0
139  nsf=0
140  nef=0
141  IF(jstart .NE. 0) GO TO 120
142 C
143 C IF THIS IS THE FIRST STEP,PERFORM
144 C OTHER INITIALIZATIONS
145  iwm(letf) = 0
146  iwm(lctf) = 0
147  k=1
148  kold=0
149  hold=0.0d0
150  jstart=1
151  psi(1)=h
152  cjold = 1.0d0/h
153  cj = cjold
154  s = 100.d0
155  jcalc = -1
156  delnrm=1.0d0
157  iphase = 0
158  ns=0
159 120 CONTINUE
160 C
161 C
162 C
163 C
164 C
165 C-----------------------------------------------------------------------
166 C BLOCK 2
167 C COMPUTE COEFFICIENTS OF FORMULAS FOR
168 C THIS STEP.
169 C-----------------------------------------------------------------------
170 200 CONTINUE
171  kp1=k+1
172  kp2=k+2
173  km1=k-1
174  xold=x
175  IF(h.NE.hold.OR.k .NE. kold) ns = 0
176  ns=min(ns+1,kold+2)
177  nsp1=ns+1
178  IF(kp1 .LT. ns)GO TO 230
179 C
180  beta(1)=1.0d0
181  alpha(1)=1.0d0
182  temp1=h
183  gamma(1)=0.0d0
184  sigma(1)=1.0d0
185  DO 210 i=2,kp1
186  temp2=psi(i-1)
187  psi(i-1)=temp1
188  beta(i)=beta(i-1)*psi(i-1)/temp2
189  temp1=temp2+h
190  alpha(i)=h/temp1
191  sigma(i)=(i-1)*sigma(i-1)*alpha(i)
192  gamma(i)=gamma(i-1)+alpha(i-1)/h
193 210 CONTINUE
194  psi(kp1)=temp1
195 230 CONTINUE
196 C
197 C COMPUTE ALPHAS, ALPHA0
198  alphas = 0.0d0
199  alpha0 = 0.0d0
200  DO 240 i = 1,k
201  alphas = alphas - 1.0d0/i
202  alpha0 = alpha0 - alpha(i)
203 240 CONTINUE
204 C
205 C COMPUTE LEADING COEFFICIENT CJ
206  cjlast = cj
207  cj = -alphas/h
208 C
209 C COMPUTE VARIABLE STEPSIZE ERROR COEFFICIENT CK
210  ck = abs(alpha(kp1) + alphas - alpha0)
211  ck = max(ck,alpha(kp1))
212 C
213 C DECIDE WHETHER NEW JACOBIAN IS NEEDED
214  temp1 = (1.0d0 - xrate)/(1.0d0 + xrate)
215  temp2 = 1.0d0/temp1
216  IF (cj/cjold .LT. temp1 .OR. cj/cjold .GT. temp2) jcalc = -1
217  IF (cj .NE. cjlast) s = 100.d0
218 C
219 C CHANGE PHI TO PHI STAR
220  IF(kp1 .LT. nsp1) GO TO 280
221  DO 270 j=nsp1,kp1
222  DO 260 i=1,neq
223 260 phi(i,j)=beta(j)*phi(i,j)
224 270 CONTINUE
225 280 CONTINUE
226 C
227 C UPDATE TIME
228  x=x+h
229 C
230 C
231 C
232 C
233 C
234 C-----------------------------------------------------------------------
235 C BLOCK 3
236 C PREDICT THE SOLUTION AND DERIVATIVE,
237 C AND SOLVE THE CORRECTOR EQUATION
238 C-----------------------------------------------------------------------
239 C
240 C FIRST,PREDICT THE SOLUTION AND DERIVATIVE
241 300 CONTINUE
242  DO 310 i=1,neq
243  y(i)=phi(i,1)
244 310 yprime(i)=0.0d0
245  DO 330 j=2,kp1
246  DO 320 i=1,neq
247  y(i)=y(i)+phi(i,j)
248 320 yprime(i)=yprime(i)+gamma(j)*phi(i,j)
249 330 CONTINUE
250  pnorm = ddanrm(neq,y,wt,rpar,ipar)
251 C
252 C
253 C
254 C SOLVE THE CORRECTOR EQUATION USING A
255 C MODIFIED NEWTON SCHEME.
256  convgd= .true.
257  m=0
258  iwm(lnre)=iwm(lnre)+1
259  ires = 0
260  CALL res(x,y,yprime,delta,ires,rpar,ipar)
261  IF (ires .LT. 0) GO TO 380
262 C
263 C
264 C IF INDICATED,REEVALUATE THE
265 C ITERATION MATRIX PD = DG/DY + CJ*DG/DYPRIME
266 C (WHERE G(X,Y,YPRIME)=0). SET
267 C JCALC TO 0 AS AN INDICATOR THAT
268 C THIS HAS BEEN DONE.
269  IF(jcalc .NE. -1)GO TO 340
270  iwm(lnje)=iwm(lnje)+1
271  jcalc=0
272  CALL ddajac(neq,x,y,yprime,delta,cj,h,
273  * ier,wt,e,wm,iwm,res,ires,uround,jac,rpar,
274  * ipar,ntemp)
275  cjold=cj
276  s = 100.d0
277  IF (ires .LT. 0) GO TO 380
278  IF(ier .NE. 0)GO TO 380
279  nsf=0
280 C
281 C
282 C INITIALIZE THE ERROR ACCUMULATION VECTOR E.
283 340 CONTINUE
284  DO 345 i=1,neq
285 345 e(i)=0.0d0
286 C
287 C
288 C CORRECTOR LOOP.
289 350 CONTINUE
290 C
291 C MULTIPLY RESIDUAL BY TEMP1 TO ACCELERATE CONVERGENCE
292  temp1 = 2.0d0/(1.0d0 + cj/cjold)
293  DO 355 i = 1,neq
294 355 delta(i) = delta(i) * temp1
295 C
296 C COMPUTE A NEW ITERATE (BACK-SUBSTITUTION).
297 C STORE THE CORRECTION IN DELTA.
298  CALL ddaslv(neq,delta,wm,iwm)
299 C
300 C UPDATE Y,E,AND YPRIME
301  DO 360 i=1,neq
302  y(i)=y(i)-delta(i)
303  e(i)=e(i)-delta(i)
304 360 yprime(i)=yprime(i)-cj*delta(i)
305 C
306 C TEST FOR CONVERGENCE OF THE ITERATION
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
310  oldnrm = delnrm
311  GO TO 367
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
316 C
317 C THE CORRECTOR HAS NOT YET CONVERGED.
318 C UPDATE M AND TEST WHETHER THE
319 C MAXIMUM NUMBER OF ITERATIONS HAVE
320 C BEEN TRIED.
321  m=m+1
322  IF(m.GE.maxit)GO TO 370
323 C
324 C EVALUATE THE RESIDUAL
325 C AND GO BACK TO DO ANOTHER ITERATION
326  iwm(lnre)=iwm(lnre)+1
327  ires = 0
328  CALL res(x,y,yprime,delta,ires,
329  * rpar,ipar)
330  IF (ires .LT. 0) GO TO 380
331  GO TO 350
332 C
333 C
334 C THE CORRECTOR FAILED TO CONVERGE IN MAXIT
335 C ITERATIONS. IF THE ITERATION MATRIX
336 C IS NOT CURRENT,RE-DO THE STEP WITH
337 C A NEW ITERATION MATRIX.
338 370 CONTINUE
339  IF(jcalc.EQ.0)GO TO 380
340  jcalc=-1
341  GO TO 300
342 C
343 C
344 C THE ITERATION HAS CONVERGED. IF NONNEGATIVITY OF SOLUTION IS
345 C REQUIRED, SET THE SOLUTION NONNEGATIVE, IF THE PERTURBATION
346 C TO DO IT IS SMALL ENOUGH. IF THE CHANGE IS TOO LARGE, THEN
347 C CONSIDER THE CORRECTOR ITERATION TO HAVE FAILED.
348 375 IF(nonneg .EQ. 0) GO TO 390
349  DO 377 i = 1,neq
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
353  DO 378 i = 1,neq
354 378 e(i) = e(i) - delta(i)
355  GO TO 390
356 C
357 C
358 C EXITS FROM BLOCK 3
359 C NO CONVERGENCE WITH CURRENT ITERATION
360 C MATRIX,OR SINGULAR ITERATION MATRIX
361 380 convgd= .false.
362 390 jcalc = 1
363  IF(.NOT.convgd)GO TO 600
364 C
365 C
366 C
367 C
368 C
369 C-----------------------------------------------------------------------
370 C BLOCK 4
371 C ESTIMATE THE ERRORS AT ORDERS K,K-1,K-2
372 C AS IF CONSTANT STEPSIZE WAS USED. ESTIMATE
373 C THE LOCAL ERROR AT ORDER K AND TEST
374 C WHETHER THE CURRENT STEP IS SUCCESSFUL.
375 C-----------------------------------------------------------------------
376 C
377 C ESTIMATE ERRORS AT ORDERS K,K-1,K-2
378  enorm = ddanrm(neq,e,wt,rpar,ipar)
379  erk = sigma(k+1)*enorm
380  terk = (k+1)*erk
381  est = erk
382  knew=k
383  IF(k .EQ. 1)GO TO 430
384  DO 405 i = 1,neq
385 405 delta(i) = phi(i,kp1) + e(i)
386  erkm1=sigma(k)*ddanrm(neq,delta,wt,rpar,ipar)
387  terkm1 = k*erkm1
388  IF(k .GT. 2)GO TO 410
389  IF(terkm1 .LE. 0.5d0*terk)GO TO 420
390  GO TO 430
391 410 CONTINUE
392  DO 415 i = 1,neq
393 415 delta(i) = phi(i,k) + delta(i)
394  erkm2=sigma(k-1)*ddanrm(neq,delta,wt,rpar,ipar)
395  terkm2 = (k-1)*erkm2
396  IF(max(terkm1,terkm2).GT.terk)GO TO 430
397 C LOWER THE ORDER
398 420 CONTINUE
399  knew=k-1
400  est = erkm1
401 C
402 C
403 C CALCULATE THE LOCAL ERROR FOR THE CURRENT STEP
404 C TO SEE IF THE STEP WAS SUCCESSFUL
405 430 CONTINUE
406  err = ck * enorm
407  IF(err .GT. 1.0d0)GO TO 600
408 C
409 C
410 C
411 C
412 C
413 C-----------------------------------------------------------------------
414 C BLOCK 5
415 C THE STEP IS SUCCESSFUL. DETERMINE
416 C THE BEST ORDER AND STEPSIZE FOR
417 C THE NEXT STEP. UPDATE THE DIFFERENCES
418 C FOR THE NEXT STEP.
419 C-----------------------------------------------------------------------
420  idid=1
421  iwm(lnst)=iwm(lnst)+1
422  kdiff=k-kold
423  kold=k
424  hold=h
425 C
426 C
427 C ESTIMATE THE ERROR AT ORDER K+1 UNLESS:
428 C ALREADY DECIDED TO LOWER ORDER, OR
429 C ALREADY USING MAXIMUM ORDER, OR
430 C STEPSIZE NOT CONSTANT, OR
431 C ORDER RAISED IN PREVIOUS STEP
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
437  DO 510 i=1,neq
438 510 delta(i)=e(i)-phi(i,kp2)
439  erkp1 = (1.0d0/(k+2))*ddanrm(neq,delta,wt,rpar,ipar)
440  terkp1 = (k+2)*erkp1
441  IF(k.GT.1)GO TO 520
442  IF(terkp1.GE.0.5d0*terk)GO TO 550
443  GO TO 530
444 520 IF(terkm1.LE.min(terk,terkp1))GO TO 540
445  IF(terkp1.GE.terk.OR.k.EQ.iwm(lmxord))GO TO 550
446 C
447 C RAISE ORDER
448 530 k=kp1
449  est = erkp1
450  GO TO 550
451 C
452 C LOWER ORDER
453 540 k=km1
454  est = erkm1
455  GO TO 550
456 C
457 C IF IPHASE = 0, INCREASE ORDER BY ONE AND MULTIPLY STEPSIZE BY
458 C FACTOR TWO
459 545 k = kp1
460  hnew = h*2.0d0
461  h = hnew
462  GO TO 575
463 C
464 C
465 C DETERMINE THE APPROPRIATE STEPSIZE FOR
466 C THE NEXT STEP.
467 550 hnew=h
468  temp2=k+1
469  r=(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
470  IF(r .LT. 2.0d0) GO TO 555
471  hnew = 2.0d0*h
472  GO TO 560
473 555 IF(r .GT. 1.0d0) GO TO 560
474  r = max(0.5d0,min(0.9d0,r))
475  hnew = h*r
476 560 h=hnew
477 C
478 C
479 C UPDATE DIFFERENCES FOR NEXT STEP
480 575 CONTINUE
481  IF(kold.EQ.iwm(lmxord))GO TO 585
482  DO 580 i=1,neq
483 580 phi(i,kp2)=e(i)
484 585 CONTINUE
485  DO 590 i=1,neq
486 590 phi(i,kp1)=phi(i,kp1)+e(i)
487  DO 595 j1=2,kp1
488  j=kp1-j1+1
489  DO 595 i=1,neq
490 595 phi(i,j)=phi(i,j)+phi(i,j+1)
491  RETURN
492 C
493 C
494 C
495 C
496 C
497 C-----------------------------------------------------------------------
498 C BLOCK 6
499 C THE STEP IS UNSUCCESSFUL. RESTORE X,PSI,PHI
500 C DETERMINE APPROPRIATE STEPSIZE FOR
501 C CONTINUING THE INTEGRATION, OR EXIT WITH
502 C AN ERROR FLAG IF THERE HAVE BEEN MANY
503 C FAILURES.
504 C-----------------------------------------------------------------------
505 600 iphase = 1
506 C
507 C RESTORE X,PHI,PSI
508  x=xold
509  IF(kp1.LT.nsp1)GO TO 630
510  DO 620 j=nsp1,kp1
511  temp1=1.0d0/beta(j)
512  DO 610 i=1,neq
513 610 phi(i,j)=temp1*phi(i,j)
514 620 CONTINUE
515 630 CONTINUE
516  DO 640 i=2,kp1
517 640 psi(i-1)=psi(i)-h
518 C
519 C
520 C TEST WHETHER FAILURE IS DUE TO CORRECTOR ITERATION
521 C OR ERROR TEST
522  IF(convgd)GO TO 660
523  iwm(lctf)=iwm(lctf)+1
524 C
525 C
526 C THE NEWTON ITERATION FAILED TO CONVERGE WITH
527 C A CURRENT ITERATION MATRIX. DETERMINE THE CAUSE
528 C OF THE FAILURE AND TAKE APPROPRIATE ACTION.
529  IF(ier.EQ.0)GO TO 650
530 C
531 C THE ITERATION MATRIX IS SINGULAR. REDUCE
532 C THE STEPSIZE BY A FACTOR OF 4. IF
533 C THIS HAPPENS THREE TIMES IN A ROW ON
534 C THE SAME STEP, RETURN WITH AN ERROR FLAG
535  nsf=nsf+1
536  r = 0.25d0
537  h=h*r
538  IF (nsf .LT. 3 .AND. abs(h) .GE. hmin) GO TO 690
539  idid=-8
540  GO TO 675
541 C
542 C
543 C THE NEWTON ITERATION FAILED TO CONVERGE FOR A REASON
544 C OTHER THAN A SINGULAR ITERATION MATRIX. IF IRES = -2, THEN
545 C RETURN. OTHERWISE, REDUCE THE STEPSIZE AND TRY AGAIN, UNLESS
546 C TOO MANY FAILURES HAVE OCCURRED.
547 650 CONTINUE
548  IF (ires .GT. -2) GO TO 655
549  idid = -11
550  GO TO 675
551 655 ncf = ncf + 1
552  r = 0.25d0
553  h = h*r
554  IF (ncf .LT. 10 .AND. abs(h) .GE. hmin) GO TO 690
555  idid = -7
556  IF (ires .LT. 0) idid = -10
557  IF (nef .GE. 3) idid = -9
558  GO TO 675
559 C
560 C
561 C THE NEWTON SCHEME CONVERGED,AND THE CAUSE
562 C OF THE FAILURE WAS THE ERROR ESTIMATE
563 C EXCEEDING THE TOLERANCE.
564 660 nef=nef+1
565  iwm(letf)=iwm(letf)+1
566  IF (nef .GT. 1) GO TO 665
567 C
568 C ON FIRST ERROR TEST FAILURE, KEEP CURRENT ORDER OR LOWER
569 C ORDER BY ONE. COMPUTE NEW STEPSIZE BASED ON DIFFERENCES
570 C OF THE SOLUTION.
571  k = knew
572  temp2 = k + 1
573  r = 0.90d0*(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
574  r = max(0.25d0,min(0.9d0,r))
575  h = h*r
576  IF (abs(h) .GE. hmin) GO TO 690
577  idid = -6
578  GO TO 675
579 C
580 C ON SECOND ERROR TEST FAILURE, USE THE CURRENT ORDER OR
581 C DECREASE ORDER BY ONE. REDUCE THE STEPSIZE BY A FACTOR OF
582 C FOUR.
583 665 IF (nef .GT. 2) GO TO 670
584  k = knew
585  h = 0.25d0*h
586  IF (abs(h) .GE. hmin) GO TO 690
587  idid = -6
588  GO TO 675
589 C
590 C ON THIRD AND SUBSEQUENT ERROR TEST FAILURES, SET THE ORDER TO
591 C ONE AND REDUCE THE STEPSIZE BY A FACTOR OF FOUR.
592 670 k = 1
593  h = 0.25d0*h
594  IF (abs(h) .GE. hmin) GO TO 690
595  idid = -6
596  GO TO 675
597 C
598 C
599 C
600 C
601 C FOR ALL CRASHES, RESTORE Y TO ITS LAST VALUE,
602 C INTERPOLATE TO FIND YPRIME AT LAST X, AND RETURN
603 675 CONTINUE
604  CALL ddatrp(x,x,y,yprime,neq,k,phi,psi)
605  RETURN
606 C
607 C
608 C GO BACK AND TRY THIS STEP AGAIN
609 690 GO TO 200
610 C
611 C------END OF SUBROUTINE DDASTP------
612  END
static T abs(T x)
Definition: pr-output.cc:1696
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:227
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:204