← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/sdastp.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

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