-
Notifications
You must be signed in to change notification settings - Fork 0
/
opkda1.f90
9795 lines (9732 loc) · 343 KB
/
opkda1.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
DOUBLE PRECISION FUNCTION DUMACH ()
!***BEGIN PROLOGUE DUMACH
!***PURPOSE Compute the unit roundoff of the machine.
!***CATEGORY R1
!***TYPE DOUBLE PRECISION (RUMACH-S, DUMACH-D)
!***KEYWORDS MACHINE CONSTANTS
!***AUTHOR Hindmarsh, Alan C., (LLNL)
!***DESCRIPTION
! *Usage:
! DOUBLE PRECISION A, DUMACH
! A = DUMACH()
!
! *Function Return Values:
! A : the unit roundoff of the machine.
!
! *Description:
! The unit roundoff is defined as the smallest positive machine
! number u such that 1.0 + u .ne. 1.0. This is computed by DUMACH
! in a machine-independent manner.
!
!***REFERENCES (NONE)
!***ROUTINES CALLED DUMSUM
!***REVISION HISTORY (YYYYMMDD)
! 19930216 DATE WRITTEN
! 19930818 Added SLATEC-format prologue. (FNF)
! 20030707 Added DUMSUM to force normal storage of COMP. (ACH)
!***END PROLOGUE DUMACH
!
DOUBLE PRECISION U, COMP
!***FIRST EXECUTABLE STATEMENT DUMACH
U = 1.0D0
10 U = U*0.5D0
CALL DUMSUM(1.0D0, U, COMP)
IF (COMP .NE. 1.0D0) GOTO 10
DUMACH = U*2.0D0
RETURN
!----------------------- End of Function DUMACH ------------------------
END
SUBROUTINE DUMSUM(A,B,C)
! Routine to force normal storing of A + B, for DUMACH.
DOUBLE PRECISION A, B, C
C = A + B
RETURN
END
SUBROUTINE DCFODE (METH, ELCO, TESCO)
!***BEGIN PROLOGUE DCFODE
!***SUBSIDIARY
!***PURPOSE Set ODE integrator coefficients.
!***TYPE DOUBLE PRECISION (SCFODE-S, DCFODE-D)
!***AUTHOR Hindmarsh, Alan C., (LLNL)
!***DESCRIPTION
!
! DCFODE is called by the integrator routine to set coefficients
! needed there. The coefficients for the current method, as
! given by the value of METH, are set for all orders and saved.
! The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2.
! (A smaller value of the maximum order is also allowed.)
! DCFODE is called once at the beginning of the problem,
! and is not called again unless and until METH is changed.
!
! The ELCO array contains the basic method coefficients.
! The coefficients el(i), 1 .le. i .le. nq+1, for the method of
! order nq are stored in ELCO(i,nq). They are given by a genetrating
! polynomial, i.e.,
! l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq.
! For the implicit Adams methods, l(x) is given by
! dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0.
! For the BDF methods, l(x) is given by
! l(x) = (x+1)*(x+2)* ... *(x+nq)/K,
! where K = factorial(nq)*(1 + 1/2 + ... + 1/nq).
!
! The TESCO array contains test constants used for the
! local error test and the selection of step size and/or order.
! At order nq, TESCO(k,nq) is used for the selection of step
! size at order nq - 1 if k = 1, at order nq if k = 2, and at order
! nq + 1 if k = 3.
!
!***SEE ALSO DLSODE
!***ROUTINES CALLED (NONE)
!***REVISION HISTORY (YYMMDD)
! 791129 DATE WRITTEN
! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
! 890503 Minor cosmetic changes. (FNF)
! 930809 Renamed to allow single/double precision versions. (ACH)
!***END PROLOGUE DCFODE
!**End
INTEGER METH
INTEGER I, IB, NQ, NQM1, NQP1
DOUBLE PRECISION ELCO, TESCO
DOUBLE PRECISION AGAMQ, FNQ, FNQM1, PC, PINT, RAGQ,RQFAC, RQ1FAC, TSIGN, XPIN
DIMENSION ELCO(13,12), TESCO(3,12)
DIMENSION PC(12)
!
!***FIRST EXECUTABLE STATEMENT DCFODE
GOTO (100, 200), METH
!
100 ELCO(1,1) = 1.0D0
ELCO(2,1) = 1.0D0
TESCO(1,1) = 0.0D0
TESCO(2,1) = 2.0D0
TESCO(1,2) = 1.0D0
TESCO(3,12) = 0.0D0
PC(1) = 1.0D0
RQFAC = 1.0D0
DO 140 NQ = 2,12
!-----------------------------------------------------------------------
! The PC array will contain the coefficients of the polynomial
! p(x) = (x+1)*(x+2)*...*(x+nq-1).
! Initially, p(x) = 1.
!-----------------------------------------------------------------------
RQ1FAC = RQFAC
RQFAC = RQFAC/NQ
NQM1 = NQ - 1
FNQM1 = NQM1
NQP1 = NQ + 1
! Form coefficients of p(x)*(x+nq-1). ----------------------------------
PC(NQ) = 0.0D0
DO 110 IB = 1,NQM1
I = NQP1 - IB
110 PC(I) = PC(I-1) + FNQM1*PC(I)
PC(1) = FNQM1*PC(1)
! Compute integral, -1 to 0, of p(x) and x*p(x). -----------------------
PINT = PC(1)
XPIN = PC(1)/2.0D0
TSIGN = 1.0D0
DO 120 I = 2,NQ
TSIGN = -TSIGN
PINT = PINT + TSIGN*PC(I)/I
120 XPIN = XPIN + TSIGN*PC(I)/(I+1)
! Store coefficients in ELCO and TESCO. --------------------------------
ELCO(1,NQ) = PINT*RQ1FAC
ELCO(2,NQ) = 1.0D0
DO 130 I = 2,NQ
130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/I
AGAMQ = RQFAC*XPIN
RAGQ = 1.0D0/AGAMQ
TESCO(2,NQ) = RAGQ
IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/NQP1
TESCO(3,NQM1) = RAGQ
140 CONTINUE
RETURN
!
200 PC(1) = 1.0D0
RQ1FAC = 1.0D0
DO 230 NQ = 1,5
!-----------------------------------------------------------------------
! The PC array will contain the coefficients of the polynomial
! p(x) = (x+1)*(x+2)*...*(x+nq).
! Initially, p(x) = 1.
!-----------------------------------------------------------------------
FNQ = NQ
NQP1 = NQ + 1
! Form coefficients of p(x)*(x+nq). ------------------------------------
PC(NQP1) = 0.0D0
DO 210 IB = 1,NQ
I = NQ + 2 - IB
210 PC(I) = PC(I-1) + FNQ*PC(I)
PC(1) = FNQ*PC(1)
! Store coefficients in ELCO and TESCO. --------------------------------
DO 220 I = 1,NQP1
220 ELCO(I,NQ) = PC(I)/PC(2)
ELCO(2,NQ) = 1.0D0
TESCO(1,NQ) = RQ1FAC
TESCO(2,NQ) = NQP1/ELCO(1,NQ)
TESCO(3,NQ) = (NQ+2)/ELCO(1,NQ)
RQ1FAC = RQ1FAC/FNQ
230 CONTINUE
RETURN
!----------------------- END OF SUBROUTINE DCFODE ----------------------
END
SUBROUTINE DINTDY (T, K, YH, NYH, DKY, IFLAG)
!***BEGIN PROLOGUE DINTDY
!***SUBSIDIARY
!***PURPOSE Interpolate solution derivatives.
!***TYPE DOUBLE PRECISION (SINTDY-S, DINTDY-D)
!***AUTHOR Hindmarsh, Alan C., (LLNL)
!***DESCRIPTION
!
! DINTDY computes interpolated values of the K-th derivative of the
! dependent variable vector y, and stores it in DKY. This routine
! is called within the package with K = 0 and T = TOUT, but may
! also be called by the user for any K up to the current order.
! (See detailed instructions in the usage documentation.)
!
! The computed values in DKY are gotten by interpolation using the
! Nordsieck history array YH. This array corresponds uniquely to a
! vector-valued polynomial of degree NQCUR or less, and DKY is set
! to the K-th derivative of this polynomial at T.
! The formula for DKY is:
! q
! DKY(i) = sum c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1)
! j=K
! where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR.
! The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are
! communicated by COMMON. The above sum is done in reverse order.
! IFLAG is returned negative if either K or T is out of bounds.
!
!***SEE ALSO DLSODE
!***ROUTINES CALLED XERRWD
!***COMMON BLOCKS DLS001
!***REVISION HISTORY (YYMMDD)
! 791129 DATE WRITTEN
! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
! 890503 Minor cosmetic changes. (FNF)
! 930809 Renamed to allow single/double precision versions. (ACH)
! 010418 Reduced size of Common block /DLS001/. (ACH)
! 031105 Restored 'own' variables to Common block /DLS001/, to
! enable interrupt/restart feature. (ACH)
! 050427 Corrected roundoff decrement in TP. (ACH)
!***END PROLOGUE DINTDY
!**End
INTEGER K, NYH, IFLAG
DOUBLE PRECISION T, YH, DKY
DIMENSION YH(NYH,*), DKY(*)
INTEGER IOWND, IOWNS, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, LYH, LEWT, LACOR, LSAVF
INTEGER LWM, LIWM, METH, MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
DOUBLE PRECISION ROWNS,CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
COMMON /DLS001/ ROWNS(209), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, IOWND(6), IOWNS(6),&
ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,&
MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
DOUBLE PRECISION C, R, S, TP
CHARACTER*80 MSG
!
!***FIRST EXECUTABLE STATEMENT DINTDY
IFLAG = 0
IF (K .LT. 0 .OR. K .GT. NQ) GOTO 80
TP = TN - HU - 100.0D0*UROUND*SIGN(ABS(TN) + ABS(HU), HU)
IF ((T-TP)*(T-TN) .GT. 0.0D0) GOTO 90
!
S = (T - TN)/H
IC = 1
IF (K .EQ. 0) GOTO 15
JJ1 = L - K
DO 10 JJ = JJ1,NQ
10 IC = IC*JJ
15 C = IC
DO 20 I = 1,N
20 DKY(I) = C*YH(I,L)
IF (K .EQ. NQ) GOTO 55
JB2 = NQ - K
DO 50 JB = 1,JB2
J = NQ - JB
JP1 = J + 1
IC = 1
IF (K .EQ. 0) GOTO 35
JJ1 = JP1 - K
DO 30 JJ = JJ1,J
30 IC = IC*JJ
35 C = IC
DO 40 I = 1,N
40 DKY(I) = C*YH(I,JP1) + S*DKY(I)
50 CONTINUE
IF (K .EQ. 0) RETURN
55 R = H**(-K)
DO 60 I = 1,N
60 DKY(I) = R*DKY(I)
RETURN
!
80 MSG = 'DINTDY- K (=I1) illegal '
CALL XERRWD (MSG, 30, 51, 0, 1, K, 0, 0, 0.0D0, 0.0D0)
IFLAG = -1
RETURN
90 MSG = 'DINTDY- T (=R1) illegal '
CALL XERRWD (MSG, 30, 52, 0, 0, 0, 0, 1, T, 0.0D0)
MSG=' T not in interval TCUR - HU (= R1) to TCUR (=R2) '
CALL XERRWD (MSG, 60, 52, 0, 0, 0, 0, 2, TP, TN)
IFLAG = -2
RETURN
!----------------------- END OF SUBROUTINE DINTDY ----------------------
END
SUBROUTINE DPREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC)
!***BEGIN PROLOGUE DPREPJ
!***SUBSIDIARY
!***PURPOSE Compute and process Newton iteration matrix.
!***TYPE DOUBLE PRECISION (SPREPJ-S, DPREPJ-D)
!***AUTHOR Hindmarsh, Alan C., (LLNL)
!***DESCRIPTION
!
! DPREPJ is called by DSTODE to compute and process the matrix
! P = I - h*el(1)*J , where J is an approximation to the Jacobian.
! Here J is computed by the user-supplied routine JAC if
! MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5.
! If MITER = 3, a diagonal approximation to J is used.
! J is stored in WM and replaced by P. If MITER .ne. 3, P is then
! subjected to LU decomposition in preparation for later solution
! of linear systems with P as coefficient matrix. This is done
! by DGEFA if MITER = 1 or 2, and by DGBFA if MITER = 4 or 5.
!
! In addition to variables described in DSTODE and DLSODE prologues,
! communication with DPREPJ uses the following:
! Y = array containing predicted values on entry.
! FTEM = work array of length N (ACOR in DSTODE).
! SAVF = array containing f evaluated at predicted y.
! WM = real work space for matrices. On output it contains the
! inverse diagonal matrix if MITER = 3 and the LU decomposition
! of P if MITER is 1, 2 , 4, or 5.
! Storage of matrix elements starts at WM(3).
! WM also contains the following matrix-related data:
! WM(1) = SQRT(UROUND), used in numerical Jacobian increments.
! WM(2) = H*EL0, saved for later use if MITER = 3.
! IWM = integer work space containing pivot information, starting at
! IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band
! parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
! EL0 = EL(1) (input).
! IERPJ = output error flag, = 0 if no trouble, .gt. 0 if
! P matrix found to be singular.
! JCUR = output flag = 1 to indicate that the Jacobian matrix
! (or approximation) is now current.
! This routine also uses the COMMON variables EL0, H, TN, UROUND,
! MITER, N, NFE, and NJE.
!
!***SEE ALSO DLSODE
!***ROUTINES CALLED DGBFA, DGEFA, DVNORM
!***COMMON BLOCKS DLS001
!***REVISION HISTORY (YYMMDD)
! 791129 DATE WRITTEN
! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
! 890504 Minor cosmetic changes. (FNF)
! 930809 Renamed to allow single/double precision versions. (ACH)
! 010418 Reduced size of Common block /DLS001/. (ACH)
! 031105 Restored 'own' variables to Common block /DLS001/, to
! enable interrupt/restart feature. (ACH)
!***END PROLOGUE DPREPJ
!**End
EXTERNAL F, JAC
INTEGER NEQ, NYH, IWM
DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WM
DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*), WM(*), IWM(*)
INTEGER IOWND, IOWNS, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, LYH, LEWT, LACOR
INTEGER LSAVF, LWM, LIWM, METH, MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
DOUBLE PRECISION ROWNS, CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
COMMON /DLS001/ ROWNS(209), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, IOWND(6), IOWNS(6),&
ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,&
MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP, MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1
DOUBLE PRECISION CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ, DVNORM
!
!***FIRST EXECUTABLE STATEMENT DPREPJ
NJE = NJE + 1
IERPJ = 0
JCUR = 1
HL0 = H*EL0
GOTO (100, 200, 300, 400, 500), MITER
! If MITER = 1, call JAC and multiply by scalar. -----------------------
100 LENP = N*N
DO 110 I = 1,LENP
110 WM(I+2) = 0.0D0
CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N)
CON = -HL0
DO 120 I = 1,LENP
120 WM(I+2) = WM(I+2)*CON
GOTO 240
! If MITER = 2, make N calls to F to approximate J. --------------------
200 FAC = DVNORM (N, SAVF, EWT)
R0 = 1000.0D0*ABS(H)*UROUND*N*FAC
IF (R0 .EQ. 0.0D0) R0 = 1.0D0
SRUR = WM(1)
J1 = 2
DO 230 J = 1,N
YJ = Y(J)
R = MAX(SRUR*ABS(YJ),R0/EWT(J))
Y(J) = Y(J) + R
FAC = -HL0/R
CALL F (NEQ, TN, Y, FTEM)
DO 220 I = 1,N
220 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
Y(J) = YJ
J1 = J1 + N
230 CONTINUE
NFE = NFE + N
! Add identity matrix. -------------------------------------------------
240 J = 3
NP1 = N + 1
DO 250 I = 1,N
WM(J) = WM(J) + 1.0D0
250 J = J + NP1
! Do LU decomposition on P. --------------------------------------------
CALL DGEFA (WM(3), N, N, IWM(21), IER)
IF (IER .NE. 0) IERPJ = 1
RETURN
! If MITER = 3, construct a diagonal approximation to J and P. ---------
300 WM(2) = HL0
R = EL0*0.1D0
DO 310 I = 1,N
310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
CALL F (NEQ, TN, Y, WM(3))
NFE = NFE + 1
DO 320 I = 1,N
R0 = H*SAVF(I) - YH(I,2)
DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I))
WM(I+2) = 1.0D0
IF (ABS(R0) .LT. UROUND/EWT(I)) GOTO 320
IF (ABS(DI) .EQ. 0.0D0) GOTO 330
WM(I+2) = 0.1D0*R0/DI
320 CONTINUE
RETURN
330 IERPJ = 1
RETURN
! If MITER = 4, call JAC and multiply by scalar. -----------------------
400 ML = IWM(1)
MU = IWM(2)
ML3 = ML + 3
MBAND = ML + MU + 1
MEBAND = MBAND + ML
LENP = MEBAND*N
DO 410 I = 1,LENP
410 WM(I+2) = 0.0D0
CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND)
CON = -HL0
DO 420 I = 1,LENP
420 WM(I+2) = WM(I+2)*CON
GOTO 570
! If MITER = 5, make MBAND calls to F to approximate J. ----------------
500 ML = IWM(1)
MU = IWM(2)
MBAND = ML + MU + 1
MBA = MIN(MBAND,N)
MEBAND = MBAND + ML
MEB1 = MEBAND - 1
SRUR = WM(1)
FAC = DVNORM (N, SAVF, EWT)
R0 = 1000.0D0*ABS(H)*UROUND*N*FAC
IF (R0 .EQ. 0.0D0) R0 = 1.0D0
DO 560 J = 1,MBA
DO 530 I = J,N,MBAND
YI = Y(I)
R = MAX(SRUR*ABS(YI),R0/EWT(I))
530 Y(I) = Y(I) + R
CALL F (NEQ, TN, Y, FTEM)
DO 550 JJ = J,N,MBAND
Y(JJ) = YH(JJ,1)
YJJ = Y(JJ)
R = MAX(SRUR*ABS(YJJ),R0/EWT(JJ))
FAC = -HL0/R
I1 = MAX(JJ-MU,1)
I2 = MIN(JJ+ML,N)
II = JJ*MEB1 - ML + 2
DO 540 I = I1,I2
540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC
550 CONTINUE
560 CONTINUE
NFE = NFE + MBA
! Add identity matrix. -------------------------------------------------
570 II = MBAND + 2
DO 580 I = 1,N
WM(II) = WM(II) + 1.0D0
580 II = II + MEBAND
! Do LU decomposition of P. --------------------------------------------
CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
IF (IER .NE. 0) IERPJ = 1
RETURN
!----------------------- END OF SUBROUTINE DPREPJ ----------------------
END
SUBROUTINE DSOLSY (WM, IWM, X, TEM)
!***BEGIN PROLOGUE DSOLSY
!***SUBSIDIARY
!***PURPOSE ODEPACK linear system solver.
!***TYPE DOUBLE PRECISION (SSOLSY-S, DSOLSY-D)
!***AUTHOR Hindmarsh, Alan C., (LLNL)
!***DESCRIPTION
!
! This routine manages the solution of the linear system arising from
! a chord iteration. It is called if MITER .ne. 0.
! If MITER is 1 or 2, it calls DGESL to accomplish this.
! If MITER = 3 it updates the coefficient h*EL0 in the diagonal
! matrix, and then computes the solution.
! If MITER is 4 or 5, it calls DGBSL.
! Communication with DSOLSY uses the following variables:
! WM = real work space containing the inverse diagonal matrix if
! MITER = 3 and the LU decomposition of the matrix otherwise.
! Storage of matrix elements starts at WM(3).
! WM also contains the following matrix-related data:
! WM(1) = SQRT(UROUND) (not used here),
! WM(2) = HL0, the previous value of h*EL0, used if MITER = 3.
! IWM = integer work space containing pivot information, starting at
! IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band
! parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
! X = the right-hand side vector on input, and the solution vector
! on output, of length N.
! TEM = vector of work space of length N, not used in this version.
! IERSL = output flag (in COMMON). IERSL = 0 if no trouble occurred.
! IERSL = 1 if a singular matrix arose with MITER = 3.
! This routine also uses the COMMON variables EL0, H, MITER, and N.
!
!***SEE ALSO DLSODE
!***ROUTINES CALLED DGBSL, DGESL
!***COMMON BLOCKS DLS001
!***REVISION HISTORY (YYMMDD)
! 791129 DATE WRITTEN
! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
! 890503 Minor cosmetic changes. (FNF)
! 930809 Renamed to allow single/double precision versions. (ACH)
! 010418 Reduced size of Common block /DLS001/. (ACH)
! 031105 Restored 'own' variables to Common block /DLS001/, to
! enable interrupt/restart feature. (ACH)
!***END PROLOGUE DSOLSY
!**End
INTEGER IWM
DOUBLE PRECISION WM, X, TEM
DIMENSION WM(*), IWM(*), X(*), TEM(*)
INTEGER IOWND, IOWNS, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, MAXORD
INTEGER MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
DOUBLE PRECISION ROWNS, CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
COMMON /DLS001/ ROWNS(209), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, IOWND(6), IOWNS(6), ICF, IERPJ, IERSL, JCUR, &
JSTART, KFLAG, L, LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
INTEGER I, MEBAND, ML, MU
DOUBLE PRECISION DI, HL0, PHL0, R
!
!***FIRST EXECUTABLE STATEMENT DSOLSY
IERSL = 0
GOTO (100, 100, 300, 400, 400), MITER
100 CALL DGESL (WM(3), N, N, IWM(21), X, 0)
RETURN
!
300 PHL0 = WM(2)
HL0 = H*EL0
WM(2) = HL0
IF (HL0 .EQ. PHL0) GOTO 330
R = HL0/PHL0
DO 320 I = 1,N
DI = 1.0D0 - R*(1.0D0 - 1.0D0/WM(I+2))
IF (ABS(DI) .EQ. 0.0D0) GOTO 390
320 WM(I+2) = 1.0D0/DI
330 DO 340 I = 1,N
340 X(I) = WM(I+2)*X(I)
RETURN
390 IERSL = 1
RETURN
!
400 ML = IWM(1)
MU = IWM(2)
MEBAND = 2*ML + MU + 1
CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0)
RETURN
!----------------------- END OF SUBROUTINE DSOLSY ----------------------
END
SUBROUTINE DSRCOM (RSAV, ISAV, JOB)
!***BEGIN PROLOGUE DSRCOM
!***SUBSIDIARY
!***PURPOSE Save/restore ODEPACK COMMON blocks.
!***TYPE DOUBLE PRECISION (SSRCOM-S, DSRCOM-D)
!***AUTHOR Hindmarsh, Alan C., (LLNL)
!***DESCRIPTION
!
! This routine saves or restores (depending on JOB) the contents of
! the COMMON block DLS001, which is used internally
! by one or more ODEPACK solvers.
!
! RSAV = real array of length 218 or more.
! ISAV = integer array of length 37 or more.
! JOB = flag indicating to save or restore the COMMON blocks:
! JOB = 1 if COMMON is to be saved (written to RSAV/ISAV)
! JOB = 2 if COMMON is to be restored (read from RSAV/ISAV)
! A call with JOB = 2 presumes a prior call with JOB = 1.
!
!***SEE ALSO DLSODE
!***ROUTINES CALLED (NONE)
!***COMMON BLOCKS DLS001
!***REVISION HISTORY (YYMMDD)
! 791129 DATE WRITTEN
! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
! 890503 Minor cosmetic changes. (FNF)
! 921116 Deleted treatment of block /EH0001/. (ACH)
! 930801 Reduced Common block length by 2. (ACH)
! 930809 Renamed to allow single/double precision versions. (ACH)
! 010418 Reduced Common block length by 209+12. (ACH)
! 031105 Restored 'own' variables to Common block /DLS001/, to
! enable interrupt/restart feature. (ACH)
! 031112 Added SAVE statement for data-loaded constants.
!***END PROLOGUE DSRCOM
!**End
INTEGER ISAV, JOB
INTEGER ILS
INTEGER I, LENILS, LENRLS
DOUBLE PRECISION RSAV, RLS
DIMENSION RSAV(*), ISAV(*)
SAVE LENRLS, LENILS
COMMON /DLS001/ RLS(218), ILS(37)
DATA LENRLS/218/, LENILS/37/
!
!***FIRST EXECUTABLE STATEMENT DSRCOM
IF (JOB .EQ. 2) GOTO 100
!
DO 10 I = 1,LENRLS
10 RSAV(I) = RLS(I)
DO 20 I = 1,LENILS
20 ISAV(I) = ILS(I)
RETURN
!
100 CONTINUE
DO 110 I = 1,LENRLS
110 RLS(I) = RSAV(I)
DO 120 I = 1,LENILS
120 ILS(I) = ISAV(I)
RETURN
!----------------------- END OF SUBROUTINE DSRCOM ----------------------
END
SUBROUTINE DSTODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, F, JAC, PJAC, SLVS)
!***BEGIN PROLOGUE DSTODE
!***SUBSIDIARY
!***PURPOSE Performs one step of an ODEPACK integration.
!***TYPE DOUBLE PRECISION (SSTODE-S, DSTODE-D)
!***AUTHOR Hindmarsh, Alan C., (LLNL)
!***DESCRIPTION
!
! DSTODE performs one step of the integration of an initial value
! problem for a system of ordinary differential equations.
! Note: DSTODE is independent of the value of the iteration method
! indicator MITER, when this is .ne. 0, and hence is independent
! of the type of chord method used, or the Jacobian structure.
! Communication with DSTODE is done with the following variables:
!
! NEQ = integer array containing problem size in NEQ(1), and
! passed as the NEQ argument in all calls to F and JAC.
! Y = an array of length .ge. N used as the Y argument in
! all calls to F and JAC.
! YH = an NYH by LMAX array containing the dependent variables
! and their approximate scaled derivatives, where
! LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
! j-th derivative of y(i), scaled by h**j/factorial(j)
! (j = 0,1,...,NQ). on entry for the first step, the first
! two columns of YH must be set from the initial values.
! NYH = a constant integer .ge. N, the first dimension of YH.
! YH1 = a one-dimensional array occupying the same space as YH.
! EWT = an array of length N containing multiplicative weights
! for local error measurements. Local errors in Y(i) are
! compared to 1.0/EWT(i) in various error tests.
! SAVF = an array of working storage, of length N.
! Also used for input of YH(*,MAXORD+2) when JSTART = -1
! and MAXORD .lt. the current order NQ.
! ACOR = a work array of length N, used for the accumulated
! corrections. On a successful return, ACOR(i) contains
! the estimated one-step local error in Y(i).
! WM,IWM = real and integer work arrays associated with matrix
! operations in chord iteration (MITER .ne. 0).
! PJAC = name of routine to evaluate and preprocess Jacobian matrix
! and P = I - h*el0*JAC, if a chord method is being used.
! SLVS = name of routine to solve linear system in chord iteration.
! CCMAX = maximum relative change in h*el0 before PJAC is called.
! H = the step size to be attempted on the next step.
! H is altered by the error control algorithm during the
! problem. H can be either positive or negative, but its
! sign must remain constant throughout the problem.
! HMIN = the minimum absolute value of the step size h to be used.
! HMXI = inverse of the maximum absolute value of h to be used.
! HMXI = 0.0 is allowed and corresponds to an infinite hmax.
! HMIN and HMXI may be changed at any time, but will not
! take effect until the next change of h is considered.
! TN = the independent variable. TN is updated on each step taken.
! JSTART = an integer used for input only, with the following
! values and meanings:
! 0 perform the first step.
! .gt.0 take a new step continuing from the last.
! -1 take the next step with a new value of H, MAXORD,
! N, METH, MITER, and/or matrix parameters.
! -2 take the next step with a new value of H,
! but with other inputs unchanged.
! On return, JSTART is set to 1 to facilitate continuation.
! KFLAG = a completion code with the following meanings:
! 0 the step was succesful.
! -1 the requested error could not be achieved.
! -2 corrector convergence could not be achieved.
! -3 fatal error in PJAC or SLVS.
! A return with KFLAG = -1 or -2 means either
! abs(H) = HMIN or 10 consecutive failures occurred.
! On a return with KFLAG negative, the values of TN and
! the YH array are as of the beginning of the last
! step, and H is the last step size attempted.
! MAXORD = the maximum order of integration method to be allowed.
! MAXCOR = the maximum number of corrector iterations allowed.
! MSBP = maximum number of steps between PJAC calls (MITER .gt. 0).
! MXNCF = maximum number of convergence failures allowed.
! METH/MITER = the method flags. See description in driver.
! N = the number of first-order differential equations.
! The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD,
! MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON.
!
!***SEE ALSO DLSODE
!***ROUTINES CALLED DCFODE, DVNORM
!***COMMON BLOCKS DLS001
!***REVISION HISTORY (YYMMDD)
! 791129 DATE WRITTEN
! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
! 890503 Minor cosmetic changes. (FNF)
! 930809 Renamed to allow single/double precision versions. (ACH)
! 010418 Reduced size of Common block /DLS001/. (ACH)
! 031105 Restored 'own' variables to Common block /DLS001/, to
! enable interrupt/restart feature. (ACH)
!***END PROLOGUE DSTODE
!**End
EXTERNAL F, JAC, PJAC, SLVS
INTEGER NEQ, NYH, IWM
DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM
DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*),ACOR(*), WM(*), IWM(*)
INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L
INTEGER LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP, R, RH, RHDN, RHSM, RHUP, TOLD, DVNORM
COMMON /DLS001/ CONIT, CRATE, EL(13), ELCO(13,12), HOLD, RMAX, TESCO(3,12), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, LYH, LEWT, LACOR, LSAVF, LWM, &
LIWM, METH, MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
!
!***FIRST EXECUTABLE STATEMENT DSTODE
KFLAG = 0
TOLD = TN
NCF = 0
IERPJ = 0
IERSL = 0
JCUR = 0
ICF = 0
DELP = 0.0D0
IF (JSTART .GT. 0) GOTO 200
IF (JSTART .EQ. -1) GOTO 100
IF (JSTART .EQ. -2) GOTO 160
!-----------------------------------------------------------------------
! On the first call, the order is set to 1, and other variables are
! initialized. RMAX is the maximum ratio by which H can be increased
! in a single step. It is initially 1.E4 to compensate for the small
! initial H, but then is normally equal to 10. If a failure
! occurs (in corrector convergence or error test), RMAX is set to 2
! for the next increase.
!-----------------------------------------------------------------------
LMAX = MAXORD + 1
NQ = 1
L = 2
IALTH = 2
RMAX = 10000.0D0
RC = 0.0D0
EL0 = 1.0D0
CRATE = 0.7D0
HOLD = H
MEO = METH
NSLP = 0
IPUP = MITER
IRET = 3
GOTO 140
!-----------------------------------------------------------------------
! The following block handles preliminaries needed when JSTART = -1.
! IPUP is set to MITER to force a matrix update.
! If an order increase is about to be considered (IALTH = 1),
! IALTH is reset to 2 to postpone consideration one more step.
! If the caller has changed METH, DCFODE is called to reset
! the coefficients of the method.
! If the caller has changed MAXORD to a value less than the current
! order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
! If H is to be changed, YH must be rescaled.
! If H or METH is being changed, IALTH is reset to L = NQ + 1
! to prevent further changes in H for that many steps.
!-----------------------------------------------------------------------
100 IPUP = MITER
LMAX = MAXORD + 1
IF (IALTH .EQ. 1) IALTH = 2
IF (METH .EQ. MEO) GOTO 110
CALL DCFODE (METH, ELCO, TESCO)
MEO = METH
IF (NQ .GT. MAXORD) GOTO 120
IALTH = L
IRET = 1
GOTO 150
110 IF (NQ .LE. MAXORD) GOTO 160
120 NQ = MAXORD
L = LMAX
DO 125 I = 1,L
125 EL(I) = ELCO(I,NQ)
NQNYH = NQ*NYH
RC = RC*EL(1)/EL0
EL0 = EL(1)
CONIT = 0.5D0/(NQ+2)
DDN = DVNORM (N, SAVF, EWT)/TESCO(1,L)
EXDN = 1.0D0/L
RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
RH = MIN(RHDN,1.0D0)
IREDO = 3
IF (H .EQ. HOLD) GOTO 170
RH = MIN(RH,ABS(H/HOLD))
H = HOLD
GOTO 175
!-----------------------------------------------------------------------
! DCFODE is called to get all the integration coefficients for the
! current METH. Then the EL vector and related constants are reset
! whenever the order NQ is changed, or at the start of the problem.
!-----------------------------------------------------------------------
140 CALL DCFODE (METH, ELCO, TESCO)
150 DO 155 I = 1,L
155 EL(I) = ELCO(I,NQ)
NQNYH = NQ*NYH
RC = RC*EL(1)/EL0
EL0 = EL(1)
CONIT = 0.5D0/(NQ+2)
GOTO (160, 170, 200), IRET
!-----------------------------------------------------------------------
! If H is being changed, the H ratio RH is checked against
! RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
! L = NQ + 1 to prevent a change of H for that many steps, unless
! forced by a convergence or error test failure.
!-----------------------------------------------------------------------
160 IF (H .EQ. HOLD) GOTO 200
RH = H/HOLD
H = HOLD
IREDO = 3
GOTO 175
170 RH = MAX(RH,HMIN/ABS(H))
175 RH = MIN(RH,RMAX)
RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
R = 1.0D0
DO 180 J = 2,L
R = R*RH
DO 180 I = 1,N
180 YH(I,J) = YH(I,J)*R
H = H*RH
RC = RC*RH
IALTH = L
IF (IREDO .EQ. 0) GOTO 690
!-----------------------------------------------------------------------
! This section computes the predicted values by effectively
! multiplying the YH array by the Pascal Triangle matrix.
! RC is the ratio of new to old values of the coefficient H*EL(1).
! When RC differs from 1 by more than CCMAX, IPUP is set to MITER
! to force PJAC to be called, if a Jacobian is involved.
! In any case, PJAC is called at least every MSBP steps.
!-----------------------------------------------------------------------
200 IF (ABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER
IF (NST .GE. NSLP+MSBP) IPUP = MITER
TN = TN + H
I1 = NQNYH + 1
DO 215 JB = 1,NQ
I1 = I1 - NYH
!dir$ ivdep
DO 210 I = I1,NQNYH
210 YH1(I) = YH1(I) + YH1(I+NYH)
215 CONTINUE
!-----------------------------------------------------------------------
! Up to MAXCOR corrector iterations are taken. A convergence test is
! made on the R.M.S. norm of each correction, weighted by the error
! weight vector EWT. The sum of the corrections is accumulated in the
! vector ACOR(i). The YH array is not altered in the corrector loop.
!-----------------------------------------------------------------------
220 M = 0
DO 230 I = 1,N
230 Y(I) = YH(I,1)
CALL F (NEQ, TN, Y, SAVF)
NFE = NFE + 1
IF (IPUP .LE. 0) GOTO 250
!-----------------------------------------------------------------------
! If indicated, the matrix P = I - h*el(1)*J is reevaluated and
! preprocessed before starting the corrector iteration. IPUP is set
! to 0 as an indicator that this has been done.
!-----------------------------------------------------------------------
CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC)
IPUP = 0
RC = 1.0D0
NSLP = NST
CRATE = 0.7D0
IF (IERPJ .NE. 0) GOTO 430
250 DO 260 I = 1,N
260 ACOR(I) = 0.0D0
270 IF (MITER .NE. 0) GOTO 350
!-----------------------------------------------------------------------
! In the case of functional iteration, update Y directly from
! the result of the last function evaluation.
!-----------------------------------------------------------------------
DO 290 I = 1,N
SAVF(I) = H*SAVF(I) - YH(I,2)
290 Y(I) = SAVF(I) - ACOR(I)
DEL = DVNORM (N, Y, EWT)
DO 300 I = 1,N
Y(I) = YH(I,1) + EL(1)*SAVF(I)
300 ACOR(I) = SAVF(I)
GOTO 400
!-----------------------------------------------------------------------
! In the case of the chord method, compute the corrector error,
! and solve the linear system with that as right-hand side and
! P as coefficient matrix.
!-----------------------------------------------------------------------
350 DO 360 I = 1,N
360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
CALL SLVS (WM, IWM, Y, SAVF)
IF (IERSL .LT. 0) GOTO 430
IF (IERSL .GT. 0) GOTO 410
DEL = DVNORM (N, Y, EWT)
DO 380 I = 1,N
ACOR(I) = ACOR(I) + Y(I)
380 Y(I) = YH(I,1) + EL(1)*ACOR(I)
!-----------------------------------------------------------------------
! Test for convergence. If M.gt.0, an estimate of the convergence
! rate constant is stored in CRATE, and this is used in the test.
!-----------------------------------------------------------------------
400 IF (M .NE. 0) CRATE = MAX(0.2D0*CRATE,DEL/DELP)
DCON = DEL*MIN(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT)
IF (DCON .LE. 1.0D0) GOTO 450
M = M + 1
IF (M .EQ. MAXCOR) GOTO 410
IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GOTO 410
DELP = DEL
CALL F (NEQ, TN, Y, SAVF)
NFE = NFE + 1
GOTO 270
!-----------------------------------------------------------------------
! The corrector iteration failed to converge.
! If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
! the next try. Otherwise the YH array is retracted to its values
! before prediction, and H is reduced, if possible. If H cannot be
! reduced or MXNCF failures have occurred, exit with KFLAG = -2.
!-----------------------------------------------------------------------
410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GOTO 430
ICF = 1
IPUP = MITER
GOTO 220
430 ICF = 2
NCF = NCF + 1
RMAX = 2.0D0
TN = TOLD
I1 = NQNYH + 1
DO 445 JB = 1,NQ
I1 = I1 - NYH
!dir$ ivdep
DO 440 I = I1,NQNYH
440 YH1(I) = YH1(I) - YH1(I+NYH)
445 CONTINUE
IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GOTO 680
IF (ABS(H) .LE. HMIN*1.00001D0) GOTO 670
IF (NCF .EQ. MXNCF) GOTO 670
RH = 0.25D0
IPUP = MITER
IREDO = 1
GOTO 170
!-----------------------------------------------------------------------
! The corrector has converged. JCUR is set to 0
! to signal that the Jacobian involved may need updating later.
! The local error test is made and control passes to statement 500
! if it fails.
!-----------------------------------------------------------------------
450 JCUR = 0
IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
IF (M .GT. 0) DSM = DVNORM (N, ACOR, EWT)/TESCO(2,NQ)
IF (DSM .GT. 1.0D0) GOTO 500
!-----------------------------------------------------------------------
! After a successful step, update the YH array.
! Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1.
! If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
! use in a possible order increase on the next step.
! If a change in H is considered, an increase or decrease in order
! by one is considered also. A change in H is made only if it is by a
! factor of at least 1.1. If not, IALTH is set to 3 to prevent
! testing for that many steps.
!-----------------------------------------------------------------------
KFLAG = 0
IREDO = 0
NST = NST + 1
HU = H
NQU = NQ
DO 470 J = 1,L
DO 470 I = 1,N
470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
IALTH = IALTH - 1
IF (IALTH .EQ. 0) GOTO 520
IF (IALTH .GT. 1) GOTO 700
IF (L .EQ. LMAX) GOTO 700
DO 490 I = 1,N
490 YH(I,LMAX) = ACOR(I)
GOTO 700
!-----------------------------------------------------------------------
! The error test failed. KFLAG keeps track of multiple failures.
! Restore TN and the YH array to their previous values, and prepare
! to try the step again. Compute the optimum step size for this or
! one lower order. After 2 or more failures, H is forced to decrease
! by a factor of 0.2 or less.
!-----------------------------------------------------------------------
500 KFLAG = KFLAG - 1
TN = TOLD
I1 = NQNYH + 1
DO 515 JB = 1,NQ
I1 = I1 - NYH
!dir$ ivdep
DO 510 I = I1,NQNYH
510 YH1(I) = YH1(I) - YH1(I+NYH)
515 CONTINUE
RMAX = 2.0D0
IF (ABS(H) .LE. HMIN*1.00001D0) GOTO 660
IF (KFLAG .LE. -3) GOTO 640
IREDO = 2
RHUP = 0.0D0
GOTO 540
!-----------------------------------------------------------------------
! Regardless of the success or failure of the step, factors
! RHDN, RHSM, and RHUP are computed, by which H could be multiplied
! at order NQ - 1, order NQ, or order NQ + 1, respectively.
! In the case of failure, RHUP = 0.0 to avoid an order increase.
! The largest of these is determined and the new order chosen