Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/DVODE/stiffoptions.f90 @ 3186

Last change on this file since 3186 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk
File size: 68.9 KB
Line 
1! DVODE_F90 demonstration program
2! Enforces nonnegativity on the Class F problems
3
4! Toronto Stiff Test Set
5
6! This program solves each of the thirty problems for several error
7! tolerances using each of the dense, banded, and sparse solution
8! options and with the numerical Jacobians formed using each of the
9! original VODE algorithm and Doug Salane's JACSP. The ODEs are
10! solved both in scaled and unscaled form. Scalar error tolerances
11! are used for all problems even though this is not appropriate in
12! some cases. The program takes 10-15 minutes to run because several
13! thousand integrations are performed in all. Details of each
14! integration are written to the output file 'stiffoptions.dat'.
15! Summaries of the integration statistics may be found by searching
16! on the string 'Individual'. There are several blocks of such summaries.
17! When you run the program, the last line in the output file should
18! read: 'No errors occurred.' Please contact one of the authors if
19! any errors do occur.
20
21! Reference:
22! ALGORITHM 648, COLLECTED ALGORITHMS FROM ACM.
23! TRANSACTIONS ON MATHEMATICAL SOFTWARE,
24! VOL. 13, NO. 1, P. 28
25
26! Caution:
27! The original test routines are altered in this version.
28! Some subroutine arguments, COMMON variables, and local
29! variables have been moved to the following global header.
30! You should obtain the original test suite from netlib
31! if you plan to use the routines for other purposes.
32
33    MODULE STIFFSET
34
35! Note:
36! In this version ID and IID are defined in the main
37! program demostiff (not in the Toronto routines).
38! The other parameters are obtained by calling IVALU
39! with a modified argument list.
40
41      IMPLICIT NONE
42      INTEGER IWT, N, ID, IID, NFCN, NJAC, NLUD
43      DOUBLE PRECISION W, DY
44      DIMENSION W(20), DY(400)
45
46    CONTAINS
47
48      SUBROUTINE DERIVS(NEQ,T,Y,YDOT)
49!       Define the system derivatives.
50        IMPLICIT NONE
51        INTEGER NEQ
52        DOUBLE PRECISION T, Y, YDOT
53        DIMENSION Y(NEQ), YDOT(NEQ)
54        CALL FCN(T,Y,YDOT)
55        RETURN
56      END SUBROUTINE DERIVS
57
58      SUBROUTINE JACD(NEQ,T,Y,ML,MU,PD,NROWPD)
59!       Load the Jacobian as a dense matrix.
60        IMPLICIT NONE
61        INTEGER NEQ, ML, MU, NROWPD, I, J
62        DOUBLE PRECISION T, Y, PD
63        DIMENSION Y(NEQ), PD(NROWPD,NEQ)
64        CALL PDERV(T,Y)
65        DO J = 1, NEQ
66          DO I = 1, NEQ
67            PD(I,J) = DY(I+(J-1)*NROWPD)
68          END DO
69        END DO
70        RETURN
71      END SUBROUTINE JACD
72
73      SUBROUTINE JACB(NEQ,T,Y,ML,MU,PD,NROWPD)
74!       Load the Jacobian as a banded matrix.
75        IMPLICIT NONE
76        INTEGER NEQ, ML, MU, NROWPD, I, J, K
77        DOUBLE PRECISION T, Y, PD
78        DIMENSION Y(NEQ), PD(NROWPD,NEQ)
79        CALL PDERV(T,Y)
80        DO J = 1, NEQ
81           DO I = 1, NEQ
82              K = I - J + MU + 1
83              PD(K,J) = DY(I+(J-1)*NEQ)
84           END DO
85        END DO
86      END SUBROUTINE JACB
87
88      SUBROUTINE JACS(NEQ,T,Y,IA,JA,NZ,P)
89!       Load the Jacobian as a sparse matrix.
90        IMPLICIT NONE
91        INTEGER NEQ, IA, JA, NZ, ML, MU, COL, ROW, I, NROWPD
92!       INTEGER NZSAVE
93        DOUBLE PRECISION T, Y, P
94        DIMENSION Y(*), IA(*), JA(*), P(*)
95        IF (NZ<=0) THEN
96          NZ = NEQ*NEQ
97          RETURN
98        END IF
99        ML = NEQ - 1
100        MU = NEQ - 1
101        NROWPD = NEQ
102        CALL JACD(NEQ,T,Y,ML,MU,P,NROWPD)
103        IA(1) = 1
104        DO I = 1, NEQ
105          IA(I+1) = IA(I) + NEQ
106        END DO
107        I = 0
108        DO COL = 1, NEQ
109          I = NEQ*(COL-1)
110          DO ROW = 1, NEQ
111            I = I + 1
112            JA(I) = ROW
113          END DO
114        END DO
115        RETURN
116      END SUBROUTINE JACS
117
118      SUBROUTINE IVALU(XSTART,XEND,HBEGIN,HMAX,Y)
119
120!      ROUTINE TO PROVIDE THE INITIAL VALUES REQUIRED TO SPECIFY
121!      THE MATHEMATICAL PROBLEM AS WELL AS VARIOUS PROBLEM
122!      PARAMETERS REQUIRED BY THE TESTING PACKAGE. THE APPROPRIATE
123!      SCALING VECTOR IS ALSO INITIALISED IN CASE THIS OPTION IS
124!      SELECTED.
125
126!      PARAMETERS (OUTPUT)
127!         N      - DIMENSION OF THE PROBLEM
128!         XSTART - INITIAL VALUE OF THE INDEPENDENT VARIABLE
129!         XEND   - FINAL VALUE OF THE INDEPENDENT VARIABLE
130!         HBEGIN - APPROPRIATE STARTING STEPSIZE
131!         Y      - VECTOR OF INITIAL CONDITIONS FOR THE DEPENDENT
132!                  VARIABLES
133!         W      - VECTOR OF WEIGHTS USED TO SCALE THE PROBLEM IF
134!                  THIS OPTION IS SELECTED.
135
136!      PARAMETER  (INPUT)
137!         IWT    - FLAG TO INDICATE IF SCALED OPTION IS SELECTED
138!         ID     - FLAG IDENTIFYING WHICH EQUATION IS BEING SOLVED
139
140        IMPLICIT NONE
141
142!     .. Scalar Arguments ..
143        DOUBLE PRECISION HBEGIN, HMAX, XEND, XSTART
144!     .. Array Arguments ..
145        DOUBLE PRECISION Y(20)
146!     .. Local Scalars ..
147        DOUBLE PRECISION XS
148        INTEGER I, IOUT, ITMP
149!     .. Data statements ..
150        DATA XS/0.D0/
151
152!     .. Executable Statements ..
153
154        XSTART = XS
155!     GOTO (40, 80, 120, 160, 20, 20, 20, 20, 20, 20, 200, 220, 220,    &
156!     220, 220, 20, 20, 20, 20, 20, 360, 400, 400, 400, 400, 20, 20, 20,&
157!     20, 20, 540, 580, 600, 640, 660, 680, 20, 20, 20, 20, 700, 740,   &
158!     760, 780, 800, 20, 20, 20, 20, 20, 840, 860, 880, 900, 920) ID
159        GOTO (20,30,40,50,10,10,10,10,10,10,60,70,70,70,70,10,10,10,10,10, &
160          130,140,140,140,140,10,10,10,10,10,200,210,220,230,240,250,10,10,10, &
161          10,260,270,280,290,300,10,10,10,10,10,310,320,330,340,350) ID
16210      IOUT = 6
163        WRITE (IOUT,FMT=90000) ID
164        STOP
165
166!     PROBLEM CLASS A - LINEAR WITH REAL EIGENVALUES
167
16820      CONTINUE
169!     PROBLEM A1
170        N = 4
171        W(1) = 0.100D+01
172        W(2) = 0.100D+01
173        W(3) = 0.100D+01
174        W(4) = 0.100D+01
175        XEND = 20.D0
176        HBEGIN = 1.0D-2
177        HMAX = 20.D0
178        DO I = 1, N
179          Y(I) = 1.0D0
180        END DO
181        GOTO 360
182
18330      CONTINUE
184!     PROBLEM A2
185        N = 9
186        W(1) = 0.100D+00
187        W(2) = 0.200D+00
188        W(3) = 0.300D+00
189        W(4) = 0.400D+00
190        W(5) = 0.500D+00
191        W(6) = 0.600D+00
192        W(7) = 0.700D+00
193        W(8) = 0.800D+00
194        W(9) = 0.900D+00
195        XEND = 120.D0
196        HBEGIN = 5.D-4
197        HMAX = 120.D0
198        DO I = 1, N
199          Y(I) = 0.D0
200        END DO
201        GOTO 360
202
20340      CONTINUE
204!     PROBLEM A3
205        N = 4
206        W(1) = 0.100D+01
207        W(2) = 0.100D+01
208        W(3) = 0.782D+01
209        W(4) = 0.100D+01
210        HBEGIN = 1.D-5
211        XEND = 20.D0
212        HMAX = 20.D0
213        DO I = 1, N
214          Y(I) = 1.D0
215        END DO
216        GOTO 360
217
21850      CONTINUE
219!     PROBLEM A4
220        N = 10
221        W(1) = 0.100D+01
222        W(2) = 0.100D+01
223        W(3) = 0.100D+01
224        W(4) = 0.100D+01
225        W(5) = 0.100D+01
226        W(6) = 0.100D+01
227        W(7) = 0.100D+01
228        W(8) = 0.100D+01
229        W(9) = 0.100D+01
230        W(10) = 0.100D+01
231        XEND = 1.D0
232        HBEGIN = 1.D-5
233        HMAX = 1.D0
234        DO I = 1, N
235          Y(I) = 1.D0
236        END DO
237        GOTO 360
238
239!     PROBLEM CLASS B - LINEAR WITH NON-REAL EIGENVALUES
240
24160      CONTINUE
242!     PROBLEM B1
243        N = 4
244        W(1) = 0.100D+01
245        W(2) = 0.859D+01
246        W(3) = 0.100D+01
247        W(4) = 0.322D+02
248        XEND = 20.D0
249        HBEGIN = 7.D-3
250        HMAX = 20.D0
251        Y(1) = 1.D0
252        Y(2) = 0.D0
253        Y(3) = 1.D0
254        Y(4) = 0.D0
255        GOTO 360
256
25770      CONTINUE
258!     PROBLEM B2, B3, B4, B5
259        N = 6
260        ITMP = IID - 1
261        GOTO (80,90,100,110) ITMP
26280      CONTINUE
263        W(1) = 0.100D+01
264        W(2) = 0.100D+01
265        W(3) = 0.100D+01
266        W(4) = 0.100D+01
267        W(5) = 0.100D+01
268        W(6) = 0.100D+01
269        GOTO 120
27090      CONTINUE
271        W(1) = 0.100D+01
272        W(2) = 0.100D+01
273        W(3) = 0.100D+01
274        W(4) = 0.100D+01
275        W(5) = 0.100D+01
276        W(6) = 0.100D+01
277        GOTO 120
278100     CONTINUE
279        W(1) = 0.112D+01
280        W(2) = 0.100D+01
281        W(3) = 0.100D+01
282        W(4) = 0.100D+01
283        W(5) = 0.100D+01
284        W(6) = 0.100D+01
285        GOTO 120
286110     CONTINUE
287        W(1) = 0.131D+01
288        W(2) = 0.112D+01
289        W(3) = 0.100D+01
290        W(4) = 0.100D+01
291        W(5) = 0.100D+01
292        W(6) = 0.100D+01
293120     CONTINUE
294        XEND = 20.D0
295        HBEGIN = 1.D-2
296        HMAX = 20.D0
297        DO I = 1, N
298          Y(I) = 1.D0
299        END DO
300        GOTO 360
301
302!     PROBLEM CLASS C - NON-LINEAR COUPLING FROM
303!                       STEADY STATE TO TRANSIENT
304
305130     CONTINUE
306!     PROBLEM C1
307        N = 4
308        W(1) = 0.102D+01
309        W(2) = 0.103D+01
310        W(3) = 0.100D+01
311        W(4) = 0.100D+01
312        XEND = 20.D0
313        HBEGIN = 1.D-2
314        HMAX = 20.D0
315        DO I = 1, N
316          Y(I) = 1.D0
317        END DO
318        GOTO 360
319
320140     CONTINUE
321!     PROBLEM C2, C3, C4, C5
322        N = 4
323        ITMP = IID - 1
324        GOTO (150,160,170,180) ITMP
325150     CONTINUE
326        W(1) = 0.200D+01
327        W(2) = 0.100D+01
328        W(3) = 0.100D+01
329        W(4) = 0.100D+01
330        GOTO 190
331160     CONTINUE
332        W(1) = 0.200D+01
333        W(2) = 0.100D+01
334        W(3) = 0.100D+01
335        W(4) = 0.100D+01
336        GOTO 190
337170     CONTINUE
338        W(1) = 0.200D+01
339        W(2) = 0.400D+01
340        W(3) = 0.200D+02
341        W(4) = 0.420D+03
342        GOTO 190
343180     CONTINUE
344        W(1) = 0.200D+01
345        W(2) = 0.800D+01
346        W(3) = 0.136D+03
347        W(4) = 0.371D+05
348190     CONTINUE
349        XEND = 20.D0
350        HBEGIN = 1.D-2
351        HMAX = 20.D0
352        DO I = 1, N
353          Y(I) = 1.D0
354        END DO
355        GOTO 360
356
357!     PROBLEM CLASS D - NON-LINEAR WITH REAL EIGENVALUES
358
359200     CONTINUE
360!     PROBLEM D1
361        N = 3
362        W(1) = 0.223D+02
363        W(2) = 0.271D+02
364        W(3) = 0.400D+03
365        XEND = 400.D0
366        HBEGIN = 1.7D-2
367        HMAX = 400.D0
368        DO I = 1, N
369          Y(I) = 0.D0
370        END DO
371        GOTO 360
372
373210     CONTINUE
374!     PROBLEM D2
375        N = 3
376        W(1) = 0.100D+01
377        W(2) = 0.365D+00
378        W(3) = 0.285D+02
379        XEND = 40.D0
380        HBEGIN = 1.D-5
381        HMAX = 40.D0
382        Y(1) = 1.D0
383        Y(2) = 0.D0
384        Y(3) = 0.D0
385        GOTO 360
386
387220     CONTINUE
388!     PROBLEM D3
389        N = 4
390        W(1) = 0.100D+01
391        W(2) = 0.100D+01
392        W(3) = 0.360D+00
393        W(4) = 0.485D+00
394        XEND = 20.D0
395        HBEGIN = 2.5D-5
396        HMAX = 20.D0
397        DO I = 1, 2
398          Y(I) = 1.D0
399          Y(I+2) = 0.D0
400        END DO
401        GOTO 360
402
403230     CONTINUE
404!     PROBLEM D4
405        N = 3
406        W(1) = 0.100D+01
407        W(2) = 0.142D+01
408        W(3) = 0.371D-05
409        XEND = 50.D0
410        HBEGIN = 2.9D-4
411        HMAX = 50.D0
412        Y(1) = 1.D0
413        Y(2) = 1.D0
414        Y(3) = 0.D0
415        GOTO 360
416
417240     CONTINUE
418!     PROBLEM D5
419        N = 2
420        W(1) = 0.992D+00
421        W(2) = 0.984D+00
422        XEND = 1.D2
423        HBEGIN = 1.D-4
424        HMAX = 1.D2
425        Y(1) = 0.D0
426        Y(2) = 0.D0
427        GOTO 360
428
429250     CONTINUE
430!     PROBLEM D6
431        N = 3
432        W(1) = 0.100D+01
433        W(2) = 0.148D+00
434        W(3) = 0.577D-07
435        XEND = 1.D0
436        HBEGIN = 3.3D-8
437        HMAX = 1.D0
438        Y(1) = 1.D0
439        Y(2) = 0.D0
440        Y(3) = 0.D0
441        GOTO 360
442
443!     PROBLEM CLASS E - NON-LINEAR WITH NON-REAL EIGENVALUES
444
445260     CONTINUE
446!     PROBLEM E1
447        N = 4
448        W(1) = 0.100D-07
449        W(2) = 0.223D-06
450        W(3) = 0.132D-04
451        W(4) = 0.171D-02
452        XEND = 1.D0
453        HBEGIN = 6.8D-3
454        HMAX = 1.D0
455        DO I = 1, N
456          Y(I) = 0.D0
457        END DO
458        GOTO 360
459
460270     CONTINUE
461!     PROBLEM E2
462        N = 2
463        W(1) = 0.202D+01
464        W(2) = 0.764D+01
465        XEND = 1.D1
466        HBEGIN = 1.D-3
467        HMAX = 1.D1
468        Y(1) = 2.D0
469        Y(2) = 0.D0
470        GOTO 360
471
472280     CONTINUE
473!     PROBLEM E3
474        N = 3
475        W(1) = 0.163D+01
476        W(2) = 0.160D+01
477        W(3) = 0.263D+02
478        XEND = 5.D2
479        HBEGIN = .2D-1
480        HMAX = 5.D2
481        Y(1) = 1.D0
482        Y(2) = 1.D0
483        Y(3) = 0.D0
484        GOTO 360
485
486290     CONTINUE
487!     PROBLEM E4
488        N = 4
489        W(1) = 0.288D+02
490        W(2) = 0.295D+02
491        W(3) = 0.155D+02
492        W(4) = 0.163D+02
493        XEND = 1.D3
494        HBEGIN = 1.D-3
495        HMAX = 1.D3
496        Y(1) = 0.D0
497        Y(2) = -2.D0
498        Y(3) = -1.D0
499        Y(4) = -1.D0
500        GOTO 360
501
502300     CONTINUE
503!     PROBLEM E5
504        N = 4
505        W(1) = 0.176D-02
506        W(2) = 0.146D-09
507        W(3) = 0.827D-11
508        W(4) = 0.138D-09
509        XEND = 1.D3
510        HBEGIN = 5.D-5
511        HMAX = 1.D3
512        Y(1) = 1.76D-3
513        DO I = 2, N
514          Y(I) = 0.D0
515        END DO
516        GOTO 360
517
518!     PROBLEM CLASS F - CHEMICAL KINETICS EQUATIONS
519
520310     CONTINUE
521!     PROBLEM F1
522        N = 4
523        W(1) = 0.121D+04
524        W(2) = 0.835D-01
525        W(3) = 0.121D+04
526        W(4) = 0.100D+00
527        HMAX = 1.D3
528        HBEGIN = 1.D-4
529        XEND = 1.D3
530        Y(1) = 761.D0
531        Y(2) = 0.D0
532        Y(3) = 600.D0
533        Y(4) = .1D0
534        GOTO 360
535
536320     CONTINUE
537!     PROBLEM F2
538        N = 2
539        W(1) = 0.100D+01
540        W(2) = 0.253D-02
541        HMAX = 240.D0
542        HBEGIN = 1.D-2
543        XEND = 240.D0
544        Y(1) = 1.0D0
545        Y(2) = 0.D0
546        GOTO 360
547
548330     CONTINUE
549!     PROBLEM F3
550        N = 5
551        W(1) = 0.400D-05
552        W(2) = 0.100D-05
553        W(3) = 0.374D-08
554        W(4) = 0.765D-06
555        W(5) = 0.324D-05
556        HBEGIN = 1.D-6
557        HBEGIN = 1.D-10
558        HMAX = 100.D0
559        XEND = 100.D0
560        Y(1) = 4.D-6
561        Y(2) = 1.D-6
562        Y(3) = 0.0D0
563        Y(4) = 0.0D0
564        Y(5) = 0.0D0
565        GOTO 360
566
567340     CONTINUE
568!     PROBLEM F4
569        N = 3
570        W(1) = 0.118D+06
571        W(2) = 0.177D+04
572        W(3) = 0.313D+05
573        HBEGIN = 1.D-3
574!       HMAX = 50.D0
575        HMAX = 1.D0
576        XEND = 300.D0
577        Y(1) = 4.D0
578        Y(2) = 1.1D0
579        Y(3) = 4.D0
580        GOTO 360
581
582350     CONTINUE
583!     PROBLEM F5
584        N = 4
585        W(1) = 0.336D-06
586        W(2) = 0.826D-02
587        W(3) = 0.619D-02
588        W(4) = 0.955D-05
589        HBEGIN = 1.D-7
590        HMAX = 100.D0
591        XEND = 100.D0
592        Y(1) = 3.365D-7
593        Y(2) = 8.261D-3
594        Y(3) = 1.642D-3
595        Y(4) = 9.380D-6
596360     CONTINUE
597        IF (IWT<0) GOTO 370
598        DO I = 1, N
599          Y(I) = Y(I)/W(I)
600        END DO
601370     CONTINUE
602        RETURN
603
60490000   FORMAT (' AN INVALID INTERNAL PROBLEM ID OF ',I4, &
605          ' WAS FOUND BY THE IVALU ROUTINE',' RUN TERMINATED. CHECK THE DATA.' &
606          )
607      END SUBROUTINE IVALU
608
609      SUBROUTINE EVALU(Y)
610
611!     ROUTINE TO PROVIDE THE 'TRUE' SOLUTION OF THE DIFFERENTIAL
612!     EQUATION EVALUATED AT THE ENDPOINT OF THE INTEGRATION.
613
614!     1986 REVISION:  SOME VERY SMALL CONSTANTS HAVE BEEN RECAST IN THE
615!     (NOT SO SMALL CONST)/(1.E38) TO AVOID COMPILE-TIME UNDERFLOW ERROR
616!     IT IS ASSUMED 1E+38 WON'T OVERFLOW.
617!     PARAMETER  (OUTPUT)
618!        Y      - THE TRUE SOLUTION VECTOR EVALUATED AT THE ENDPOINT
619
620!     PARAMETERS (INPUT)
621!        N      - DIMENSION OF THE PROBLEM
622!        W      - VECTOR OF WEIGHTS USED TO SCALE THE PROBLEM
623!                 IF THIS OPTION IS SELECTED
624!        IWT    - FLAG USED TO SIGNAL WHEN THE SCALED PROBLEM IS
625!                 BEING SOLVED
626!        ID     - FLAG USED TO INDICATE WHICH EQUATION IS BEING
627!                 SOLVED
628
629        IMPLICIT NONE
630
631!     .. Parameters ..
632        DOUBLE PRECISION TENE38
633        PARAMETER (TENE38=1.D38)
634!     .. Array Arguments ..
635        DOUBLE PRECISION Y(20)
636!     .. Local Scalars ..
637        INTEGER I
638
639!     .. Executable Statements ..
640
641!     GOTO (20, 40, 60, 80, 620, 620, 620, 620, 620, 620, 100, 120, 140,&
642!     160, 180, 620, 620, 620, 620, 620, 200, 220, 240, 260, 280, 620,  &
643!     620, 620, 620, 620, 300, 320, 340, 360, 380, 400, 620, 620, 620,  &
644!     620, 420, 440, 460, 480, 500, 620, 620, 620, 620, 620, 520, 540,  &
645!     560, 580, 600, 620, 620, 620, 620, 620) ID
646        GOTO (10,20,30,40,310,310,310,310,310,310,50,60,70,80,90,310,310,310, &
647          310,310,100,110,120,130,140,310,310,310,310,310,150,160,170,180,190, &
648          200,310,310,310,310,210,220,230,240,250,310,310,310,310,310,260,270, &
649          280,290,300,310,310,310,310,310) ID
650        GOTO 310
651
652!     PROBLEM CLASS A
653
654!     PROBLEM A1
65510      Y(1) = 4.539992969929191D-05
656        Y(2) = 2.061153036149920D-09
657        Y(3) = 2.823006338263857D-18/TENE38
658        Y(4) = 5.235792540515384D-14/TENE38
659        GOTO 310
660
661!     PROBLEM A2
66220      Y(1) = 9.999912552999704D-02
663        Y(2) = 1.999982511586291D-01
664        Y(3) = 2.999975543202422D-01
665        Y(4) = 3.999971057541257D-01
666        Y(5) = 4.999969509963023D-01
667        Y(6) = 5.999971057569546D-01
668        Y(7) = 6.999975543256127D-01
669        Y(8) = 7.999982511659962D-01
670        Y(9) = 8.999991255386128D-01
671        GOTO 310
672
673!     PROBLEM A3
67430      Y(1) = -1.353352661867235D-03
675        Y(2) = 1.368526917891521D-02
676        Y(3) = 1.503725348455117D+00
677        Y(4) = 1.353352832366099D-01
678        GOTO 310
679
680!     PROBLEM A4
68140      Y(1) = 3.678794411714325D-01
682        Y(2) = 1.265870722340194D-14
683        Y(3) = 1.911533219339204D-04/TENE38
684        Y(4) = 2.277441666729596D-17/TENE38
685        Y(5) = 0.0D0
686        Y(6) = 0.0D0
687        Y(7) = 0.0D0
688        Y(8) = 0.0D0
689        Y(9) = 0.0D0
690        Y(10) = 0.0D0
691        GOTO 310
692
693!     PROBLEM CLASS B
694
695!     PROBLEM B1
69650      Y(1) = 1.004166730990124D-09
697        Y(2) = 1.800023280346500D-08
698        Y(3) = 0.0D0
699        Y(4) = -6.042962877027475D-03/TENE38/TENE38
700        GOTO 310
701
702!     PROBLEM B2
70360      Y(1) = 6.181330838820067D-31
704        Y(2) = 8.963657877626303D-31
705        Y(3) = 2.738406773453261D-27
706        Y(4) = 2.061153063164016D-09
707        Y(5) = 4.539992973654118D-05
708        Y(6) = 1.353352832365270D-01
709        GOTO 310
710
711!     PROBLEM B3
71270      Y(1) = -1.076790816984970D-28
713        Y(2) = 5.455007683862160D-28
714        Y(3) = 2.738539964946867D-27
715        Y(4) = 2.061153071123456D-09
716        Y(5) = 4.539992974611305D-05
717        Y(6) = 1.353352832365675D-01
718        GOTO 310
719
720!     PROBLEM B4
72180      Y(1) = 1.331242472678293D-22
722        Y(2) = -2.325916064237926D-22
723        Y(3) = 1.517853928534857D-35
724        Y(4) = 2.061152428936651D-09
725        Y(5) = 4.539992963392291D-05
726        Y(6) = 1.353352832363442D-01
727        GOTO 310
728
729!     PROBLEM B5
73090      Y(1) = -3.100634584292190D-14
731        Y(2) = 3.862788998076547D-14
732        Y(3) = 1.804851385304217D-35
733        Y(4) = 2.061153622425655D-09
734        Y(5) = 4.539992976246673D-05
735        Y(6) = 1.353352832366126D-01
736        GOTO 310
737
738!     PROBLEM CLASS C
739
740!     PROBLEM C1
741100     Y(1) = 4.003223925456179D-04
742        Y(2) = 4.001600000000000D-04
743        Y(3) = 4.000000000000000D-04
744        Y(4) = 2.000000000000000D-02
745        GOTO 310
746
747!     PROBLEM C2
748110     Y(1) = 1.999999997938994D+00
749        Y(2) = 3.999999990839974D-02
750        Y(3) = 4.001599991537078D-02
751        Y(4) = 4.003201271914461D-02
752        GOTO 310
753
754!     PROBLEM C3
755120     Y(1) = 1.999999997939167D+00
756        Y(2) = 3.999999990840744D-01
757        Y(3) = 4.159999990793773D-01
758        Y(4) = 4.333055990159567D-01
759        GOTO 310
760
761!     PROBLEM C4
762130     Y(1) = 1.999999997938846D+00
763        Y(2) = 3.999999990839318D+00
764        Y(3) = 1.999999991637941D+01
765        Y(4) = 4.199999965390368D+02
766        GOTO 310
767
768!     PROBLEM C5
769140     Y(1) = 1.999999997938846D+00
770        Y(2) = 7.999999981678634D+00
771        Y(3) = 1.359999993817714D+02
772        Y(4) = 3.712799965967762D+04
773        GOTO 310
774
775!     PROBLEM CLASS D
776
777!     PROBLEM D1
778150     Y(1) = 2.224222010616901D+01
779        Y(2) = 2.711071334484136D+01
780        Y(3) = 3.999999999999999D+02
781        GOTO 310
782
783!     PROBLEM D2
784160     Y(1) = 7.158270687193941D-01
785        Y(2) = 9.185534764557338D-02
786        Y(3) = 2.841637457458413D+01
787        GOTO 310
788
789!     PROBLEM D3
790170     Y(1) = 6.397604446889910D-01
791        Y(2) = 5.630850708287990D-03
792        Y(3) = 3.602395553110090D-01
793        Y(4) = 3.170647969903515D-01
794        GOTO 310
795
796!     PROBLEM D4
797180     Y(1) = 5.976546980673215D-01
798        Y(2) = 1.402343408546138D+00
799        Y(3) = -1.893386540441913D-06
800        GOTO 310
801
802!     PROBLEM D5
803190     Y(1) = -9.916420698713913D-01
804        Y(2) = 9.833363588544478D-01
805        GOTO 310
806
807!     PROBLEM D6
808200     Y(1) = 8.523995440749948D-01
809        Y(2) = 1.476003981941319D-01
810        Y(3) = 5.773087333950041D-08
811        GOTO 310
812
813!     PROBLEM CLASS E
814
815!     PROBLEM E1
816210     Y(1) = 1.000000000000012D-08
817        Y(2) = -1.625323873316817D-19
818        Y(3) = 2.025953375595861D-17
819        Y(4) = -1.853149807630002D-15
820        GOTO 310
821
822!     PROBLEM E2
823220     Y(1) = -1.158701266031984D+00
824        Y(2) = 4.304698089780476D-01
825        GOTO 310
826
827!     PROBLEM E3
828230     Y(1) = 4.253052197643089D-03
829        Y(2) = 5.317019548450387D-03
830        Y(3) = 2.627647748753926D+01
831        GOTO 310
832
833!     PROBLEM E4
834240     Y(1) = 1.999999977523654D+01
835        Y(2) = -2.000000022476345D+01
836        Y(3) = -2.247634567084293D-07
837        Y(4) = 2.247634567084293D-07
838        GOTO 310
839
840!     PROBLEM E5
841250     Y(1) = 1.618076919919600D-03
842        Y(2) = 1.382236955418478D-10
843        Y(3) = 8.251573436034144D-12
844        Y(4) = 1.299721221058136D-10
845        GOTO 310
846
847!     PROBLEM CLASS F
848
849!     PROBLEM F1
850260     Y(1) = 1.211129474696585D+03
851        Y(2) = 1.271123619113051D-05
852        Y(3) = 1.208637804660361D+03
853        Y(4) = 3.241981171933418D-04
854        GOTO 310
855
856!     PROBLEM F2
857270     Y(1) = 3.912699122292088D-01
858        Y(2) = 1.329964166084866D-03
859        GOTO 310
860
861!     PROBLEM F3
862280     Y(1) = 3.235910070806680D-13
863        Y(2) = 2.360679774997897D-07
864        Y(3) = 7.639319089351045D-14
865        Y(4) = 7.639319461070194D-07
866        Y(5) = 3.236067653908783D-06
867        GOTO 310
868
869!     PROBLEM F4
870290     Y(1) = 4.418303324022590D+00
871        Y(2) = 1.290244712916425D+00
872        Y(3) = 3.019282584050490D+00
873        GOTO 310
874
875!     PROBLEM F5
876300     Y(1) = 1.713564284690712D-07
877        Y(2) = 3.713563071160676D-03
878        Y(3) = 6.189271785267793D-03
879        Y(4) = 9.545143571530929D-06
880310     CONTINUE
881        IF (IWT<0) GOTO 320
882        DO I = 1, N
883          Y(I) = Y(I)/W(I)
884        END DO
885320     CONTINUE
886        RETURN
887      END SUBROUTINE EVALU
888
889      SUBROUTINE FCN(X,Y,YP)
890
891!     ROUTINE TO EVALUATE THE DERIVATIVE F(X,Y) CORRESPONDING TO
892!     THE DIFFERENTIAL EQUATION:
893!                    DY/DX = F(X,Y) .
894!     THE ROUTINE STORES THE VECTOR OF DERIVATIVES IN YP(*). THE
895!     PARTICULAR EQUATION BEING INTEGRATED IS INDICATED BY THE
896!     VALUE OF THE FLAG ID WHICH IS PASSED THROUGH COMMON. THE
897!     DIFFERENTIAL EQUATION IS SCALED BY THE WEIGHT VECTOR W(*)
898!     IF THIS OPTION HAS BEEN SELECTED (IF SO IT IS SIGNALLED
899!     BY THE FLAG IWT).
900
901        IMPLICIT NONE
902
903!     .. Scalar Arguments ..
904        DOUBLE PRECISION X
905!     .. Array Arguments ..
906        DOUBLE PRECISION Y(20), YP(20)
907!     .. Local Scalars ..
908        DOUBLE PRECISION F, Q, S, SUM, T, TEMP, XTEMP
909        INTEGER I
910!     .. Local Arrays ..
911        DOUBLE PRECISION BPARM(4), CPARM(4), VECT1(4), VECT2(4), YTEMP(20)
912!     .. Data statements ..
913        DATA BPARM/3.D0, 8.D0, 25.D0, 1.D2/
914        DATA CPARM/1.D-1, 1.D0, 1.D1, 2.D1/
915
916!     .. Executable Statements ..
917
918        NFCN = NFCN + 1
919        IF (IWT<0) GOTO 10
920        DO I = 1, N
921          YTEMP(I) = Y(I)
922          Y(I) = Y(I)*W(I)
923        END DO
92410      CONTINUE
925        GOTO (20,30,40,50,260,260,260,260,260,260,60,70,70,70,70,260,260,260, &
926          260,260,80,90,90,90,90,260,260,260,260,260,100,110,120,130,140,150, &
927          260,260,260,260,160,170,180,190,200,260,260,260,260,260,210,220,230, &
928          240,250) ID
929        GOTO 260
930
931!     PROBLEM CLASS A - LINEAR WITH REAL EIGENVALUES
932
933!     PROBLEM A1
93420      YP(1) = -.5D0*Y(1)
935        YP(2) = -1.D0*Y(2)
936        YP(3) = -1.D2*Y(3)
937        YP(4) = -9.D1*Y(4)
938        GOTO 260
939
940!     PROBLEM A2
94130      YP(1) = -1.8D3*Y(1) + 9.D2*Y(2)
942        DO I = 2, 8
943          YP(I) = Y(I-1) - 2.D0*Y(I) + Y(I+1)
944        END DO
945        YP(9) = 1.D3*Y(8) - 2.D3*Y(9) + 1.D3
946        GOTO 260
947
948!     PROBLEM A3
94940      YP(1) = -1.D4*Y(1) + 1.D2*Y(2) - 1.D1*Y(3) + 1.D0*Y(4)
950        YP(2) = -1.D3*Y(2) + 1.D1*Y(3) - 1.D1*Y(4)
951        YP(3) = -1.D0*Y(3) + 1.D1*Y(4)
952        YP(4) = -1.D-1*Y(4)
953        GOTO 260
954
955!     PROBLEM A4
95650      DO I = 1, 10
957          YP(I) = -REAL(I)**5*Y(I)
958        END DO
959        GOTO 260
960
961!     PROBLEM CLASS B - LINEAR WITH NON-REAL EIGENVALUES
962
963!     PROBLEM B1
96460      YP(1) = -Y(1) + Y(2)
965        YP(2) = -1.D2*Y(1) - Y(2)
966        YP(3) = -1.D2*Y(3) + Y(4)
967        YP(4) = -1.D4*Y(3) - 1.D2*Y(4)
968        GOTO 260
969
970!     PROBLEMS B2, B3, B4, B5
97170      YP(1) = -1.D1*Y(1) + BPARM(IID-1)*Y(2)
972        YP(2) = -BPARM(IID-1)*Y(1) - 1.D1*Y(2)
973        YP(3) = -4.D0*Y(3)
974        YP(4) = -1.D0*Y(4)
975        YP(5) = -.5D0*Y(5)
976        YP(6) = -.1D0*Y(6)
977        GOTO 260
978
979!     PROBLEM CLASS C - NON-LINEAR COUPLING FROM
980!                       STEADY STATE TO TRANSIENT
981
982!     PROBLEM C1
98380      YP(1) = -Y(1) + (Y(2)*Y(2)+Y(3)*Y(3)+Y(4)*Y(4))
984        YP(2) = -1.D1*Y(2) + 1.D1*(Y(3)*Y(3)+Y(4)*Y(4))
985        YP(3) = -4.D1*Y(3) + 4.D1*Y(4)*Y(4)
986        YP(4) = -1.D2*Y(4) + 2.D0
987        GOTO 260
988
989!     PROBLEMS C2, C3, C4, C5
99090      YP(1) = -Y(1) + 2.D0
991        YP(2) = -1.D1*Y(2) + CPARM(IID-1)*Y(1)*Y(1)
992        YP(3) = -4.D1*Y(3) + (Y(1)*Y(1)+Y(2)*Y(2))*CPARM(IID-1)*4.D0
993        YP(4) = (Y(1)*Y(1)+Y(2)*Y(2)+Y(3)*Y(3))*CPARM(IID-1)*1.D1 - 1.D2*Y(4)
994        GOTO 260
995
996!     PROBLEM CLASS D - NON-LINEAR WITH REAL EIGENVALUES
997
998!     PROBLEM D1
999100     YP(1) = .2D0*Y(2) - .2D0*Y(1)
1000        YP(2) = 1.D1*Y(1) - (6.D1-.125D0*Y(3))*Y(2) + .125D0*Y(3)
1001        YP(3) = 1.D0
1002        GOTO 260
1003
1004!     PROBLEM D2
1005110     YP(1) = -.04D0*Y(1) + .01D0*Y(2)*Y(3)
1006        YP(2) = 4.D2*Y(1) - 1.D2*Y(2)*Y(3) - 3.D3*Y(2)**2
1007        YP(3) = 3.D1*Y(2)**2
1008        GOTO 260
1009
1010!     PROBLEM D3
1011120     YP(1) = Y(3) - 1.D2*Y(1)*Y(2)
1012        YP(3) = -YP(1)
1013        YP(4) = -Y(4) + 1.D4*Y(2)**2
1014        YP(2) = YP(1) - YP(4) + Y(4) - 1.D4*Y(2)**2
1015        GOTO 260
1016
1017!     PROBLEM D4
1018130     YP(1) = -.013D0*Y(1) - 1.D3*Y(1)*Y(3)
1019        YP(2) = -2.5D3*Y(2)*Y(3)
1020        YP(3) = YP(1) + YP(2)
1021        GOTO 260
1022
1023!     PROBLEM D5
1024140     XTEMP = .01D0 + Y(1) + Y(2)
1025        YP(1) = .01D0 - XTEMP*(1.D0+(Y(1)+1.D3)*(Y(1)+1.D0))
1026        YP(2) = .01D0 - XTEMP*(1.D0+Y(2)**2)
1027        GOTO 260
1028
1029!     PROBLEM D6
1030150     YP(1) = -Y(1) + 1.D8*Y(3)*(1.D0-Y(1))
1031        YP(2) = -1.D1*Y(2) + 3.D7*Y(3)*(1.D0-Y(2))
1032        YP(3) = -YP(1) - YP(2)
1033        GOTO 260
1034
1035!     PROBLEM CLASS E - NON-LINEAR WITH NON-REAL EIGENVALUES
1036
1037!     PROBLEM E1
1038160     YP(1) = Y(2)
1039        YP(2) = Y(3)
1040        YP(3) = Y(4)
1041        YP(4) = (Y(1)**2-SIN(Y(1))-1.D8)*Y(1) + (Y(2)*Y(3)/(Y(1)**2+1.D0)-4.D6 &
1042          )*Y(2) + (1.D0-6.D4)*Y(3) + (1.D1*EXP(-Y(4)**2)-4.D2)*Y(4) + 1.D0
1043        GOTO 260
1044
1045!     PROBLEM E2
1046170     YP(1) = Y(2)
1047        YP(2) = 5.D0*Y(2) - 5.D0*Y(1)*Y(1)*Y(2) - Y(1)
1048        GOTO 260
1049
1050!     PROBLEM E3
1051180     YP(1) = -55.D0*Y(1) - Y(3)*Y(1) + 65.D0*Y(2)
1052        YP(2) = .785D-1*Y(1) - .785D-1*Y(2)
1053        YP(3) = .1D0*Y(1)
1054        GOTO 260
1055
1056!     PROBLEM E4
1057190     SUM = Y(1) + Y(2) + Y(3) + Y(4)
1058        DO I = 1, 4
1059          VECT2(I) = -Y(I) + .5D0*SUM
1060        END DO
1061        VECT1(1) = .5D0*(VECT2(1)**2-VECT2(2)**2)
1062        VECT1(2) = VECT2(1)*VECT2(2)
1063        VECT1(3) = VECT2(3)**2
1064        VECT1(4) = VECT2(4)**2
1065        TEMP = -1.D1*VECT2(1) - 1.D1*VECT2(2)
1066        VECT2(2) = 1.D1*VECT2(1) - 1.D1*VECT2(2)
1067        VECT2(1) = TEMP
1068        VECT2(3) = 1.D3*VECT2(3)
1069        VECT2(4) = 1.D-2*VECT2(4)
1070        SUM = 0.D0
1071        DO I = 1, 4
1072          SUM = SUM + VECT1(I) - VECT2(I)
1073        END DO
1074        DO I = 1, 4
1075          YP(I) = VECT2(I) - VECT1(I) + .5D0*SUM
1076        END DO
1077        GOTO 260
1078
1079!     PROBLEM E5
1080200     XTEMP = -7.89D-10*Y(1)
1081        YP(1) = XTEMP - 1.1D7*Y(1)*Y(3)
1082        YP(2) = -XTEMP - 1.13D9*Y(2)*Y(3)
1083        YP(4) = 1.1D7*Y(1)*Y(3) - 1.13D3*Y(4)
1084        YP(3) = YP(2) - YP(4)
1085        GOTO 260
1086
1087!     PROBLEM CLASS F - CHEMICAL KINETICS EQUATIONS
1088
1089!     PROBLEM F1
1090210     TEMP = 6.D-3*EXP(20.7D0-1.5D4/Y(1))
1091        YP(1) = 1.3D0*(Y(3)-Y(1)) + 1.04D4*TEMP*Y(2)
1092        YP(2) = 1.88D3*(Y(4)-Y(2)*(1.D0+TEMP))
1093        YP(3) = 1752.D0 - 269.D0*Y(3) + 267.D0*Y(1)
1094        YP(4) = .1D0 + 320.D0*Y(2) - 321.D0*Y(4)
1095        GOTO 260
1096
1097!     PROBLEM F2
1098220     YP(1) = -Y(1) - Y(1)*Y(2) + 294.D0*Y(2)
1099        YP(2) = Y(1)*(1.D0-Y(2))/98.D0 - 3.D0*Y(2)
1100        GOTO 260
1101
1102!     PROBLEM F3
1103230     YP(1) = -1.0D7*Y(2)*Y(1) + 1.D1*Y(3)
1104        YP(2) = -1.0D7*Y(2)*Y(1) - 1.D7*Y(2)*Y(5) + 1.D1*Y(3) + 1.D1*Y(4)
1105        YP(3) = 1.0D7*Y(2)*Y(1) - 1.001D4*Y(3) + 1.D-3*Y(4)
1106        YP(4) = 1.D4*Y(3) - 1.0001D1*Y(4) + 1.D7*Y(2)*Y(5)
1107        YP(5) = 1.D1*Y(4) - 1.D7*Y(2)*Y(5)
1108        GOTO 260
1109
1110!     PROBLEM F4
1111240     S = 77.27D0
1112        T = 0.161D0
1113        Q = 8.375D-6
1114        F = 1.D0
1115        YP(1) = S*(Y(2)-Y(1)*Y(2)+Y(1)-Q*Y(1)*Y(1))
1116        YP(2) = (-Y(2)-Y(1)*Y(2)+F*Y(3))/S
1117        YP(3) = T*(Y(1)-Y(3))
1118        GOTO 260
1119
1120!     PROBLEM F5
1121250     YP(1) = -3.D11*Y(1)*Y(2) + 1.2D8*Y(4) - 9.D11*Y(1)*Y(3)
1122        YP(2) = -3.D11*Y(1)*Y(2) + 2.D7*Y(4)
1123        YP(3) = -9.D11*Y(1)*Y(3) + 1.D8*Y(4)
1124        YP(4) = 3.D11*Y(1)*Y(2) - 1.2D8*Y(4) + 9.D11*Y(1)*Y(3)
1125260     CONTINUE
1126        IF (IWT<0) GOTO 270
1127        DO I = 1, N
1128          YP(I) = YP(I)/W(I)
1129          Y(I) = YTEMP(I)
1130        END DO
1131270     CONTINUE
1132        RETURN
1133      END SUBROUTINE FCN
1134
1135      SUBROUTINE PDERV(X,Y)
1136
1137!     ROUTINE TO EVALUATE THE JACOBIAN MATRIX OF PARTIAL DERIVATIVES
1138!     CORRESPONDING TO THE DIFFERENTIAL EQUATION:
1139!                   DY/DX = F(X,Y).
1140!     THE N**2 ELEMENTS OF THE ARRAY DY(*) ARE ASSIGNED THE VALUE OF
1141!     THE JACOBIAN MATRIX WITH ELEMENT I+(J-1)*N BEING ASSIGNED THE
1142!     VALUE OF DF(I)/DY(J). THE PARTICULAR EQUATION BEING INTEGRATED
1143!     IS INDICATED BY THE VALUE OF THE FLAG ID WHICH IS PASSED THROUGH
1144!     COMMON. IF A SCALED DIFFERENTIAL EQUATION IS BEING SOLVED (AS
1145!     SIGNALLED IWT) THE ELEMENTS OF THE JACOBIAN ARE SCALED ACCORDING-
1146!     LY BY THE WEIGHT VECTOR W(*).
1147
1148        IMPLICIT NONE
1149
1150!     .. Scalar Arguments ..
1151        DOUBLE PRECISION X
1152!     .. Array Arguments ..
1153        DOUBLE PRECISION Y(20)
1154!     .. Local Scalars ..
1155        DOUBLE PRECISION F, Q, S, SUM, T, TEMP, XTEMP1, XTEMP2, XTEMP3
1156        INTEGER I, ITMP, J, L
1157!     .. Local Arrays ..
1158        DOUBLE PRECISION BPARM(4), CPARM(4), VECT2(4), YTEMP(20)
1159!     .. Data statements ..
1160        DATA BPARM/3.D0, 8.D0, 25.D0, 1.D2/
1161        DATA CPARM/1.D-1, 1.D0, 1.D1, 2.D1/
1162
1163!     .. Executable Statements ..
1164
1165        NJAC = NJAC + 1
1166        IF (IWT<0) GOTO 10
1167        DO I = 1, N
1168          YTEMP(I) = Y(I)
1169          Y(I) = Y(I)*W(I)
1170        END DO
117110      CONTINUE
1172        GOTO (20,30,40,50,260,260,260,260,260,260,60,70,70,70,70,260,260,260, &
1173          260,260,80,90,90,90,90,260,260,260,260,260,100,110,120,130,140,150, &
1174          260,260,260,260,160,170,180,190,200,260,260,260,260,260,210,220,230, &
1175          240,250) ID
1176        GOTO 260
1177
1178!     PROBLEM CLASS A - LINEAR WITH REAL EIGENVALUES
1179
1180!     PROBLEM A1
118120      DO I = 1, 16
1182          DY(I) = 0.D0
1183        END DO
1184        DY(1) = -.5D0
1185        DY(6) = -1.D0
1186        DY(11) = -1.D2
1187        DY(16) = -9.D1
1188        GOTO 260
1189
1190!     PROBLEM A2
119130      DO I = 1, 81
1192          DY(I) = 0.D0
1193        END DO
1194        DO I = 2, 62, 10
1195          DY(I) = 1.D0
1196          DY(I+9) = -2.D0
1197          DY(I+18) = 1.D0
1198        END DO
1199        DY(1) = -1.8D3
1200        DY(10) = 9.D2
1201        DY(72) = 1.D3
1202        DY(81) = -2.D3
1203        GOTO 260
1204
1205!     PROBLEM A3
120640      DO I = 1, 16
1207          DY(I) = 0.D0
1208        END DO
1209        DY(1) = -1.D4
1210        DY(5) = 1.D2
1211        DY(6) = -1.D3
1212        DY(9) = -1.D1
1213        DY(10) = 1.D1
1214        DY(11) = -1.D0
1215        DY(13) = 1.D0
1216        DY(14) = -1.D1
1217        DY(15) = 1.D1
1218        DY(16) = -1.D-1
1219        GOTO 260
1220
1221!     PROBLEM A4
122250      DO I = 1, 100
1223          DY(I) = 0.D0
1224        END DO
1225        DO I = 1, 10
1226          DY((I-1)*10+I) = -REAL(I)**5
1227        END DO
1228        GOTO 260
1229
1230!     PROBLEM CLASS B - LINEAR WITH NON-REAL EIGENVALUES
1231
1232!     PROBLEM B1
123360      DO I = 1, 16
1234          DY(I) = 0.D0
1235        END DO
1236        DY(1) = -1.D0
1237        DY(2) = -1.D2
1238        DY(5) = 1.D0
1239        DY(6) = -1.D0
1240        DY(11) = -1.D2
1241        DY(12) = -1.D4
1242        DY(15) = 1.D0
1243        DY(16) = -1.D2
1244        GOTO 260
1245
1246!     PROBLEMS B2, B3, B4, B5
124770      DO I = 1, 36
1248          DY(I) = 0.D0
1249        END DO
1250        DY(1) = -1.D1
1251        DY(2) = -BPARM(IID-1)
1252        DY(7) = BPARM(IID-1)
1253        DY(8) = -1.D1
1254        DY(15) = -4.D0
1255        DY(22) = -1.D0
1256        DY(29) = -.5D0
1257        DY(36) = -.1D0
1258        GOTO 260
1259
1260!     PROBLEM CLASS C - NON-LINEAR COUPLING FROM
1261!                       STEADY STATE TO TRANSIENT
1262
1263!     PROBLEM C1
126480      DO I = 1, 16
1265          DY(I) = 0.D0
1266        END DO
1267        DY(1) = -1.D0
1268        DY(5) = 2.D0*Y(2)
1269        DY(6) = -1.D1
1270        DY(9) = 2.D0*Y(3)
1271        DY(10) = 2.D1*Y(3)
1272        DY(11) = -4.D1
1273        DY(13) = 2.D0*Y(4)
1274        DY(14) = 2.D1*Y(4)
1275        DY(15) = 8.D1*Y(4)
1276        DY(16) = -1.D2
1277        GOTO 260
1278
1279!     PROBLEMS C2, C3, C4, C5
128090      DO I = 1, 16
1281          DY(I) = 0.D0
1282        END DO
1283        DY(1) = -1.D0
1284        DY(2) = 2.D0*Y(1)*CPARM(IID-1)
1285        DY(3) = 8.D0*Y(1)*CPARM(IID-1)
1286        DY(4) = 2.D1*Y(1)*CPARM(IID-1)
1287        DY(6) = -1.D1
1288        DY(7) = 8.D0*Y(2)*CPARM(IID-1)
1289        DY(8) = 2.D1*Y(2)*CPARM(IID-1)
1290        DY(11) = -4.D1
1291        DY(12) = 2.D1*Y(3)*CPARM(IID-1)
1292        DY(16) = -1.D2
1293        GOTO 260
1294
1295!     PROBLEM CLASS D - NON-LINEAR WITH REAL EIGENVALUES
1296
1297!     PROBLEM D1
1298100     DY(1) = -.2D0
1299        DY(2) = 1.D1
1300        DY(3) = 0.D0
1301        DY(4) = .2D0
1302        DY(5) = -6.D1 + .125D0*Y(3)
1303        DY(6) = 0.D0
1304        DY(7) = 0.D0
1305        DY(8) = .125D0*Y(2) + .125D0
1306        DY(9) = 0.D0
1307        GOTO 260
1308
1309!     PROBLEM D2
1310110     DY(1) = -4.D-2
1311        DY(2) = 4.D2
1312        DY(3) = 0.D0
1313        DY(4) = 1.D-2*Y(3)
1314        DY(5) = -1.D2*Y(3) - 6.D3*Y(2)
1315        DY(6) = 6.D1*Y(2)
1316        DY(7) = .1D-1*Y(2)
1317        DY(8) = -1.D2*Y(2)
1318        DY(9) = 0.D0
1319        GOTO 260
1320
1321!     PROBLEM D3
1322120     DY(1) = -1.D2*Y(2)
1323        DY(2) = DY(1)
1324        DY(3) = -DY(1)
1325        DY(4) = 0.D0
1326        DY(5) = -1.D2*Y(1)
1327        DY(7) = -DY(5)
1328        DY(8) = 2.D4*Y(2)
1329        DY(6) = DY(5) - DY(8)
1330        DY(6) = DY(6) - 2.D4*Y(2)
1331        DY(9) = 1.D0
1332        DY(10) = 1.D0
1333        DY(11) = -1.D0
1334        DY(12) = 0.D0
1335        DY(13) = 0.D0
1336        DY(14) = 2.D0
1337        DY(15) = 0.D0
1338        DY(16) = -1.D0
1339        GOTO 260
1340
1341!     PROBLEM D4
1342130     DY(1) = -.013D0 - 1.D3*Y(3)
1343        DY(2) = 0.D0
1344        DY(4) = 0.D0
1345        DY(5) = -2.5D3*Y(3)
1346        DY(7) = -1.D3*Y(1)
1347        DY(8) = -2.5D3*Y(2)
1348        DO I = 3, 9, 3
1349          DY(I) = DY(I-1) + DY(I-2)
1350        END DO
1351        GOTO 260
1352
1353!     PROBLEM D5
1354140     XTEMP1 = Y(1) + 1.D3
1355        XTEMP2 = Y(1) + 1.D0
1356        XTEMP3 = .01D0 + Y(1) + Y(2)
1357        DY(2) = -(1.D0+Y(2)**2)
1358        DY(3) = -(1.D0+XTEMP1*XTEMP2)
1359        DY(1) = -(-DY(3)+XTEMP3*(XTEMP1+XTEMP2))
1360        DY(4) = -(2.D0*XTEMP3*Y(2)-DY(2))
1361        GOTO 260
1362
1363!     PROBLEM D6
1364150     DY(1) = -1.D0 - 1.D8*Y(3)
1365        DY(2) = 0.D0
1366        DY(4) = 0.D0
1367        DY(5) = -1.D1 - 3.D7*Y(3)
1368        DY(7) = 1.D8*(1.D0-Y(1))
1369        DY(8) = 3.D7*(1.D0-Y(2))
1370        DO I = 3, 9, 3
1371          DY(I) = -DY(I-2) - DY(I-1)
1372        END DO
1373        GOTO 260
1374
1375!     PROBLEM CLASS E - NON-LINEAR WITH NON-REAL EIGENVALUES
1376
1377!     PROBLEM E1
1378160     DO I = 1, 16
1379          DY(I) = 0.D0
1380        END DO
1381        DY(5) = 1.D0
1382        DY(10) = 1.D0
1383        DY(15) = 1.D0
1384        XTEMP1 = Y(1)
1385        XTEMP2 = Y(2)/(XTEMP1**2+1.D0)**2
1386        DY(4) = 3.D0*XTEMP1**2 - XTEMP1*COS(XTEMP1) - SIN(XTEMP1) - 1.D8 - &
1387          2.D0*XTEMP1*Y(2)*Y(3)*XTEMP2
1388        DY(8) = 2.D0*Y(3)*Y(2)/(1.D0+Y(1)**2) - 4.D6
1389        DY(12) = Y(2)*Y(2)/(1.D0+Y(1)**2) + 1.D0 - 6.D4
1390        DY(16) = 1.D1*EXP(-Y(4)**2)*(1.D0-2.D0*Y(4)**2) - 4.D2
1391        GOTO 260
1392
1393!     PROBLEM E2
1394170     DY(1) = 0.D0
1395        DY(2) = -1.D1*Y(1)*Y(2) - 1.D0
1396        DY(3) = 1.D0
1397        DY(4) = 5.D0 - 5.D0*Y(1)*Y(1)
1398        GOTO 260
1399
1400!     PROBLEM E3
1401180     DY(1) = -55.D0 - Y(3)
1402        DY(2) = .785D-1
1403        DY(3) = 0.1D0
1404        DY(4) = 65.D0
1405        DY(5) = -.785D-1
1406        DY(6) = 0.D0
1407        DY(7) = -Y(1)
1408        DY(8) = 0.D0
1409        DY(9) = 0.D0
1410        GOTO 260
1411
1412!     PROBLEM E4
1413190     SUM = Y(1) + Y(2) + Y(3) + Y(4)
1414        DO I = 1, 4
1415          VECT2(I) = -Y(I) + .5D0*SUM
1416        END DO
1417        DO I = 1, 16
1418          DY(I) = 0.D0
1419        END DO
1420        DY(1) = VECT2(1) + 1.D1
1421        DY(2) = VECT2(2) - 1.D1
1422        DY(5) = -DY(2)
1423        DY(6) = DY(1)
1424        DY(11) = 2.D0*VECT2(3) - 1.D3
1425        DY(16) = 2.D0*VECT2(4) - 1.D-2
1426        DO I = 1, 4
1427          SUM = 0.D0
1428          DO J = 1, 4
1429            L = I + (J-1)*4
1430            SUM = SUM + DY(L)
1431          END DO
1432          DO J = 1, 4
1433            L = I + (J-1)*4
1434            DY(L) = -DY(L) + .5D0*SUM
1435          END DO
1436        END DO
1437        DO J = 1, 4
1438          SUM = 0.D0
1439          DO I = 1, 4
1440            L = I + (J-1)*4
1441            SUM = SUM + DY(L)
1442          END DO
1443          DO I = 1, 4
1444            L = I + (J-1)*4
1445            DY(L) = -DY(L) + .5D0*SUM
1446          END DO
1447        END DO
1448        GOTO 260
1449
1450!     PROBLEM E5
1451200     DY(1) = -7.89D-10 - 1.1D7*Y(3)
1452        DY(2) = 7.89D-10
1453        DY(4) = 1.1D7*Y(3)
1454        DY(5) = 0.D0
1455        DY(6) = -1.13D9*Y(3)
1456        DY(8) = 0.D0
1457        DY(9) = -1.1D7*Y(1)
1458        DY(10) = -1.13D9*Y(2)
1459        DY(12) = -DY(9)
1460        DY(13) = 0.D0
1461        DY(14) = 0.D0
1462        DY(16) = -1.13D3
1463        DO I = 3, 15, 4
1464          DY(I) = DY(I-1) - DY(I+1)
1465        END DO
1466        GOTO 260
1467
1468!     PROBLEM CLASS F - CHEMICAL KINETICS EQUATIONS
1469
1470!     PROBLEM F1
1471210     TEMP = 90.D0*EXP(20.7D0-1.5D4/Y(1))/Y(1)**2
1472        DY(1) = -1.3D0 + 1.04D4*TEMP*Y(2)
1473        DY(2) = -1.88D3*Y(2)*TEMP
1474        DY(3) = 267.D0
1475        DY(4) = 0.D0
1476        TEMP = 6.D-3*EXP(20.7D0-1.5D4/Y(1))
1477        DY(5) = 1.04D4*TEMP
1478        DY(6) = -1.88D3*(1.D0+TEMP)
1479        DY(7) = 0.D0
1480        DY(8) = 320.D0
1481        DY(9) = 1.3D0
1482        DY(10) = 0.D0
1483        DY(11) = -269.D0
1484        DY(12) = 0.0D0
1485        DY(13) = 0.0D0
1486        DY(14) = 1.88D3
1487        DY(15) = 0.0D0
1488        DY(16) = -321.0D0
1489        GOTO 260
1490
1491!     PROBLEM F2
1492220     DY(1) = -1.D0 - Y(2)
1493        DY(2) = (1.D0-Y(2))/98.D0
1494        DY(3) = -Y(1) + 294.D0
1495        DY(4) = -Y(1)/98.D0 - 3.D0
1496        GOTO 260
1497
1498!     PROBLEM F3
1499230     DY(1) = -1.D7*Y(2)
1500        DY(2) = -1.D7*Y(2)
1501        DY(3) = 1.D7*Y(2)
1502        DY(4) = 0.0D0
1503        DY(5) = 0.0D0
1504        DY(6) = -1.D7*Y(1)
1505        DY(7) = -1.D7*Y(1) - 1.D7*Y(5)
1506        DY(8) = 1.D7*Y(1)
1507        DY(9) = 1.D7*Y(5)
1508        DY(10) = -1.D7*Y(5)
1509        DY(11) = 1.D1
1510        DY(12) = 1.D1
1511        DY(13) = -1.001D4
1512        DY(14) = 1.D4
1513        DY(15) = 0.0D0
1514        DY(16) = 0.0D0
1515        DY(17) = 1.D1
1516        DY(18) = 1.D-3
1517        DY(19) = -1.0001D1
1518        DY(20) = 1.D1
1519        DY(21) = 0.0D0
1520        DY(22) = -1.D7*Y(2)
1521        DY(23) = 0.0D0
1522        DY(24) = 1.D7*Y(2)
1523        DY(25) = -1.0D7*Y(2)
1524        GOTO 260
1525
1526!     PROBLEM F4
1527240     S = 77.27D0
1528        T = 0.161D0
1529        Q = 8.375D-6
1530        F = 1.D0
1531        DY(1) = S*(-Y(2)+1.D0-2.D0*Q*Y(1))
1532        DY(2) = -Y(2)/S
1533        DY(3) = T
1534        DY(4) = S*(1.D0-Y(1))
1535        DY(5) = (-1.D0-Y(1))/S
1536        DY(6) = 0.D0
1537        DY(7) = 0.D0
1538        DY(8) = F/S
1539        DY(9) = -T
1540        GOTO 260
1541
1542!     PROBLEM F5
1543250     DY(1) = -3.D11*Y(2) - 9.D11*Y(3)
1544        DY(2) = -3.D11*Y(2)
1545        DY(3) = -9.D11*Y(3)
1546        DY(4) = 3.D11*Y(2) + 9.D11*Y(3)
1547        DY(5) = -3.D11*Y(1)
1548        DY(6) = -3.D11*Y(1)
1549        DY(7) = 0.0D0
1550        DY(8) = 3.D11*Y(1)
1551        DY(9) = -9.D11*Y(1)
1552        DY(10) = 0.0D0
1553        DY(11) = -9.D11*Y(1)
1554        DY(12) = 9.D11*Y(1)
1555        DY(13) = 1.2D8
1556        DY(14) = 2.D7
1557        DY(15) = 1.D8
1558        DY(16) = -1.2D8
1559260     CONTINUE
1560        IF (IWT<0) GOTO 270
1561        DO I = 1, N
1562          Y(I) = YTEMP(I)
1563          DO J = 1, N
1564            ITMP = I + (J-1)*N
1565            DY(ITMP) = DY(ITMP)*W(J)/W(I)
1566          END DO
1567        END DO
1568270     CONTINUE
1569        RETURN
1570      END SUBROUTINE PDERV
1571
1572    END MODULE STIFFSET
1573
1574!******************************************************************
1575
1576    PROGRAM DEMOSTIFF
1577
1578      USE STIFFSET
1579      USE DVODE_F90_M
1580
1581      IMPLICIT NONE
1582      INTEGER ITASK, ISTATE, ISTATS, NEQ, I, CLASS, PROBLEM, MYID, ITEST, ISTR, &
1583        IJAC, JTYPE, CONSTANTJ, ITOL, IADIM, IA, JADIM, JA, ROW, COL, ML, MU,   &
1584        IJACSP , NUM_S, NUM_J, NUM_D, Total_Solutions, ISCALE, ITIME, COMPS,    &
1585        METHOD, WHICH_PROBS, WHICH_TOLS, COUNTS
1586      DOUBLE PRECISION RSTATS, T, TOUT, HBEGIN, HBOUND, TBEGIN, TEND, Y, EPS,   &
1587        YINIT, YFINAL, RELERR_TOLERANCES, ABSERR_TOLERANCES, AERROR, ERRMAX,    &
1588        ERSTATS, LBOUNDS, UBOUNDS, ERSTATS2
1589      REAL DVTIME, DVTIME1, DVTIME2, EXECUTION_TIMES 
1590      DIMENSION Y(20), RSTATS(22), ISTATS(31), YINIT(20), YFINAL(20), MYID(55)
1591      DIMENSION RELERR_TOLERANCES(20), ABSERR_TOLERANCES(20), AERROR(20),       &
1592        IA(21), JA(400), ERSTATS(12,55,3), NUM_S(12,55,3), NUM_J(12,55,3),      &
1593        NUM_D(12,55,3), EXECUTION_TIMES(12), COMPS(20), LBOUNDS(20),            &
1594        UBOUNDS(20), ERSTATS2(2,12,55,6), WHICH_PROBS(55), WHICH_TOLS(12),      &
1595        COUNTS(2,12,55,6,3)
1596      LOGICAL J_IS_CONSTANT, ANALYTIC_JACOBIAN, ERRORS, SUPPLY_STRUCTURE,       &
1597        TOLVEC, USE_JACSP, BOUNDS, USEW, USEHBEGIN
1598      INTRINSIC ABS, MAX, MAXVAL, EPSILON
1599      TYPE (VODE_OPTS) :: OPTIONS
1600      DATA MYID/1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 11, 12, 13, 14, 15,  0, 0, 0, 0,  &
1601        0, 21, 22, 23, 24, 25, 0, 0, 0, 0, 0, 31, 32, 33, 34, 35, 36, 0, 0, 0,  &
1602        0, 41, 42, 43, 44, 45, 0, 0, 0, 0, 0, 51, 52, 53, 54, 55/
1603
1604      OPEN (UNIT=6,FILE='stiffoptions.dat')
1605
1606!     Initialize the cumulative solutions total.
1607      Total_Solutions = 0
1608
1609!     Use vector tolerances (same solution in either case)?
1610      TOLVEC = .FALSE.
1611     
1612!     Initialize the flag to indicate whether any errors occurred.
1613      ERRORS = .FALSE.
1614
1615!     Execution times per tolerance.
1616      EXECUTION_TIMES(1:12) = 0.0E0
1617      ERSTATS2(1:2,1:12,1:55,1:6) = 0.0D0
1618      WHICH_PROBS(1:55) = 0
1619      WHICH_TOLS(1:12) = 0
1620      COUNTS(1:2,1:12,1:55,1:6,1:3) = 0
1621
1622      WRITE(6,90032)
1623
1624!     Enter the scale equations loop.
1625!        1 = the ODEs will be scaled
1626!        2 = the ODEs will not be scaled
1627      DO ISCALE = 1 , 2
1628      USEW = .TRUE.
1629      IF (ISCALE == 2) USEW = .FALSE.
1630
1631!     Enter the numerical Jacobian formation loop.
1632!        1 = use Doug's JACSP for Jacobian
1633!        2 = use original VODE Jacobian algorithm
1634      DO IJACSP = 1, 2
1635      USE_JACSP = .TRUE.
1636      IF (IJACSP == 2) USE_JACSP = .FALSE. 
1637
1638!     Initialize the max error array for this value of IJACSP.
1639      ERSTATS(1:12,1:15,1:3) = 0.0D0
1640
1641!     Initialize the stats arrays.
1642      NUM_S(1:12,1:15,1:3) = 0
1643      NUM_J(1:12,1:15,1:3) = 0
1644      NUM_D(1:12,1:15,1:3) = 0
1645
1646!     Enter the sparsity structure loop.
1647!        1 = numerically determined sparse pointer arrays
1648!        2 = supplied pointer arrays
1649!     Make this a one-trip loop (1 to 1 or 2 to 2) for stats to be meaningful:
1650      DO ISTR = 2, 2
1651      SUPPLY_STRUCTURE = .FALSE.
1652      IF (ISTR == 2) SUPPLY_STRUCTURE = .TRUE.
1653
1654!     Enter the error tolerance loop.
1655!        Error tolerance = 10^(-ITOL)
1656      DO ITOL = 5, 12, 1
1657
1658      WHICH_TOLS(ITOL) = 1
1659
1660!     Enter the linear algebra type loop.
1661!        1 = dense Jacobian
1662!        2 = sparse Jacobian
1663!        3 = banded Jacobian
1664      DO JTYPE = 1, 3
1665
1666!     Enter the constant Jacobian loop.
1667!        1 = Jacobian is not constant
1668!        2 = Jacobian is constant for classes 0 and 1
1669!     Make this a one-trip loop (1 to 1 or 2 to 2) for stats to be meaningful:
1670      DO CONSTANTJ = 1, 1
1671
1672!     Enter the numerical or analytical Jacobian loop.
1673!        1 = numerical Jacobian
1674!        2 = analytical Jacobian
1675!     Make this a one-trip loop (1 to 1 or 2 to 2) for stats to be meaningful:
1676      DO IJAC = 1, 1
1677      ANALYTIC_JACOBIAN = .TRUE.
1678      IF (IJAC == 1) ANALYTIC_JACOBIAN = .FALSE.
1679
1680      IF (JTYPE==1 .AND. IJACSP==1) METHOD = 1
1681      IF (JTYPE==2 .AND. IJACSP==1) METHOD = 2
1682      IF (JTYPE==3 .AND. IJACSP==1) METHOD = 3
1683      IF (JTYPE==1 .AND. IJACSP==2) METHOD = 4
1684      IF (JTYPE==2 .AND. IJACSP==2) METHOD = 5
1685      IF (JTYPE==3 .AND. IJACSP==2) METHOD = 6
1686
1687!     Enter the problem test loop.
1688      DO ITEST = 1, 55
1689
1690!       JACS supplies the full matrix in the analytic Jacobian case.
1691        IF (IJAC == 2 .AND. ISTR == 1) GOTO 20
1692        ID = MYID(ITEST)
1693        IF (ID == 0) GOTO 20
1694        WHICH_PROBS(ITEST) = 1
1695        IID = MOD(ID,10)
1696        WRITE (6,90010)
1697        CLASS = ID/10
1698        PROBLEM = ID - 10*CLASS
1699        J_IS_CONSTANT = .FALSE.
1700        IF ((CLASS==0 .OR. CLASS==1) .AND. CONSTANTJ==2) J_IS_CONSTANT = .TRUE.
1701
1702!       Output header.
1703        IF (CLASS == 0) THEN
1704          WRITE (6,90000) PROBLEM
1705        ELSE IF (CLASS == 1) THEN
1706          WRITE (6,90001) PROBLEM
1707        ELSE IF (CLASS == 2) THEN
1708          WRITE (6,90002) PROBLEM
1709        ELSE IF (CLASS == 3) THEN
1710          WRITE (6,90003) PROBLEM
1711        ELSE IF (CLASS == 4) THEN
1712          WRITE (6,90004) PROBLEM
1713        ELSE IF (CLASS == 5) THEN
1714          WRITE (6,90005) PROBLEM
1715        END IF
1716        WRITE(6,*)    ' Results for ITOL = ', ITOL, ' follow.'
1717        PRINT *,      ' Results for ITOL = ', ITOL, ' follow.'
1718        IF (IJACSP == 1) THEN
1719           WRITE(6,*) ' Results for JACSP Jacobian follow.'
1720           PRINT *,   ' Results for JACSP Jacobian follow.'
1721        ELSE
1722           WRITE(6,*) ' Results for VODE Jacobian follow.'
1723           PRINT *,   ' Results for VODE Jacobian follow.'
1724        END IF
1725        IF (ISCALE == 1) THEN
1726           WRITE(6,*) ' The ODEs will be scaled.'
1727           PRINT *,   ' The ODEs will be scaled.'
1728        ELSE
1729           WRITE(6,*) ' The ODEs will not be scaled.'
1730           PRINT *,   ' The ODEs will not be scaled.'
1731        END IF
1732        IF (ISTR == 1) THEN
1733           WRITE(6,*) ' Results for numerically determined sparse pointers follow.'
1734           PRINT *,   ' Results for numerically determined sparse pointers follow.'
1735        ELSE
1736           WRITE(6,*) ' Results for user supplied sparse pointers follow.'
1737           PRINT *,   ' Results for user supplied sparse pointers follow.'
1738        END IF
1739        IF (JTYPE == 1) THEN
1740           WRITE(6,*) ' Results for dense Jacobian option follow.'
1741           PRINT *,   ' Results for dense Jacobian option follow.'
1742        END IF
1743        IF (JTYPE == 2) THEN
1744           WRITE(6,*) ' Results for sparse Jacobian option follow.'
1745           PRINT *,   ' Results for sparse Jacobian option follow.'
1746        END IF
1747        IF (JTYPE == 3) THEN
1748           WRITE(6,*) ' Results for banded Jacobian option follow.'
1749           PRINT *,   ' Results for banded Jacobian option follow.'
1750        END IF
1751        IF (CONSTANTJ == 1) THEN
1752           WRITE(6,*) ' Results for non-constant Jacobian option follow.'
1753           PRINT *,   ' Results for non-constant Jacobian option follow.'
1754        ELSE
1755           WRITE(6,*) ' Results for constant Jacobian option follow.'
1756           PRINT *,   ' Results for constant Jacobian option follow.'
1757        END IF
1758        IF (IJAC == 1) THEN
1759           WRITE(6,*) ' Results for numerical Jacobian option follow.'
1760           PRINT *,   ' Results for numerical Jacobian option follow.'
1761        ELSE
1762           WRITE(6,*) ' Results for analytical Jacobian option follow.'
1763           PRINT *,   ' Results for analytical Jacobian option follow.'
1764        END IF
1765
1766!       Initialize the problem.
1767!       Scale the ODEs?
1768!         USEW = .TRUE.
1769          IWT = -1
1770          IF (USEW) IWT = 1
1771!       Get the problem information.
1772          CALL IVALU(TBEGIN,TEND,HBEGIN,HBOUND,YINIT)
1773!       Use the IVALU starting step size?
1774          USEHBEGIN = .TRUE.
1775          IF (.NOT. USEHBEGIN) HBEGIN = 0.0D0
1776!       Define the number of ODEs.
1777          NEQ = N
1778!       Use the full matrix for the banded solution.
1779          ML = NEQ - 1
1780          MU = NEQ - 1
1781!       Initial and final integration times.
1782          T = TBEGIN
1783          TOUT = TEND
1784!       Initial solution.
1785          Y(1:NEQ) = YINIT(1:NEQ)
1786!       Error tolerances.
1787          EPS = 1.0D0 / 10.0D0**ITOL
1788          RELERR_TOLERANCES(1:NEQ) = EPS
1789          ABSERR_TOLERANCES(1:NEQ) = EPS
1790          WRITE (6,90007) ID,TBEGIN,TEND,HBEGIN,HBOUND,IWT,N,EPS,Y(1:NEQ)
1791          BOUNDS = .FALSE.
1792!         Force nonnegativity for the chemical kinetics problems in Class F.
1793          IF (ITEST >= 51) THEN
1794             BOUNDS = .TRUE.
1795             COMPS(1) = 1
1796             COMPS(2) = 2
1797             COMPS(3) = 3
1798             DO I = 1, NEQ
1799                COMPS(I) = I
1800             END DO
1801             LBOUNDS(1:NEQ) = 0.0D0
1802             UBOUNDS(1:NEQ) = 1.0D50
1803          END IF
1804        IF (BOUNDS) WRITE(6,90040)
1805!       Initialize the integration flags.
1806          ITASK = 1
1807          ISTATE = 1
1808
1809!       Use the full matrix for sparsity structure.
1810        IF (JTYPE == 2 .AND. SUPPLY_STRUCTURE) THEN
1811           IADIM = NEQ + 1
1812           JADIM = NEQ * NEQ
1813           IA(1) = 1
1814           DO I = 1, NEQ
1815             IA(I+1) = IA(I) + NEQ
1816           END DO
1817           I = 0
1818           DO COL = 1, NEQ
1819             I = NEQ*(COL-1)
1820             DO ROW = 1, NEQ
1821               I = I + 1
1822               JA(I) = ROW
1823             END DO
1824           END DO
1825        END IF
1826
1827        CALL CPU_TIME(DVTIME1)
1828
1829!       Set the integration options.
1830        IF (JTYPE == 1) THEN
1831!          Dense solution ...
1832           IF (TOLVEC) THEN
1833!             Supply vectors of tolerances.
1834              IF (BOUNDS) THEN
1835                 OPTIONS = SET_OPTS(DENSE_J=.TRUE.,                          &
1836                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1837                                 RELERR_VECTOR=RELERR_TOLERANCES(1:NEQ),     &
1838                                 ABSERR_VECTOR=ABSERR_TOLERANCES(1:NEQ),     &
1839                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1840                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1841                                 JACOBIAN_BY_JACSP=USE_JACSP,                &
1842                                 CONSTRAINED=COMPS(1:NEQ),                   &
1843                                 CLOWER=LBOUNDS(1:NEQ),CUPPER=UBOUNDS(1:NEQ))
1844              ELSE
1845                 OPTIONS = SET_OPTS(DENSE_J=.TRUE.,                          &
1846                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1847                                 RELERR_VECTOR=RELERR_TOLERANCES(1:NEQ),     &
1848                                 ABSERR_VECTOR=ABSERR_TOLERANCES(1:NEQ),     &
1849                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1850                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1851                                 JACOBIAN_BY_JACSP=USE_JACSP)
1852              END IF
1853           ELSE
1854!             Supply scalar tolerances.
1855              IF (BOUNDS) THEN
1856                 OPTIONS = SET_OPTS(DENSE_J=.TRUE.,                          &
1857                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1858                                 RELERR=EPS,ABSERR=EPS,                      &
1859                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1860                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1861                                 JACOBIAN_BY_JACSP=USE_JACSP,                &
1862                                 CONSTRAINED=COMPS(1:NEQ),                   &
1863                                 CLOWER=LBOUNDS(1:NEQ),CUPPER=UBOUNDS(1:NEQ))
1864              ELSE
1865                 OPTIONS = SET_OPTS(DENSE_J=.TRUE.,                          &
1866                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1867                                 RELERR=EPS,ABSERR=EPS,                      &
1868                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1869                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1870                                 JACOBIAN_BY_JACSP=USE_JACSP)
1871              END IF
1872           END IF
1873        END IF
1874
1875        IF (JTYPE == 2) THEN
1876!          Sparse solution ...
1877           IF (TOLVEC) THEN
1878!             Supply vectors of tolerances.
1879              IF (BOUNDS) THEN
1880                 OPTIONS = SET_OPTS(SPARSE_J=.TRUE.,                         &
1881                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1882                                 RELERR_VECTOR=RELERR_TOLERANCES(1:NEQ),     &
1883                                 ABSERR_VECTOR=ABSERR_TOLERANCES(1:NEQ),     &
1884                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1885                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1886                                 USER_SUPPLIED_SPARSITY=SUPPLY_STRUCTURE,    &
1887                                 JACOBIAN_BY_JACSP=USE_JACSP,                &
1888                                 MA28_RPS=.TRUE.,CONSTRAINED=COMPS(1:NEQ),   &
1889                                 CLOWER=LBOUNDS(1:NEQ),CUPPER=UBOUNDS(1:NEQ))
1890              ELSE
1891                 OPTIONS = SET_OPTS(SPARSE_J=.TRUE.,                         &
1892                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1893                                 RELERR_VECTOR=RELERR_TOLERANCES(1:NEQ),     &
1894                                 ABSERR_VECTOR=ABSERR_TOLERANCES(1:NEQ),     &
1895                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1896                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1897                                 USER_SUPPLIED_SPARSITY=SUPPLY_STRUCTURE,    &
1898                                 JACOBIAN_BY_JACSP=USE_JACSP,                &
1899                                 MA28_RPS=.TRUE.)
1900              END IF
1901           ELSE
1902!             Supply scalar tolerances.
1903              IF (BOUNDS) THEN
1904                 OPTIONS = SET_OPTS(SPARSE_J=.TRUE.,                         &
1905                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1906                                 RELERR=EPS,ABSERR=EPS,                      &
1907                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1908                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1909                                 USER_SUPPLIED_SPARSITY=SUPPLY_STRUCTURE,    &
1910                                 JACOBIAN_BY_JACSP=USE_JACSP,                &
1911                                 MA28_RPS=.TRUE.,CONSTRAINED=COMPS(1:NEQ),   &
1912                                 CLOWER=LBOUNDS(1:NEQ),CUPPER=UBOUNDS(1:NEQ))
1913              ELSE
1914                 OPTIONS = SET_OPTS(SPARSE_J=.TRUE.,                         &
1915                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1916                                 RELERR=EPS,ABSERR=EPS,                      &
1917                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1918                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1919                                 USER_SUPPLIED_SPARSITY=SUPPLY_STRUCTURE,    &
1920                                 JACOBIAN_BY_JACSP=USE_JACSP,                &
1921                                 MA28_RPS=.TRUE.)
1922              END IF
1923           END IF
1924        END IF
1925
1926        IF (JTYPE == 3) THEN
1927!          Banded solution ...
1928           IF (TOLVEC) THEN
1929!             Supply vectors of tolerances.
1930              IF (BOUNDS) THEN
1931                 OPTIONS = SET_OPTS(BANDED_J=.TRUE.,                         &
1932                                 LOWER_BANDWIDTH=ML,UPPER_BANDWIDTH=MU,      &
1933                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1934                                 RELERR_VECTOR=RELERR_TOLERANCES(1:NEQ),     &
1935                                 ABSERR_VECTOR=ABSERR_TOLERANCES(1:NEQ),     &
1936                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1937                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1938                                 JACOBIAN_BY_JACSP=USE_JACSP,                &
1939                                 CONSTRAINED=COMPS(1:NEQ),                   &
1940                                 CLOWER=LBOUNDS(1:NEQ),CUPPER=UBOUNDS(1:NEQ))
1941              ELSE
1942                 OPTIONS = SET_OPTS(BANDED_J=.TRUE.,                         &
1943                                 LOWER_BANDWIDTH=ML,UPPER_BANDWIDTH=MU,      &
1944                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1945                                 RELERR_VECTOR=RELERR_TOLERANCES(1:NEQ),     &
1946                                 ABSERR_VECTOR=ABSERR_TOLERANCES(1:NEQ),     &
1947                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1948                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1949                                 JACOBIAN_BY_JACSP=USE_JACSP)
1950              END IF
1951           ELSE
1952!             Supply scalar tolerances.
1953              IF (BOUNDS) THEN
1954                 OPTIONS = SET_OPTS(BANDED_J=.TRUE.,                         &
1955                                 LOWER_BANDWIDTH=ML,UPPER_BANDWIDTH=MU,      &
1956                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1957                                 RELERR=EPS,ABSERR=EPS,                      &
1958                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1959                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1960                                 JACOBIAN_BY_JACSP=USE_JACSP,                &
1961                                 CONSTRAINED=COMPS(1:NEQ),                   &
1962                                 CLOWER=LBOUNDS(1:NEQ),CUPPER=UBOUNDS(1:NEQ))
1963              ELSE
1964                 OPTIONS = SET_OPTS(BANDED_J=.TRUE.,                         &
1965                                 LOWER_BANDWIDTH=ML,UPPER_BANDWIDTH=MU,      &
1966                                 USER_SUPPLIED_JACOBIAN=ANALYTIC_JACOBIAN,   &
1967                                 RELERR=EPS,ABSERR=EPS,                      &
1968                                 MXSTEP=100000,H0=HBEGIN,HMAX=HBOUND,        &
1969                                 CONSTANT_JACOBIAN=J_IS_CONSTANT,            &
1970                                 JACOBIAN_BY_JACSP=USE_JACSP)
1971              END IF
1972           END IF
1973        END IF
1974
1975!       Supply the sparse structure arrays directly.
1976        IF (JTYPE==2 .AND. SUPPLY_STRUCTURE) THEN
1977           CALL USERSETS_IAJA(IA(1:IADIM),IADIM,JA(1:JADIM),JADIM)
1978        END IF
1979
1980!       Perform the integration.
1981        IF (JTYPE == 1) THEN
1982           CALL DVODE_F90(DERIVS,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JACD)
1983        END IF
1984        IF (JTYPE == 2) THEN
1985           CALL DVODE_F90(DERIVS,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JACS)
1986        END IF
1987        IF (JTYPE == 3) THEN
1988           CALL DVODE_F90(DERIVS,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JACB)
1989        END IF
1990
1991        CALL CPU_TIME(DVTIME2)
1992        DVTIME = DVTIME2 - DVTIME1
1993        EXECUTION_TIMES(ITOL) = EXECUTION_TIMES(ITOL) + DVTIME
1994
1995        Total_Solutions = Total_Solutions + 1
1996
1997!       Check if an error occurred.
1998        IF (ISTATE < 0) THEN
1999          WRITE (6,90006) ISTATE
2000          ERRORS = .TRUE.
2001        END IF
2002
2003!       Gather the integration statistics.
2004        CALL GET_STATS(RSTATS,ISTATS)
2005        WRITE (6,90009) ISTATS(11), ISTATS(12), ISTATS(13)
2006        NUM_S(ITOL,ITEST,JTYPE) = ISTATS(11)
2007        NUM_D(ITOL,ITEST,JTYPE) = ISTATS(12)
2008        NUM_J(ITOL,ITEST,JTYPE) = ISTATS(13)
2009
2010!       Calculate the error in the final solution.
2011        CALL EVALU(YFINAL)
2012        DO I = 1, NEQ
2013          AERROR(I) = ABS(Y(I)-YFINAL(I))
2014        END DO
2015        ERRMAX = MAXVAL(AERROR(1:NEQ))
2016        ERSTATS(ITOL,ITEST,JTYPE) = ERRMAX
2017        ERSTATS2(ISCALE,ITOL,ITEST,METHOD) = ERRMAX
2018        COUNTS(ISCALE,ITOL,ITEST,METHOD,1) = ISTATS(11)
2019        COUNTS(ISCALE,ITOL,ITEST,METHOD,2) = ISTATS(12)
2020        COUNTS(ISCALE,ITOL,ITEST,METHOD,3) = ISTATS(13)
2021        WRITE (6,90008) (I,Y(I),YFINAL(I),AERROR(I),I=1,NEQ)
2022
2023!       Release the work arrays and determine how much storage was required.
2024        CALL RELEASE_ARRAYS
2025
202620    CONTINUE
2027      END DO ! End of ITEST Loop
2028      END DO ! End of IJAC Loop
2029      END DO ! End of JTYPE Loop
2030      END DO ! End of CONSTANTJ Loop
2031      END DO ! End of ITOL Loop
2032      END DO ! End of ISTR Loop
2033
2034!     Print the maximum errors for this value of IJACSP.
2035      DO ITOL = 5, 12, 1
2036         WRITE(6,90016) ITOL
2037         WRITE(6,90020)
2038         IF (IJACSP == 1) WRITE(6,90023)
2039         IF (IJACSP == 2) WRITE(6,90022)
2040         IF (ISCALE == 1) WRITE(6,90030)
2041         IF (ISCALE == 2) WRITE(6,90031)
2042         WRITE(6,90021)
2043         WRITE(6,90015)
2044         DO ITEST = 1, 55
2045            IF (WHICH_PROBS(ITEST) == 1) THEN
2046               ID = MYID(ITEST)
2047               IF (ID /= 0) WRITE(6,90014) ITOL, ITEST,         &
2048                            (ERSTATS(ITOL,ITEST,JTYPE),JTYPE=1,3)
2049            END IF
2050         END DO
2051      END DO 
2052
2053!     Print the integration stats for this value of IJACSP.
2054      DO ITOL = 5, 12, 1
2055         WRITE(6,90016) ITOL
2056         WRITE(6,90024)
2057         IF (IJACSP == 1) WRITE(6,90023)
2058         IF (IJACSP == 2) WRITE(6,90022)
2059         IF (ISCALE == 1) WRITE(6,90030)
2060         IF (ISCALE == 2) WRITE(6,90031)
2061         WRITE(6,90029)
2062         DO ITEST = 1, 55
2063            IF (WHICH_PROBS(ITEST) == 1) THEN
2064               ID = MYID(ITEST)
2065               IF (ID /= 0) THEN
2066                  WRITE(6,90026) ITEST, (NUM_S(ITOL,ITEST,JTYPE),JTYPE=1,3)
2067                  WRITE(6,90027)        (NUM_J(ITOL,ITEST,JTYPE),JTYPE=1,3)
2068                  WRITE(6,90028)        (NUM_D(ITOL,ITEST,JTYPE),JTYPE=1,3)
2069                  WRITE(6,*) ' '
2070               END IF
2071            END IF
2072         END DO
2073      END DO
2074
2075      END DO ! End of IJACSP Loop
2076      END DO ! End of ISCALE Loop
2077
2078      WRITE(6,90045)
2079
2080      DO ITEST = 1, 55
2081         IF (WHICH_PROBS(ITEST) == 1) THEN
2082            ID = MYID(ITEST)
2083            IID = MOD(ID,10)
2084            WRITE (6,90010)
2085            CLASS = ID/10
2086            PROBLEM = ID - 10*CLASS
2087            IF (CLASS == 0) THEN
2088              WRITE (6,90000) PROBLEM
2089            ELSE IF (CLASS == 1) THEN
2090              WRITE (6,90001) PROBLEM
2091            ELSE IF (CLASS == 2) THEN
2092              WRITE (6,90002) PROBLEM
2093            ELSE IF (CLASS == 3) THEN
2094              WRITE (6,90003) PROBLEM
2095            ELSE IF (CLASS == 4) THEN
2096              WRITE (6,90004) PROBLEM
2097            ELSE IF (CLASS == 5) THEN
2098              WRITE (6,90005) PROBLEM
2099            END IF
2100            DO ITOL = 5, 12, 1
2101               IF (WHICH_TOLS(ITOL) == 1) THEN
2102                  DO ISCALE = 1,  2
2103                     IF (ISCALE == 1) WRITE(6,90035) ITEST, ITOL
2104                     IF (ISCALE == 2) WRITE(6,90036) ITEST, ITOL
2105                     WRITE(6,90037) ERSTATS2(ISCALE,ITOL,ITEST,1), &
2106                                    ERSTATS2(ISCALE,ITOL,ITEST,4)
2107                     WRITE(6,90038) ERSTATS2(ISCALE,ITOL,ITEST,2), &
2108                                    ERSTATS2(ISCALE,ITOL,ITEST,5)
2109                     WRITE(6,90039) ERSTATS2(ISCALE,ITOL,ITEST,3), &
2110                                    ERSTATS2(ISCALE,ITOL,ITEST,6)
2111                     WRITE(6,90044)
2112                     WRITE(6,90041) COUNTS(ISCALE,ITOL,ITEST,1,1), &
2113                                    COUNTS(ISCALE,ITOL,ITEST,4,1), &
2114                                    COUNTS(ISCALE,ITOL,ITEST,1,2), &
2115                                    COUNTS(ISCALE,ITOL,ITEST,4,2), &
2116                                    COUNTS(ISCALE,ITOL,ITEST,1,3), &
2117                                    COUNTS(ISCALE,ITOL,ITEST,4,3)
2118                     WRITE(6,90042) COUNTS(ISCALE,ITOL,ITEST,2,1), &
2119                                    COUNTS(ISCALE,ITOL,ITEST,5,1), &
2120                                    COUNTS(ISCALE,ITOL,ITEST,2,2), &
2121                                    COUNTS(ISCALE,ITOL,ITEST,5,2), &
2122                                    COUNTS(ISCALE,ITOL,ITEST,2,3), &
2123                                    COUNTS(ISCALE,ITOL,ITEST,5,3)
2124                     WRITE(6,90043) COUNTS(ISCALE,ITOL,ITEST,3,1), &
2125                                    COUNTS(ISCALE,ITOL,ITEST,6,1), &
2126                                    COUNTS(ISCALE,ITOL,ITEST,3,2), &
2127                                    COUNTS(ISCALE,ITOL,ITEST,6,2), &
2128                                    COUNTS(ISCALE,ITOL,ITEST,3,3), &
2129                                    COUNTS(ISCALE,ITOL,ITEST,6,3)
2130                     WRITE (6,90010)
2131                  END DO
2132               END IF
2133            END DO
2134         END IF
2135      END DO
2136
2137!     Total number of solutions performed.
2138      WRITE(6,90025) Total_Solutions
2139      WRITE (6,90010)
2140
2141!     Execution Times.
2142      WRITE(6,90033)
2143      DO ITOL = 5, 12, 1
2144         ITIME = INT(EXECUTION_TIMES(ITOL))
2145         WRITE(6,90034) ITOL, ITIME
2146      END DO
2147      WRITE (6,90010)
2148
2149!     Where to search in output file.
2150      WRITE(6,90032)
2151
2152!     Indicate whether any DVODE_F90 made any error returns.
2153      IF (ERRORS) THEN
2154         WRITE(6,90011)
2155      ELSE
2156         WRITE(6,90012)
2157      END IF
2158
2159!     Format statements follow.
216090000 FORMAT (' Class/Problem = A',I1)
216190001 FORMAT (' Class/Problem = B',I1)
216290002 FORMAT (' Class/Problem = C',I1)
216390003 FORMAT (' Class/Problem = D',I1)
216490004 FORMAT (' Class/Problem = E',I1)
216590005 FORMAT (' Class/Problem = F',I1)
216690006 FORMAT (' An error occurred in VODE_F90. ISTATE = ',I3)
216790007 FORMAT (' Problem ID       = ',I3,/,   ' Initial time     = ',D15.5,/, &
2168              ' Final time       = ',D15.5,/,' Initial stepsize = ',D15.5,/, &
2169              ' Maximum stepsize = ',D15.5,/,' IWT flag         = ',I3,/,    &
2170              ' Number of ODEs   = ',I3,/,   ' Error tolerance  = ',D15.5,/, &
2171              ' Initial solution = ',/,(D15.5))
217290008 FORMAT (' Computed and reference solutions and absolute' &
2173              ' errors follow:',/,(I3,3D15.5))
217490009 FORMAT (' Steps = ',I10,' f-s = ',I10,' J-s = ',I10)
217590010 FORMAT (' ____________________________________________________________')
217690011 FORMAT (/,' Errors occurred. (Search on ISTATE.)')
217790012 FORMAT (/,' No errors occurred.')
217890013 FORMAT (' For ITOL = ', I3,' Dense  max. observed error = ',D15.5)
217990014 FORMAT (I5,6X,I5,3D15.5)
218090015 FORMAT (/,'Tolerance',3X,'Problem',4X,'Dense',10X,'Sparse',9X,'Banded')
218190016 FORMAT (//,20X,' Results for ITOL = ',I3)
218290017 FORMAT (' For ITOL = ', I3,' Sparse max. observed error = ',D15.5)
218390018 FORMAT (' For ITOL = ', I3,' Banded max. observed error = ',D15.5)
218490019 FORMAT (' The dense and banded results are not identical.')
218590020 FORMAT (/,' Maximum Errors for All Problems Follow:')
218690021 FORMAT (/,10X,' Maximum Errors for Individual Problems Follow:')
218790022 FORMAT (' The original VODE Jacobian was used.')
218890023 FORMAT (' The JACSP Jacobian was used.')
218990024 FORMAT (/,' Integration Statistics Follow:')
219090025 FORMAT (/,' Total number of solutions performed = ',I5)
219190026 FORMAT (I5,'   Steps: ',3I10)
219290027 FORMAT (5X,'   Jacs:  ',3I10)
219390028 FORMAT (5X,'   Ders:  ',3I10)
219490029 FORMAT (/,' Problem',12X,'Dense',4X,'Sparse',4X,'Banded')
219590030 FORMAT (' The ODEs were scaled.')
219690031 FORMAT (' The ODEs were not scaled.')
219790032 FORMAT (//,' To locate the summary results, search on the phrase',/ &
2198                 ' Summary Statistics Follow in the output file',         &
2199                 ' (stiffoptions.dat).',//)
220090033 FORMAT (/,' Tolerance', 10X, 'Total Execution Time')
220190034 FORMAT ((4X,I5,10X,I10))
220290035 FORMAT (/,' Problem No. = ',I3,/,' ITOL = ',I3,/,' Scaled ODEs', &
2203              13X,'Absolute Errors')
220490036 FORMAT (' Problem No. = ',I3,/,' ITOL = ',I3,/,' Unscaled ODEs', &
2205              18X,'Absolute Errors')
220690037 FORMAT (' Dense : JACSP - VODE: ',2D15.5)
220790038 FORMAT (' Sparse: JACSP - VODE: ',2D15.5)
220890039 FORMAT (' Banded: JACSP - VODE: ',2D15.5)
220990040 FORMAT (' Nonnegativity will be enforced for this problem.')
221090041 FORMAT (' Dense  SDJ: ',6I7)
221190042 FORMAT (' Sparse SDJ: ',6I7)
221290043 FORMAT (' Banded SDJ: ',6I7)
221390044 FORMAT (19X,'J',6X,'V',6X,'J',6X,'V',6X,'J',6X,'V')
221490045 FORMAT (///,20X,'Summary Statistics Follow')
2215
2216      STOP
2217    END PROGRAM DEMOSTIFF
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |