GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 OCCURED.
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
Generated on Mon Dec 30 2013 03:04:46 for GNU Octave by   doxygen 1.8.1.2