-
Notifications
You must be signed in to change notification settings - Fork 0
/
DVODE_F90_M.f90
18952 lines (17993 loc) · 716 KB
/
DVODE_F90_M.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
! If your usage of this software requires a license, the following
! BSD-like license is applicable. If it is not appropriate for your
! usage, please contact the authors.
!
! Copyright (c) 2013
! G.D. Byrne Illinois Institute of Technology, Chicago, ILL
! S. Thompson, Radford University, Radford, VA
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
! Redistributions of source code must retain the above copyright
! notice, this list of conditions, and the following disclaimer.
!
! Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions, and the following disclaimer
! in the documentation and/or other materials provided with the
! distribution.
!
! Neither the name of the organizations nor the names of its
! contributors may be used to endorse or promote products derived
! from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
MODULE DVODE_F90_M
! This version is the August 2013 release.
! Last change: 08/23/13
! _____________________________________________________________________
! Working Precision
IMPLICIT NONE
! Define the working precision for DVODE_F90. Change D0 to E0 in the
! next statement to convert to single precision.
INTEGER, PARAMETER, PRIVATE :: WP = KIND(1.0D0)
! ______________________________________________________________________
! Overview
! The f77 ordinary differential equation solver VODE.f is applicable to
! nonstiff systems of odes and to stiff systems having dense or banded
! Jacobians. DVODE_F90 is a Fortran 90 extension of VODE.f. While
! retaining all of the features available in VODE.f, we have
! incorporated several new options in DVODE_F90 including:
! 1. the ability to solve stiff systems with sparse Jacobians
! 2. internal management of storage and work arrays
! 3. specification of options via optional keywords
! 4. the ability to perform root finding or "event detection"
! 5. various new diagnostic and warning messages
! 6. the ability to impose solution bounds
! 7. several specialized options for dealing with sparsity
! ______________________________________________________________________
! Version Information
! This is DVODE_F90, the double precision FORTRAN 90 extension of the
! f77 DVODE.f ordinary differential equation solver. This version uses
! MA28 for sparse Jacobians. This file and related information can be
! obtained at the following support page:
!
! http:https://www.radford.edu/~thompson/vodef90web/
!
! We are indebted to Richard Cox (ORNL) for providing us with his
! implementation of MA28 in LSOD28.f (a variant of Alan Hindmarsh's
! lsodes.f). We are indebted to Alan Hindmarsh for numerous contributions.
! In particular, we borrowed liberally from the f77 solvers VODE.f,
! LSODAR.f, and LSODES.f while developing DVODE_F90. We are indebted
! to Doug Salane for providing us with his JACSP Jacobian routines.
!
! If you find a bug or encounter a problem with DVODE_F90, please
! contact one of us:
! G.D. Byrne ([email protected])
! S. Thompson ([email protected])
! A set of quick start instructions is provided below.
! ______________________________________________________________________
! Note on F90/F95 Compilers
! To date we have used DVODE_F90 successfully with all F90/F95 compilers
! to which we have access. In particular, we have used it with the Lahey
! F90 and Lahey-Fujitsu F95 compilers, the Compaq Visual F90 compiler,
! the g90 compiler, and the INTEL, Portland, Salford, and SUN compilers.
! It should be noted that compilers such as Salford's FTN95 complain
! about uninitialized arrays passed as subroutine arguments and the use of
! slices of two dimensional arrays as one dimensional vectors, and will
! not run using the strictest compiler options. It is perfectly safe to
! use the /-CHECK compiler option to avoid these FTN95 runtime checks.
! DVODE_F90 does not use any variable for numerical purposes until it
! has been assigned an appropriate value.
! ______________________________________________________________________
! Quick Start Instructions
! (1) Compile this file. Then compile, link, and execute the program
! example1.f90. The output is written to the file example1.dat.
! Verify that the last line of the output is the string
! 'No errors occurred.'
! (2) Repeat this process for the program example2.f90.
!
! Other test programs you may wish to run to verify your installation
! of DVODE_F90 are:
!
! (3) Run the test programs nonstiffoptions.f90 and stiffoptions.f90
! and verify that the last line in the output files produced is
! 'No errors occurred.' They solve the problems in the Toronto
! test suites using several different error tolerances and various
! solution options. Note that stiffoptions.f90 takes several
! minutes to run because it performs several thousand separate
! integrations.
! (4) Locate the file robertson.f90 in the demo programs and look at
! how options are set using SET_OPTS, how DVODE_F90 is called to
! obtain the solution at desired output times, and how the
! derivative and Jacobian routines are supplied. Note too the
! manner in which the solution is constrained to be nonnegative.
! (5) Locate demoharmonic.f90 and look at how root finding options
! are set and how the event residual routine is supplied to
! DVODE_F90.
! (6) The other demo programs available from the DVODE_F90 support
! page illustrate various other solution options available in
! DVODE_F90. The demo programs may be obtained from
!
! http:https://www.radford.edu/~thompson/vodef90web/index.html
! ______________________________________________________________________
! DVODE_F90 Full Documentation Prologue
! Section 1. Setting Options in DVODE_F90
! Section 2. Calling DVODE_F90
! Section 3. Choosing Error Tolerances
! Section 4. Choosing the Method of Integration
! Section 5. Interpolation of the Solution and Derivatives
! Section 6. Handling Events (Root Finding)
! Section 7. Gathering Integration Statistics
! Section 8. Determining Jacobian Sparsity Structure Arrays
! Section 9. Original DVODE Documentation Prologue
! Section 10. Example Usage
! Note: Search on the string 'Section' to locate these sections. You
! may wish to refer to the support page which has the sections broken
! into smaller pieces.
! ______________________________________________________________________
! Section 1. Setting Options in DVODE_F90
!
! You can use any of three options routines:
!
! SET_NORMAL_OPTS
! SET_INTERMEDIATE_OPTS
! SET_OPTS
! OPTIONS = SET_NORMAL_OPTS(DENSE_J, BANDED_J, SPARSE_J, &
! USER_SUPPLIED_JACOBIAN, LOWER_BANDWIDTH, UPPER_BANDWIDTH, &
! RELERR, ABSERR, ABSERR_VECTOR, NEVENTS)
! OPTIONS = SET_INTERMEDIATE_OPTS(DENSE_J, BANDED_J, SPARSE_J, &
! USER_SUPPLIED_JACOBIAN,LOWER_BANDWIDTH, UPPER_BANDWIDTH, &
! RELERR, ABSERR, ABSERR_VECTOR,TCRIT, H0, HMAX, HMIN, MAXORD, &
! MXSTEP, MXHNIL, NZSWAG, USER_SUPPLIED_SPARSITY, MA28_RPS, &
! NEVENTS, CONSTRAINED, CLOWER, CUPPER, CHANGE_ONLY_f77_OPTIONS) &
! OPTIONS = SET_OPTS(METHOD_FLAG, DENSE_J, BANDED_J, SPARSE_J, &
! USER_SUPPLIED_JACOBIAN, SAVE_JACOBIAN, CONSTANT_JACOBIAN, &
! LOWER_BANDWIDTH, UPPER_BANDWIDTH, SUB_DIAGONALS, SUP_DIAGONALS, &
! RELERR, RELERR_VECTOR, ABSERR, ABSERR_VECTOR, TCRIT, H0, HMAX, &
! HMIN, MAXORD, MXSTEP, MXHNIL, YMAGWARN, SETH, UPIVOT, NZSWAG, &
! USER_SUPPLIED_SPARSITY, NEVENTS, CONSTRAINED, CLOWER, CUPPER, &
! MA28_ELBOW_ROOM, MC19_SCALING, MA28_MESSAGES, MA28_EPS, &
! MA28_RPS, CHANGE_ONLY_f77_OPTIONS, JACOBIAN_BY_JACSP)
! Please refer to the documentation prologue for each of these functions
! to see what options may be used with each. Note that input to each is
! via keyword and all variables except error tolerances are optional.
! Defaults are used for unspecified options. If an option is available
! in SET_NORMAL OPTS, it is available and has the same meaning in
! SET_INTERMEDIATE_OPTS and SET_OPTS. Similarly, if an option is available
! in SET_INTERMEDIATE_OPTS, it is available and has the same meaning in
! SET_OPTS.
! The first two functions are provided merely for convenience.
! SET_NORMAL_OPTS is available simply to relieve you of reading the
! documentation for SET_OPTS and to use default values for all but
! the most common options. SET_INTERMEDIATE_OPTS is available to allow
! you more control of the integration while still using default values
! for less commonly used options. SET_OPTS allows you to specify any
! of the options available in DVODE_F90.
! Roughly, SET_NORMAL_OPTS is intended to provide for dense, banded,
! and numerical sparse Jacobians without the need to specify other
! specialized options. SET_INTERMEDIATE_OPTIONS is intended to allow
! more general sparse Jacobian options. SET_OPTS is intended to provide
! access to all options in DVODE_F90.
! Please note that SET_INTERMEDIATE_OPTS can be invoked using the same
! arguments as SET_NORMAL_OPTS; and SET_OPTS can be invoked using the
! same arguments as either SET_NORMAL_OPTS or SET_INTERMEDIATE_OPTS.
! If you wish you can simply delete SET_NORMAL_OPTS as well as
! SET_INTERMEDIATE_OPTS and use only SET_OPTS for all problems. If you
! do so, you need only include the options that you wish to use when
! you invoke SET_OPTIONS.
! In the following description any reference to SET_OPTS applies equally
! to SET_NORMAL_OPTS and SET_INTERMEDIATE OPTS.
! Before calling DVODE_F90 for the first time, SET_OPTS must be invoked.
! Typically, SET_OPTS is called once to set the desired integration
! options and parameters. DVODE_F90 is then called in an output loop to
! obtain the solution for the desired times. A detailed description of
! the DVODE_F90 arguments is given in a section below. Detailed descriptions
! of the options available via SET_OPTS are given in the documentation prologue.
! Although each option available in the f77 version of DVODE as well as
! several additional ones are available in DVODE_F90 via SET_OPTS,
! several of the available options are not relevant for most problems
! and need not be specified. Refer to the accompanying demonstration
! programs for specific examples of each usage. Note that after any call
! to DVODE_F90, you may call GET_STATS to gather relevant integration
! statistics. After your problem has completed, you may call
! RELEASE_ARRAYS to deallocate any internal arrays allocated by
! DVODE_F90 and to determine how much storage was used by DVODE_F90.
!
! To communicate with DVODE_F90 you will need to include the following
! statement in your calling program:
! USE DVODE_F90_M
! and include the following statement in your type declarations section:
! TYPE(VODE_OPTS) :: OPTIONS
! Below are brief summaries of typical uses of SET_OPTS.
! Nonstiff Problems:
! OPTIONS = SET_OPTS(RELERR=RTOL,ABSERR=ATOL)
! The above use of SET_OPTS will integrate your system of odes
! using the nonstiff Adams methods while using a relative error
! tolerance of RTOL and an absolute error tolerance of ATOL.
! Your subsequent call to DVODE_F90 might look like:
! CALL DVODE_F90(F,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS)
! OPTIONS = SET_OPTS(RELERR=RTOL,ABSERR=ATOL,NEVENTS=NG)
! If you wish to do root finding, SET_OPTS can be used as above.
! Here, NEVENTS is the desired number of root finding functions.
! Your subsequent call to DVODE_F90 might look like:
! CALL DVODE_F90(F,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,G_FCN=G)
! Here F is the name of your derivative subroutine and G is the
! name of your subroutine to evaluate residuals of the root
! finding functions.
! OPTIONS = SET_OPTS(RELERR=RTOL,ABSERR_VECTOR=ATOL)
! This use of SET_OPTS indicates that a scalar relative error
! tolerance and a vector of absolute error tolerances will be
! used.
! Stiff Problems, internally generated dense Jacobian:
! OPTIONS = SET_OPTS(DENSE_J=.TRUE.,RELERR=RTOL,ABSERR=ATOL)
! This use of DENSE_J=.TRUE. indicates that DVODE_F90 will
! use the stiffly stable BDF methods and will approximate
! the Jacobian, considered to be a dense matrix, using
! finite differences. Your subsequent call to DVODE_F90
! might look like:
! CALL DVODE_F90(F,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS)
! OPTIONS = SET_OPTS(DENSE_J=.TRUE.,ABSERR=ATOL,RELERR=RTOL, &
! USER_SUPPLIED_JACOBIAN=.TRUE.)
! If you know the Jacobian and wish to supply subroutine JAC
! as described in the documentation for DVODE_F90, the options
! call could look like the above.
! Your subsequent call to DVODE_F90 might look like:
! CALL DVODE_F90(F1,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JAC)
! Here, JAC is the name of the subroutine that you provide to
! evaluate the known Jacobian.
! Stiff Problems, internally generated banded Jacobian:
! OPTIONS = SET_OPTS(BANDED_J=.TRUE.,RELERR=RTOL,ABSERR=ATOL, &
! LOWER_BANDWIDTH=ML,UPPER_BANDWIDTH=MU)
! This use of BANDED_J=.TRUE. indicates that DVODE_F90 will
! use the stiffly stable BDF methods and will approximate the
! Jacobian, considered to be a banded matrix, using finite
! differences. Here ML is the lower bandwidth of the Jacobian
! and ML is the upper bandwidth of the Jacobian.
! Stiff Problems, internally generated sparse Jacobian:
! OPTIONS = SET_OPTS(SPARSE_J=.TRUE.,ABSERR=ATOL,RELERR=RTOL)
! This use of SET_OPTS indicates that the Jacobian is a sparse
! matrix. Its structure will be approximated internally by
! making calls to your derivative routine. If you know the
! structure before hand, you may provide it directly in a
! variety of ways as described in the documentation prologue
! for SET_OPTS. In addition, several other options related
! to sparsity are available.
! More complicated common usage:
! Suppose you need to solve a stiff problem with a sparse Jacobian.
! After some time, the structure of the Jacobian changes and you
! wish to have DVODE_F90 recalculate the structure before continuing
! the integration. Suppose that initially you want to use an absolute
! error tolerance of 1.0D-5 and that when the Jacobian structure is
! changed you wish to reduce the error tolerance 1.0D-7. Your calls
! might look like this.
! RTOL = ...
! ATOL = 1.0D-5
! OPTIONS = SET_OPTS(SPARSE_J=.TRUE.,ABSERR=ATOL,RELERR=RTOL)
! Output loop:
! CALL DVODE_F90(FCN,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS)
! At desired time:
! ISTATE = 3
! ATOL = 1.0D-7
! OPTIONS = SET_OPTS(SPARSE_J=.TRUE.,ABSERR=ATOL,RELERR=RTOL)
! End of output loop
! In the following we have summarized and described how some of the demonstration
! programs set options and call DVODE_F90. In each case the necessary parameters
! are defined before invoking SET_OPTS. The call to DVODE_F90 is in a loop in
! which the output time is successively updated. The actual programs are available
! from the web support page
!
! http:https://www.radford. edu/~thompson/vodef90web/index.html/
!
! Problem Summary
!
! Problem NEQ Jacobian Features Illustrated
!
! Prologue Example 1 3 Dense Basic
!
! Prologue Example 2 3 Dense Root finding
!
! Robertson 3 Dense Solution bounds
!
! Harmonic Oscillator 4 Nonstiff Root finding
!
! Flow Equations 5-1800 Sparse Automatic determination
! of sparsity arrays
!
! Diurnal Kinetics 50-5000 Sparse or banded Sparsity options
!
! Options Used and DVODE_F90 Call
!
! Prologue Example 1
!
! OPTIONS = SET_NORMAL_OPTS(DENSE_J=.TRUE., ABSERR_VECTOR=ATOL, RELERR=RTOL, &
! USER_SUPPLIED_JACOBIAN=.TRUE.)
! CALL DVODE_F90(FEX,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JEX)
!
! The problem consists of a stiff system of NEQ=3 equations. The dense
! Jacobian option (DENSE_J) is used. A vector ATOL(*) of error tolerances
! is used. A scalar relative error tolerance RTOL is used. Subroutine JEX
! is provided to evaluate the analytical Jacobian. If the last argument
! J_FCN=JEX is omitted (as in Example 2), a numerical Jacobian will
! be used.
!
! Prologue Example 2
!
! OPTIONS = SET_NORMAL_OPTS(DENSE_J=.TRUE., RELERR=RTOL, ABSERR_VECTOR=ATOL, &
! NEVENTS=NG)
! CALL DVODE_F90(FEX,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,G_FCN=GEX)
!
! The system in Example 1 is used to illustrate root finding. It is
! desired to locate the times at which two of the solution components
! attain prescribed values. NEVENTS=2 informs the solver that two such
! functions are used. Subroutine GEX is used to calculate the residuals
! for these two functions. A dense numerical Jacobian is used.
!
! Robertson Problem
!
! OPTIONS = SET_INTERMEDIATE_OPTS(DENSE_J=.TRUE., RELERR_VECTOR=RTOL, &
! ABSERR_VECTOR=ABSERR_TOLERANCES, CONSTRAINED=BOUNDED_COMPONENTS, &
! CLOWER=LOWER_BOUNDS, CUPPER=UPPER_BOUNDS)
!
! CALL DVODE_F90(DERIVS,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JACD)
! The system used in Examples 1 and 2 is solved over a much larger
! time interval. The solution is constrained to be nonegative. This
! is done by prescribing the components to be constrained (BOUNDED_COMPONENTS).
! Artificially large values are used to impose upper bounds (UPPER_BOUNDS)
! and lower bounds of zero are used to force a nonnegative solution.
!
! Harmonic Oscillator Problem
!
! OPTIONS = SET_NORMAL_OPTS(RELERR=RTOL, ABSERR=ATOL, NEVENTS=NG)
! CALL DVODE_F90(DERIVS,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,G_FCN=GEVENTS)
!
! A nonstiff system of NEQ=4 equations is solved. The nonstiff option is
! used because neither DENSE_ nor BANDED_J nor SPARSE_J is present. It is
! desired to find the times at which Y(2) or Y(3) is equal to 0. Residuals
! for the two corresponding event functions are calculated in subroutine
! GEVENTS.
!
! Flow Equations Problem
!
! OPTIONS = SET_OPTS(SPARSE_J=SPARSE, ABSERR=ATOL(1), RELERR=RTOL(1), &
! MXSTEP=100000, NZSWAG=20000)
! CALL DVODE_F90(FCN,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS)
!
! This is a stiff system of equations resulting a method of lines
! discretization. The Jacobian is sparse. Scalar absolute and relative
! error tolerances are used. The Jacobian structure and a numerical
! Jacobian are used. The solver is limited to a maximum of MXSTEP steps.
! NZSWAG is the amount by which allocated array sizes will be increased.
! The accompanying test program may be used to illutrate several other
! solution options.
!
! Diurnal Kinetics Problem
!
! OPTIONS = SET_OPTS(SPARSE_J=SPARSE, BANDED_J=BANDED, DENSE_J=DENSE, &
! ABSERR_VECTOR=ATOL(1:NEQ), RELERR=RTOL(1), MXSTEP=100000, &
! NZSWAG=50000, HMAX=MAXH, LOWER_BANDWIDTH=ML, UPPER_BANDWIDTH=MU, &
! MA28_ELBOW_ROOM=10, MC19_SCALING=.TRUE., MA28_MESSAGES=.FALSE., &
! MA28_EPS=1.0D-4, MA28_RPS=.TRUE., &
! USER_SUPPLIED_SPARSITY=SUPPLY_STRUCTURE)
! CALL USERSETS_IAJA(IA, IADIM, JA, JADIM)
! CALL DVODE_F90(FCN, NEQ, Y, T, TOUT, ITASK, ISTATE, OPTIONS)
!
! This problem can be used to illustrate most solution options. Here, dense,
! banded, or sparse Jacobians are used depending on the values of the first
! three parameters. A vector error tolerance is used and a scalar relative
! error tolerance is used. If a banded solution is desired, it is necessary
! to supply the bandwidths ML and MU. If a sparse solution is desired,
! several special options are used. The most important one is MA28_RPS to
! force the solver to update the partial pivoting sequence when singular
! iteration matrices are encountered. The sparsity pattern is determined
! numerically if SUPPLY_STRUCTURE is FALSE. Otherwise the user will supply
! the pattern by calling subroutine USERSETS_IAJA.
! ______________________________________________________________________
! Section 2. Calling DVODE_F90
!
! DVODE_F90 solves the initial value problem for stiff or nonstiff
! systems of first order ODEs,
! dy/dt = f(t,y), or, in component form,
! dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ).
! DVODE_F90 is a package based on the EPISODE and EPISODEB packages,
! and on the ODEPACK user interface standard. It was developed from
! the f77 solver DVODE developed by Brown, Byrne, and Hindmarsh.
! DVODE_F90 also provides for the solution of sparse systems in a
! fashion similar to LSODES and LSOD28. Currently, MA28 is used
! to perform the necessary sparse linear algebra. DVODE_F90 also
! contains the provision to do root finding in a fashion similar
! to the LSODAR solver from ODEPACK.
! Communication between the user and the DVODE_F90 package, for normal
! situations, is summarized here. This summary describes only a subset
! of the full set of options available. See the full description for
! details, including optional communication, nonstandard options, and
! instructions for special situations.
! CALL DVODE_F90(F,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JAC,G_FCN=GEX)
! The arguments in the call list to DVODE_F90 have the following
! meanings.
! F = The name of the user-supplied subroutine defining the
! ODE system. The system must be put in the first-order
! form dy/dt = f(t,y), where f is a vector-valued function
! of the scalar t and the vector y. Subroutine F is to
! compute the function f. It is to have the form
! SUBROUTINE F(NEQ,T,Y,YDOT)
! DOUBLE PRECISION T,Y(NEQ),YDOT(NEQ)
! where NEQ, T, and Y are input, and the array YDOT = f(t,y)
! is output. Y and YDOT are arrays of length NEQ.
! Subroutine F should not alter Y(1),...,Y(NEQ).
! If F (and JAC) are not contained in a module available to
! your calling program, you must declare F to be EXTERNAL
! in the calling program.
! NEQ = The size of the ODE system (number of first order
! ordinary differential equations).
! Y = A double precision array for the vector of dependent variables,
! of length NEQ or more. Used for both input and output on the
! first call (ISTATE = 1), and only for output on other calls.
! On the first call, Y must contain the vector of initial
! values. In the output, Y contains the computed solution
! evaluated at T.
! T = The independent variable. In the input, T is used only on
! the first call, as the initial point of the integration.
! In the output, after each call, T is the value at which a
! computed solution Y is evaluated (usually the same as TOUT).
! On an error return, T is the farthest point reached.
! TOUT = The next value of t at which a computed solution is desired.
! TOUT is Used only for input. When starting the problem
! (ISTATE = 1), TOUT may be equal to T for one call, then
! should not equal T for the next call. For the initial T,
! an input value of TOUT unequal to T is used in order to
! determine the direction of the integration (i.e. the
! algebraic sign of the step sizes) and the rough scale
! of the problem. Integration in either direction (forward
! or backward in t) is permitted. If ITASK = 2 or 5 (one-step
! modes), TOUT is ignored after the first call (i.e. the
! first call with TOUT \= T). Otherwise, TOUT is required
! on every call. If ITASK = 1, 3, or 4, the values of TOUT
! need not be monotone, but a value of TOUT which backs up
! is limited to the current internal t interval, whose
! endpoints are TCUR - HU and TCUR. (Refer to the description
! of GET_STATS for a description of TCUR and HU.)
! ITASK = An index specifying the task to be performed.
! Input only. ITASK has the following values and meanings.
! 1 means normal computation of output values of y(t) at
! t = TOUT (by overshooting and interpolating).
! 2 means take one step only and return.
! 3 means stop at the first internal mesh point at or
! beyond t = TOUT and return.
! 4 means normal computation of output values of y(t) at
! t = TOUT but without overshooting t = TCRIT.
! TCRIT must be specified in your SET_OPTS call. TCRIT
! may be equal to or beyond TOUT, but not behind it in
! the direction of integration. This option is useful
! if the problem has a singularity at or beyond t = TCRIT.
! 5 means take one step, without passing TCRIT, and return.
! TCRIT must be specified in your SET_OPTS call.
! If ITASK = 4 or 5 and the solver reaches TCRIT (within
! roundoff), it will return T = TCRIT(exactly) to indicate
! this (unless ITASK = 4 and TOUT comes before TCRIT, in
! which case answers at T = TOUT are returned first).
! ISTATE = an index used for input and output to specify the
! the state of the calculation.
! In the input, the values of ISTATE are as follows.
! 1 means this is the first call for the problem
! (initializations will be done). See note below.
! 2 means this is not the first call, and the calculation
! is to continue normally, with no change in any input
! parameters except possibly TOUT and ITASK.
! 3 means this is not the first call, and the
! calculation is to continue normally, but with
! a change in input parameters other than
! TOUT and ITASK. Desired changes require SET_OPTS
! be called prior to calling DVODE_F90 again.
! A preliminary call with TOUT = T is not counted as a
! first call here, as no initialization or checking of
! input is done. (Such a call is sometimes useful to
! include the initial conditions in the output.)
! Thus the first call for which TOUT is unequal to T
! requires ISTATE = 1 in the input.
! In the output, ISTATE has the following values and meanings.
! 1 means nothing was done, as TOUT was equal to T with
! ISTATE = 1 in the input.
! 2 means the integration was performed successfully.
! 3 means a root of one of your root finding functions
! has been located.
! A negative value of ISTATE indicates that DVODE_F90
! encountered an error as described in the printed error
! message. Since the normal output value of ISTATE is 2,
! it does not need to be reset for normal continuation.
! Also, since a negative input value of ISTATE will be
! regarded as illegal, a negative output value requires
! the user to change it, and possibly other input, before
! calling the solver again.
! OPTIONS = The options structure produced by your call to SET_OPTS.
! JAC = The name of the user-supplied routine (MITER = 1 or 4 or 6)
! If you do not specify that a stiff method is to be used
! in your call to SET_OPTS, you need not include JAC in
! your call to DVODE_F90. If you specify a stiff method and
! that a user supplied Jacobian will be supplied, JAC must
! compute the Jacobian matrix, df/dy, as a function of the
! scalar t and the vector y. It is to have the form:
! SUBROUTINE JAC(NEQ, T, Y, ML, MU, PD, NROWPD)
! DOUBLE PRECISION T, Y(NEQ), PD(NROWPD,NEQ)
! where NEQ, T, Y, ML, MU, and NROWPD are input and the array
! PD is to be loaded with partial derivatives (elements of the
! Jacobian matrix) in the output. PD must be given a first
! dimension of NROWPD. T and Y have the same meaning as in
! Subroutine F.
! In the full matrix case (MITER = 1), ML and MU are
! ignored, and the Jacobian is to be loaded into PD in
! columnwise manner, with df(i)/dy(j) loaded into PD(i,j).
! In the band matrix case (MITER = 4), the elements
! within the band are to be loaded into PD in columnwise
! manner, with diagonal lines of df/dy loaded into the rows
! of PD. Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j).
! ML and MU are the half-bandwidth parameters. (See IUSER).
! The locations in PD in the two triangular areas which
! correspond to nonexistent matrix elements can be ignored
! or loaded arbitrarily, as they are overwritten by DVODE_F90.
! In the sparse matrix case the elements of the matrix
! are determined by the sparsity structure given by the
! IA and JA pointer arrays. Refer to the documentation
! prologue for SET_OPTS for a description of the arguments
! for JAC since they differ from the dense and banded cases.
! JAC need not provide df/dy exactly. A crude
! approximation (possibly with a smaller bandwidth) will do.
! In either case, PD is preset to zero by the solver,
! so that only the nonzero elements need be loaded by JAC.
! In the sparse matrix case, JAC has a different form:
! SUBROUTINE JAC (N, T, Y, IA, JA, NZ, PD)
! Given the number of odes N, the current time T, and the
! current solution vector Y, JAC must do the following:
! If NZ = 0 on input:
! Replace NZ by the number of nonzero elements in the
! Jacobian. The diagonal of the Jacobian must be included.
! Do NOT define the arrays IA, JA, PD at this time.
! Once JAC has been called with NZ = 0 and you have
! defined the value of NZ, future calls to JAC will use
! this value of NZ.
! When a call is made with NZ unequal to 0, you must
! define the sparsity structure arrays IA and JA, and
! the sparse Jacobian PD.
! IA defines the number of nonzeros including the
! diagonal in each column of the Jacobian. Define
! IA(1) = 1 and for J = 1,..., N,
! IA(J+1) = IA(J) + number of nonzeros in column J.
! Diagonal elements must be included even if they are
! zero. You should check to ensure that IA(N+1)-1 = NZ.
! JA defines the rows in which the nonzeros occur.
! For I = 1,...,NZ, JA(I) is the row in which the Ith
! element of the Jacobian occurs. JA must also include
! the diagonal elements of the Jacobian.
! PD defines the numerical value of the Jacobian
! elements. For I = 1,...,NZ, PD(I) is the numerical
! value of the Ith element in the Jacobian. PD must
! also include the diagonal elements of the Jacobian.
! GFUN = the name of the subroutine to evaluate the residuals for
! event functions. If you do not specify that events are
! present (by specifying NEVENTS > 0 in SET_OPTS), you
! need not include GFUN in your call list for DVODE_F90.
! If GFUN is not contained in a module available to your
! calling program, you must declare GFUN to be EXTERNAL
! in your calling program.
! To continue the integration after a successful return, simply
! reset TOUT and call DVODE_F90 again. No other parameters need
! be reset unless ISTATE=3 in which case, reset it to 2 before
! calling DVODE_F90 again.
! ______________________________________________________________________
! Section 3. Choosing Error Tolerances
!
! This is the most important aspect of solving odes numerically.
! You may supply any of four keywords and values. If you wish to
! use scalar error tolerances, you should supply ABSERR and RELERR.
! For a good many problems, it is advisable to supply a vector of
! absolute error tolerances ABSERR_VECTOR = desired vector. This
! allows you to use different tolerances for each component of
! the solution. If ABSERR_VECTOR is supplied, it must be a vector
! of length NEQ where NEQ is the number of odes in your system.
! Similarly, you may supply a vector of relative error tolerances,
! RELERR_VECTOR. If no tolerances are specified, DVODE_F90 will use
! default error tolerances ABSERR=1D-6 and RELERR=1D-4; but it is
! strongly recommended that you supply values that are appropriate
! for your problem. In the event you do not supply error tolerances,
! DVODE_F90 will print a reminder that the default error tolerances
! are not appropriate for all problems.
!
! RELERR can be set to a scalar value as follows.
! Determine the number of significant digits of accuracy desired,
! which will be a positive integer, say, N.
! Then set RELERR = 10**-(N+1).
! The authors recommend that RELERR be no larger than 10**-4.
! The authors recommend a vector valued absolute error tolerance,
! which can be set as follows.
! For the I-th component of the solution vector, Y(I), determine
! the positive number FLOOR(I) at which ABS(Y(I)) becomes
! negligible for the problem at hand. FLOOR(I) is sometimes called
! the problem zero or the floor value for the I-th component and is
! problem dependent. For a given problem that is not scaled, these
! floor values may well vary by up to 9 orders of magnitude.
! Set ABSERR(I) = FLOOR(I) or to be conservative
! ABSERR_VECTOR(I) = 0.1*FLOOR(I). There is no variable FLOOR in
! DVODE_F90. If it is difficult to divine the components of ABSERR,
! (or FLOOR) make a reasonable guess, run the problem, then set that
! ABSERR_VECTOR so for I = 1, 2,...NEQ,
! ABSERR_VECTOR(I) = 1D-6*RELERR*MAX{ABS(Y(I,TOUT): for all TOUT}.
! The correct choices for RELERR and ABSERR can and do have
! significant impact on both the quality of the solution and run
! time. Counter intuitively, error tolerances that are too loose
! can and do increase run time significantly and the quality of
! the solution may be seriously compromised.
! Examples:
! 1. OPTIONS = SET_OPTS(DENSE_J=.TRUE., ABSERR=1D-8,RELERR=1D-8)
! This will yield MF = 22. Both the relative error tolerance
! and the absolute error tolerance will equal 1D-8.
! 2. OPTIONS = SET_OPTS(DENSE_J=.TRUE.,RELERR=1D-5, &
! ABSERR_VECTOR=(/1D-6,1D-8/))
! For a system with NEQ=2 odes, this will yield MF = 22. A scalar
! relative error tolerance equal to 1D-5 will be used. Component 1
! of the solution will use an absolute error tolerance of 1D-6
! while component 2 will use an absolute error tolerance of 1D-8.
! ______________________________________________________________________
! Section 4. Choosing the Method of Integration
!
! If you wish to specify a numerical value for METHOD_FLAG, it can equal
! any of the legal values of MF for DVODE or LSODES. (See below.) If
! you do not wish to specify a numerical value for METHOD_FLAG, you
! may specify any combination of the five logical keywords DENSE_J,
! BANDED_J, SPARSE_J, USER_SUPPLIED_JACOBIAN, SAVE_JACOBIAN that you
! wish. Appropriate values will be used in DVODE_F90 for any variables
! that are not present. The first three flags indicate the type of
! Jacobian, dense, banded, or sparse. If USER_SUPPLIED_JACOBIAN=.TRUE.,
! the Jacobian will be determined using the subroutine JAC you supply
! in your call to DVODE_F90. Otherwise, an internal Jacobian will be
! generated using finite differences.
! Examples:
! 1. OPTIONS = SET_OPTS(METHOD_FLAG=22,...)
! DVODE will use the value MF=22 as described in the documentation
! prologue. In this case, the stiff BDF methods will be used with
! a dense, internally generated Jacobian.
! 2. OPTIONS = SET_OPTS(METHOD_FLAG=21,...)
! This is the same an Example 1 except that a user supplied dense
! Jacobian will be used and DVODE will use MF=21.
! 3. OPTIONS = SET_OPTS(DENSE_J=.TRUE.,...)
! This will yield MF = 22 as in Example 1, provided
! USER_SUPPLIED_JACOBIAN and SAVE_JACOBIAN are not present, or if
! present are set to their default
! values of .FALSE. and .TRUE., respectively.
! 4. OPTIONS = SET_OPTS(DENSE_J=.TRUE.,&
! USER_SUPPLIED_JACOBIAN=.TRUE.,...)
! This will yield MF = 21 as in Example 2, provided SAVE_JACOBIAN
! is not present, or if present is set to its default value .FALSE.
! Notes:
! 1. If you specify more than one of DENSE_J, BANDED_J, and SPARSE_J,
! DENSE_J takes precedence over BANDED_J which in turn takes
! precedence over SPARSE_J.
! 2. By default, DVODE_F90 saves a copy of the Jacobian and reuses the
! copy when necessary. For problems in which storage requirements
! are acute, you may wish to override this default and have
! DVODE_F90 recalculate the Jacobian rather than use a saved copy.
! You can do this by specifying SAVE_JACOBIAN=.FALSE. It is
! recommended that you not do this unless necessary since it can
! have a significant impact on the efficiency of DVODE_F90. (For
! example, when solving a linear problem only one evaluation of
! the Jacobian is required with the default option.)
! 3. If you choose BANDED_J = .TRUE. or if you supply a value of MF
! that corresponds to a banded Jacobian, you must also supply the
! lower bandwidth ML and the upper bandwidth of the Jacobian MU
! by including the keywords
! LOWER_BANDWIDTH = value of ML and UPPER_BANDWIDTH = value of M
! More on Method Selection
! The keyword options available in SET_OPTS are intended to replace
! the original method indicator flag MF. However, if you wish to
! retain the flexibility of the original solver, you may specify MF
! directly in your call to SET_OPTS. This is done by using the
! keyword METHOD_FLAG=MF in your SET_OPTS where MF has any of the
! values in the following description. Refer to the demonstration
! program demosp.f90 for an example in which this is done.
! MF = The method flag. Used only for input. The legal values of
! MF are:
! 10, 11, 12, 13, 14, 15, 16, 17, 20, 21, 22, 23, 24, 25, 26,
! 27, -11, -12, -14, -15, -21, -22, -24, -25, -26, -27.
! MF is a signed two-digit integer, MF = JSV*(10*METH+MITER).
! JSV = SIGN(MF) indicates the Jacobian-saving strategy:
! JSV = 1 means a copy of the Jacobian is saved for reuse
! in the corrector iteration algorithm.
! JSV = -1 means a copy of the Jacobian is not saved
! (valid only for MITER = 1, 2, 4, or 5).
! METH indicates the basic linear multistep method:
! METH = 1 means the implicit Adams method.
! METH = 2 means the method based on backward
! differentiation formulas (BDF-s).
! MITER indicates the corrector iteration method:
! MITER = 0 means functional iteration (no Jacobian matrix
! is involved).
! MITER = 1 means chord iteration with a user-supplied
! full (NEQ by NEQ) Jacobian.
! MITER = 2 means chord iteration with an internally
! generated (difference quotient) full Jacobian
! (using NEQ extra calls to F per df/dy value).
! MITER = 4 means chord iteration with a user-supplied
! banded Jacobian.
! MITER = 5 means chord iteration with an internally
! generated banded Jacobian (using ML+MU+1 extra
! calls to F per df/dy evaluation).
! MITER = 6 means chord iteration with a user-supplied
! sparse Jacobian.
! MITER = 7 means chord iteration with an internally
! generated sparse Jacobian
! If MITER = 1, 4, or 6 the user must supply a subroutine
! JAC(the name is arbitrary) as described above under JAC.
! For other values of MITER, JAC need not be provided.
! ______________________________________________________________________
! Section 5. Interpolation of the Solution and Derivative
!
! Following a successful return from DVODE_F90, you may call
! subroutine DVINDY to interpolate the solution or derivative.
! SUBROUTINE DVINDY(T, K, DKY, IFLAG)
! DVINDY computes interpolated values of the K-th derivative of the
! dependent variable vector y, and stores it in DKY. This routine
! is called with K = 0 or K = 1 and T = TOUT. In either case, the
! results are returned in the array DKY of length at least NEQ which
! must be declared and dimensioned in your calling program. The
! computed values in DKY are obtained by interpolation using the
! Nordsieck history array.
! ______________________________________________________________________
! Section 6. Handling Events (Root Finding)
!
! DVODE_F90 contains root finding provisions. It attempts to
! locates the roots of a set of functions
! g(i) = g(i,t,y(1),...,y(NEQ)) (i = 1,...,ng).
! To use root finding include NEVENTS=NG in your call to SET_OPTS
! where NG is the number of root finding functions. You must then
! supply subroutine GFUN in your call to DVODE_F90 using
! G_FCN=GFUN as the last argument. GFUN must has the form
! SUBROUTINE GFUN(NEQ, T, Y, NG, GOUT)
! where NEQ, T, Y, and NG are input, and the array GOUT is output.
! NEQ, T, and Y have the same meaning as in the F routine, and
! GOUT is an array of length NG. For i = 1,...,NG, this routine is
! to load into GOUT(i) the value at (T,Y) of the i-th constraint
! function g(i). DVODE_F90 will find roots of the g(i) of odd
! multiplicity (i.e. sign changes) as they occur during the
! integration. GFUN must be declared EXTERNAL in the calling
! program. Note that because of numerical errors in the functions
! g(i) due to roundoff and integration error, DVODE_F90 may return
! false roots, or return the same root at two or more nearly equal
! values of t. This is particularly true for problems in which the
! integration is restarted (ISTATE = 1) at a root. If such false
! roots are suspected, you should consider smaller error tolerances
! and/or higher precision in the evaluation of the g(i). Note
! further that if a root of some g(i) defines the end of the
! problem, the input to DVODE_F90 should nevertheless allow
! integration to a point slightly past that root, so that DVODE_F90
! can locate the root by interpolation. Each time DVODE_F90 locates
! a root of one of your event functions it makes a return to the
! calling program with ISTATE = 3. When such a return is made and
! you have processed the results, simply change ISTATE = 2 and call
! DVODE_F90 again without making other changes.
! ______________________________________________________________________
! Section 7. Gathering Integration Statistics
!
! SUBROUTINE GET_STATS(RSTATS, ISTATS, NUMEVENTS, JROOTS)
! Caution:
! RSTATS and ISTATS must be declared and dimensioned in your
! main program. The minimum dimensions are:
! DIMENSION RSTATS(22), ISTATS(31)
! This subroutine returns the user portions of the original DVODE
! RUSER and IUSER arrays, and if root finding is being done, it
! returns the original LSODAR JROOT vector. NUMEVENTS and JROOTS
! are optional parameters. NUMEVENTS is the number of root functions
! and JROOTS is an integer array of length NUMEVENTS.
! Available Integration Statistics:
! HU RUSER(11) The step size in t last used (successfully).
! HCUR RUSER(12) The step size to be attempted on the next step.
! TCUR RUSER(13) The current value of the independent variable
! which the solver has actually reached, i.e. the
! current internal mesh point in t. In the output,
! TCUR will always be at least as far from the
! initial value of t as the current argument T,
! but may be farther (if interpolation was done).
! TOLSF RUSER(14) A tolerance scale factor, greater than 1.0,
! computed when a request for too much accuracy was
! detected (ISTATE = -3 if detected at the start of
! the problem, ISTATE = -2 otherwise). If ITOL is
! left unaltered but RTOL and ATOL are uniformly
! scaled up by a factor of TOLSF for the next call,
! then the solver is deemed likely to succeed.
! (The user may also ignore TOLSF and alter the
! tolerance parameters in any other way appropriate.)
! NST IUSER(11) The number of steps taken for the problem so far.
! NFE IUSER(12) The number of f evaluations for the problem so far.
! NJE IUSER(13) The number of Jacobian evaluations so far.
! NQU IUSER(14) The method order last used (successfully).
! NQCUR IUSER(15) The order to be attempted on the next step.
! IMXER IUSER(16) The index of the component of largest magnitude in
! the weighted local error vector (E(i)/EWT(i)),
! on an error return with ISTATE = -4 or -5.
! LENRW IUSER(17) The length of RUSER actually required.
! This is defined on normal returns and on an illegal
! input return for insufficient storage.
! LENIW IUSER(18) The length of IUSER actually required.
! This is defined on normal returns and on an illegal
! input return for insufficient storage.
! NLU IUSER(19) The number of matrix LU decompositions so far.
! NNI IUSER(20) The number of nonlinear (Newton) iterations so far.
! NCFN IUSER(21) The number of convergence failures of the nonlinear
! solver so far.
! NETF IUSER(22) The number of error test failures of the integrator
! so far.
! MA28AD_CALLS IUSER(23) The number of calls made to MA28AD
! MA28BD_CALLS IUSER(24) The number of calls made to MA28BD
! MA28CD_CALLS IUSER(25) The number of calls made to MA28CD
! MC19AD_CALLS IUSER(26) The number of calls made to MC19AD
! IRNCP IUSER(27) The number of compressions done on array JAN
! ICNCP IUSER(28) The number of compressions done on array ICN
! MINIRN IUSER(29) Minimum size for JAN array
! MINICN IUSER(30) Minimum size for ICN array
! MINNZ IUSER(31) Number of nonzeros in sparse Jacobian
! JROOTS JROOTS Optional array of component indices for components
! having a zero at the current time
! ______________________________________________________________________
! Section 8. Determining Jacobian Sparsity Structure Arrays
!
! If you are solving a problem with a sparse Jacobian, the arrays
! that define the sparsity structure are needed. The arrays may
! be determined in any of several ways.
! 1. If you choose the default mode by indicating SPARSE=.TRUE.,
! the sparsity arrays will be determined internally by DVODE_F90
! by making calls to your derivative subroutine. This mode is
! equivalent to using the integration method flag MF = 227.
! 2. The DVODE_F90 method flag MF is defined to be
! MF = 100*MOSS + 10*METH + MITER. If you supply MF = 227 (or 217),
! the sparse Jacobian will be determined using finite differences;
! and the sparsity arrays will be determined by calling your
! derivative subroutine.
! 3. If you supply MF = 126 (or 116), you must supply the Jacobian
! subroutine JAC to define the exact Jacobian. JAC must have the
! following form:
! SUBROUTINE JAC (N, T, Y, IA, JA, NZ, PD)
! Given the number of odes N, the current time T, and the current
! solution vector Y, JAC must do the following:
! - If NZ = 0 on input:
! Replace NZ by the number of nonzero elements in the Jacobian.
! The diagonal of the Jacobian must be included.
! Do NOT define the arrays IA, JA, PD at this time.
! Once JAC has been called with NZ = 0 and you have defined the
! value of NZ, future calls to JAC will use this value of NZ.
! - When a call is made with NZ unequal to 0, you must define the
! sparsity structure arrays IA and JA, and the sparse Jacobian
! PD.
! - IA defines the number of nonzeros including the diagonal
! in each column of the Jacobian. Define IA(1) = 1 and for
! J = 1,..., N,
! IA(J+1) = IA(J) + number of nonzeros in column J.
! Diagonal elements must be include even if they are zero.
! You should check to ensure that IA(N+1)-1 = NZ.
! - JA defines the rows in which the nonzeros occur. For
! I = 1,...,NZ, JA(I) is the row in which the Ith element
! of the Jacobian occurs. JA must also include the diagonal
! elements of the Jacobian.
! - PD defines the numerical value of the Jacobian elements.
! For I = 1,...,NZ, PD(I) is the numerical value of the
! Ith element in the Jacobian. PD must also include the
! diagonal elements of the Jacobian.
! 4. If you wish to supply the IA and JA arrays directly, use
! MF = 27. In this case, after calling SET_OPTS, you must call
! SET_IAJA supplying the arrays IAUSER and JAUSER described in
! the documentation prologue for SET_IAJA. These arrays will be
! used when approximate Jacobians are determined using finite
! differences.
! There are two user callable sparsity structure subroutines:
! USERSETS_IAJA may be used if you wish to supply the sparsity
! structure directly.
! SUBROUTINE USERSETS_IAJA(IAUSER,NIAUSER,JAUSER,NJAUSER)
! Caution:
! If it is called, USERSETS_IAJA must be called after the
! call to SET_OPTS.
! Usage:
! CALL SET_IAJA(IAUSER,NIAUSER,JAUSER,NJAUSER)
! In this case, IAUSER of length NIAUSER will be used for
! IA; and JAUSER of length NJAUSER will be used for JA.
! Arguments:
! IAUSER = user supplied IA array
! NIAUSER = length of IAUSER array
! JAUSER = user supplied JA vector
! NJAUSER = length of JAUSER array
! The second subroutine allows you to approximate the sparsity
! structure using derivative differences. It allows more flexibility
! in the determination of perturbation increments used.
! SUBROUTINE SET_IAJA(DFN,NEQ,T,Y,FMIN,NTURB,DTURB,IAUSER, &
! NIAUSER, JAUSER, NJAUSER)
! Caution:
! If it is called, SET_IAJA must be called after the call to
! SET_OPTS.
! Usage:
! SET_IAJA may be called in one of two ways:
! CALL SET_IAJA(DFN,NEQ,T,Y,FMIN,NTURB,DTURB)
! In this case IA and JA will be determined using calls
! to your derivative routine DFN.
! CALL SET_IAJA(DFN,NEQ,T,Y,FMIN,NTURB,DTURB,IAUSER,NIAUSER, &
! JAUSER, NJAUSER)
! In this case, IAUSER of length NIAUSER will be used for
! IA; and JAUSER of length NJAUSER will be used for JA.
! T, Y, FMIN, NTURB, and DTURB will be ignored (though
! they must be present in the argument list).
! Arguments:
! DFN = DVODE derivative subroutine
! NEQ = Number of odes
! T = independent variable t
! Y = solution y(t)
! FMIN = Jacobian threshold value. Elements of the Jacobian
! with magnitude smaller than FMIN will be ignored.
! FMIN will be ignored if it is less than or equal
! to ZERO.
! NTURB = Perturbation flag. If NTURB=1, component I of Y
! will be perturbed by 1.01D0.
! If NTURB=NEQ, component I of Y will be perturbed
! by ONE + DTURB(I).
! DTURB = perturbation vector of length 1 or NEQ.
! If these four optional parameters are present, IAUSER and JAUSER
! will be copied to IA and JA rather than making derivative calls
! to approximate IA and JA:
! IAUSER = user supplied IA array
! NIAUSER = length of IAUSER array
! JAUSER = user supplied JA vector
! NJAUSER = length of JAUSER array
! ______________________________________________________________________
! Section 9. Original DVODE.F Documentation Prologue
!
! SUBROUTINE DVODE(F, NEQ, Y, T, TOUT, ITASK, ISTATE, OPTS, JAC, GFUN)
! DVODE: Variable-coefficient Ordinary Differential Equation solver,
! with fixed-leading-coefficient implementation.
! Note:
! Numerous changes have been made in the documentation and the code
! from the original Fortran 77 DVODE solver. With regard to the new
! F90 version, if you choose options that correspond to options
! available in the original f77 solver, you should obtain the same
! results. In all testing, identical results have been obtained
! between this version and a simple F90 translation of the original
! solver.
! DVODE solves the initial value problem for stiff or nonstiff
! systems of first order ODEs,
! dy/dt = f(t,y), or, in component form,
! dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ).
! DVODE is a package based on the EPISODE and EPISODEB packages, and
! on the ODEPACK user interface standard, with minor modifications.
! This version is based also on elements of LSODES and LSODAR.
! Authors:
! Peter N. Brown and Alan C. Hindmarsh
! Center for Applied Scientific Computing, L-561
! Lawrence Livermore National Laboratory
! Livermore, CA 94551
! George D. Byrne
! Illinois Institute of Technology
! Chicago, IL 60616
! References:
! 1. P. N. Brown, G. D. Byrne, and A. C. Hindmarsh, "VODE: A Variable
! Coefficient ODE Solver," SIAM J. Sci. Stat. Comput., 10 (1989),
! pp. 1038-1051. Also, LLNL Report UCRL-98412, June 1988.
! 2. G. D. Byrne and A. C. Hindmarsh, "A Polyalgorithm for the