[2696] | 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 |
---|
| 162 | 10 IOUT = 6 |
---|
| 163 | WRITE (IOUT,FMT=90000) ID |
---|
| 164 | STOP |
---|
| 165 | |
---|
| 166 | ! PROBLEM CLASS A - LINEAR WITH REAL EIGENVALUES |
---|
| 167 | |
---|
| 168 | 20 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 | |
---|
| 183 | 30 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 | |
---|
| 203 | 40 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 | |
---|
| 218 | 50 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 | |
---|
| 241 | 60 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 | |
---|
| 257 | 70 CONTINUE |
---|
| 258 | ! PROBLEM B2, B3, B4, B5 |
---|
| 259 | N = 6 |
---|
| 260 | ITMP = IID - 1 |
---|
| 261 | GOTO (80,90,100,110) ITMP |
---|
| 262 | 80 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 |
---|
| 270 | 90 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 |
---|
| 278 | 100 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 |
---|
| 286 | 110 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 |
---|
| 293 | 120 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 | |
---|
| 305 | 130 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 | |
---|
| 320 | 140 CONTINUE |
---|
| 321 | ! PROBLEM C2, C3, C4, C5 |
---|
| 322 | N = 4 |
---|
| 323 | ITMP = IID - 1 |
---|
| 324 | GOTO (150,160,170,180) ITMP |
---|
| 325 | 150 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 |
---|
| 331 | 160 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 |
---|
| 337 | 170 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 |
---|
| 343 | 180 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 |
---|
| 348 | 190 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 | |
---|
| 359 | 200 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 | |
---|
| 373 | 210 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 | |
---|
| 387 | 220 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 | |
---|
| 403 | 230 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 | |
---|
| 417 | 240 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 | |
---|
| 429 | 250 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 | |
---|
| 445 | 260 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 | |
---|
| 460 | 270 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 | |
---|
| 472 | 280 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 | |
---|
| 486 | 290 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 | |
---|
| 502 | 300 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 | |
---|
| 520 | 310 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 | |
---|
| 536 | 320 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 | |
---|
| 548 | 330 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 | |
---|
| 567 | 340 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 | |
---|
| 582 | 350 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 |
---|
| 596 | 360 CONTINUE |
---|
| 597 | IF (IWT<0) GOTO 370 |
---|
| 598 | DO I = 1, N |
---|
| 599 | Y(I) = Y(I)/W(I) |
---|
| 600 | END DO |
---|
| 601 | 370 CONTINUE |
---|
| 602 | RETURN |
---|
| 603 | |
---|
| 604 | 90000 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 |
---|
| 655 | 10 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 |
---|
| 662 | 20 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 |
---|
| 674 | 30 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 |
---|
| 681 | 40 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 |
---|
| 696 | 50 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 |
---|
| 703 | 60 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 |
---|
| 712 | 70 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 |
---|
| 721 | 80 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 |
---|
| 730 | 90 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 |
---|
| 741 | 100 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 |
---|
| 748 | 110 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 |
---|
| 755 | 120 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 |
---|
| 762 | 130 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 |
---|
| 769 | 140 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 |
---|
| 778 | 150 Y(1) = 2.224222010616901D+01 |
---|
| 779 | Y(2) = 2.711071334484136D+01 |
---|
| 780 | Y(3) = 3.999999999999999D+02 |
---|
| 781 | GOTO 310 |
---|
| 782 | |
---|
| 783 | ! PROBLEM D2 |
---|
| 784 | 160 Y(1) = 7.158270687193941D-01 |
---|
| 785 | Y(2) = 9.185534764557338D-02 |
---|
| 786 | Y(3) = 2.841637457458413D+01 |
---|
| 787 | GOTO 310 |
---|
| 788 | |
---|
| 789 | ! PROBLEM D3 |
---|
| 790 | 170 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 |
---|
| 797 | 180 Y(1) = 5.976546980673215D-01 |
---|
| 798 | Y(2) = 1.402343408546138D+00 |
---|
| 799 | Y(3) = -1.893386540441913D-06 |
---|
| 800 | GOTO 310 |
---|
| 801 | |
---|
| 802 | ! PROBLEM D5 |
---|
| 803 | 190 Y(1) = -9.916420698713913D-01 |
---|
| 804 | Y(2) = 9.833363588544478D-01 |
---|
| 805 | GOTO 310 |
---|
| 806 | |
---|
| 807 | ! PROBLEM D6 |
---|
| 808 | 200 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 |
---|
| 816 | 210 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 |
---|
| 823 | 220 Y(1) = -1.158701266031984D+00 |
---|
| 824 | Y(2) = 4.304698089780476D-01 |
---|
| 825 | GOTO 310 |
---|
| 826 | |
---|
| 827 | ! PROBLEM E3 |
---|
| 828 | 230 Y(1) = 4.253052197643089D-03 |
---|
| 829 | Y(2) = 5.317019548450387D-03 |
---|
| 830 | Y(3) = 2.627647748753926D+01 |
---|
| 831 | GOTO 310 |
---|
| 832 | |
---|
| 833 | ! PROBLEM E4 |
---|
| 834 | 240 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 |
---|
| 841 | 250 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 |
---|
| 850 | 260 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 |
---|
| 857 | 270 Y(1) = 3.912699122292088D-01 |
---|
| 858 | Y(2) = 1.329964166084866D-03 |
---|
| 859 | GOTO 310 |
---|
| 860 | |
---|
| 861 | ! PROBLEM F3 |
---|
| 862 | 280 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 |
---|
| 870 | 290 Y(1) = 4.418303324022590D+00 |
---|
| 871 | Y(2) = 1.290244712916425D+00 |
---|
| 872 | Y(3) = 3.019282584050490D+00 |
---|
| 873 | GOTO 310 |
---|
| 874 | |
---|
| 875 | ! PROBLEM F5 |
---|
| 876 | 300 Y(1) = 1.713564284690712D-07 |
---|
| 877 | Y(2) = 3.713563071160676D-03 |
---|
| 878 | Y(3) = 6.189271785267793D-03 |
---|
| 879 | Y(4) = 9.545143571530929D-06 |
---|
| 880 | 310 CONTINUE |
---|
| 881 | IF (IWT<0) GOTO 320 |
---|
| 882 | DO I = 1, N |
---|
| 883 | Y(I) = Y(I)/W(I) |
---|
| 884 | END DO |
---|
| 885 | 320 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 |
---|
| 924 | 10 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 |
---|
| 934 | 20 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 |
---|
| 941 | 30 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 |
---|
| 949 | 40 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 |
---|
| 956 | 50 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 |
---|
| 964 | 60 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 |
---|
| 971 | 70 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 |
---|
| 983 | 80 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 |
---|
| 990 | 90 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 |
---|
| 999 | 100 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 |
---|
| 1005 | 110 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 |
---|
| 1011 | 120 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 |
---|
| 1018 | 130 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 |
---|
| 1024 | 140 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 |
---|
| 1030 | 150 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 |
---|
| 1038 | 160 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 |
---|
| 1046 | 170 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 |
---|
| 1051 | 180 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 |
---|
| 1057 | 190 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 |
---|
| 1080 | 200 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 |
---|
| 1090 | 210 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 |
---|
| 1098 | 220 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 |
---|
| 1103 | 230 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 |
---|
| 1111 | 240 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 |
---|
| 1121 | 250 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) |
---|
| 1125 | 260 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 |
---|
| 1131 | 270 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 |
---|
| 1171 | 10 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 |
---|
| 1181 | 20 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 |
---|
| 1191 | 30 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 |
---|
| 1206 | 40 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 |
---|
| 1222 | 50 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 |
---|
| 1233 | 60 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 |
---|
| 1247 | 70 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 |
---|
| 1264 | 80 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 |
---|
| 1280 | 90 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 |
---|
| 1298 | 100 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 |
---|
| 1310 | 110 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 |
---|
| 1322 | 120 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 |
---|
| 1342 | 130 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 |
---|
| 1354 | 140 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 |
---|
| 1364 | 150 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 |
---|
| 1378 | 160 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 |
---|
| 1394 | 170 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 |
---|
| 1401 | 180 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 |
---|
| 1413 | 190 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 |
---|
| 1451 | 200 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 |
---|
| 1471 | 210 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 |
---|
| 1492 | 220 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 |
---|
| 1499 | 230 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 |
---|
| 1527 | 240 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 |
---|
| 1543 | 250 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 |
---|
| 1559 | 260 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 |
---|
| 1568 | 270 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 | |
---|
| 2026 | 20 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. |
---|
| 2160 | 90000 FORMAT (' Class/Problem = A',I1) |
---|
| 2161 | 90001 FORMAT (' Class/Problem = B',I1) |
---|
| 2162 | 90002 FORMAT (' Class/Problem = C',I1) |
---|
| 2163 | 90003 FORMAT (' Class/Problem = D',I1) |
---|
| 2164 | 90004 FORMAT (' Class/Problem = E',I1) |
---|
| 2165 | 90005 FORMAT (' Class/Problem = F',I1) |
---|
| 2166 | 90006 FORMAT (' An error occurred in VODE_F90. ISTATE = ',I3) |
---|
| 2167 | 90007 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)) |
---|
| 2172 | 90008 FORMAT (' Computed and reference solutions and absolute' & |
---|
| 2173 | ' errors follow:',/,(I3,3D15.5)) |
---|
| 2174 | 90009 FORMAT (' Steps = ',I10,' f-s = ',I10,' J-s = ',I10) |
---|
| 2175 | 90010 FORMAT (' ____________________________________________________________') |
---|
| 2176 | 90011 FORMAT (/,' Errors occurred. (Search on ISTATE.)') |
---|
| 2177 | 90012 FORMAT (/,' No errors occurred.') |
---|
| 2178 | 90013 FORMAT (' For ITOL = ', I3,' Dense max. observed error = ',D15.5) |
---|
| 2179 | 90014 FORMAT (I5,6X,I5,3D15.5) |
---|
| 2180 | 90015 FORMAT (/,'Tolerance',3X,'Problem',4X,'Dense',10X,'Sparse',9X,'Banded') |
---|
| 2181 | 90016 FORMAT (//,20X,' Results for ITOL = ',I3) |
---|
| 2182 | 90017 FORMAT (' For ITOL = ', I3,' Sparse max. observed error = ',D15.5) |
---|
| 2183 | 90018 FORMAT (' For ITOL = ', I3,' Banded max. observed error = ',D15.5) |
---|
| 2184 | 90019 FORMAT (' The dense and banded results are not identical.') |
---|
| 2185 | 90020 FORMAT (/,' Maximum Errors for All Problems Follow:') |
---|
| 2186 | 90021 FORMAT (/,10X,' Maximum Errors for Individual Problems Follow:') |
---|
| 2187 | 90022 FORMAT (' The original VODE Jacobian was used.') |
---|
| 2188 | 90023 FORMAT (' The JACSP Jacobian was used.') |
---|
| 2189 | 90024 FORMAT (/,' Integration Statistics Follow:') |
---|
| 2190 | 90025 FORMAT (/,' Total number of solutions performed = ',I5) |
---|
| 2191 | 90026 FORMAT (I5,' Steps: ',3I10) |
---|
| 2192 | 90027 FORMAT (5X,' Jacs: ',3I10) |
---|
| 2193 | 90028 FORMAT (5X,' Ders: ',3I10) |
---|
| 2194 | 90029 FORMAT (/,' Problem',12X,'Dense',4X,'Sparse',4X,'Banded') |
---|
| 2195 | 90030 FORMAT (' The ODEs were scaled.') |
---|
| 2196 | 90031 FORMAT (' The ODEs were not scaled.') |
---|
| 2197 | 90032 FORMAT (//,' To locate the summary results, search on the phrase',/ & |
---|
| 2198 | ' Summary Statistics Follow in the output file', & |
---|
| 2199 | ' (stiffoptions.dat).',//) |
---|
| 2200 | 90033 FORMAT (/,' Tolerance', 10X, 'Total Execution Time') |
---|
| 2201 | 90034 FORMAT ((4X,I5,10X,I10)) |
---|
| 2202 | 90035 FORMAT (/,' Problem No. = ',I3,/,' ITOL = ',I3,/,' Scaled ODEs', & |
---|
| 2203 | 13X,'Absolute Errors') |
---|
| 2204 | 90036 FORMAT (' Problem No. = ',I3,/,' ITOL = ',I3,/,' Unscaled ODEs', & |
---|
| 2205 | 18X,'Absolute Errors') |
---|
| 2206 | 90037 FORMAT (' Dense : JACSP - VODE: ',2D15.5) |
---|
| 2207 | 90038 FORMAT (' Sparse: JACSP - VODE: ',2D15.5) |
---|
| 2208 | 90039 FORMAT (' Banded: JACSP - VODE: ',2D15.5) |
---|
| 2209 | 90040 FORMAT (' Nonnegativity will be enforced for this problem.') |
---|
| 2210 | 90041 FORMAT (' Dense SDJ: ',6I7) |
---|
| 2211 | 90042 FORMAT (' Sparse SDJ: ',6I7) |
---|
| 2212 | 90043 FORMAT (' Banded SDJ: ',6I7) |
---|
| 2213 | 90044 FORMAT (19X,'J',6X,'V',6X,'J',6X,'V',6X,'J',6X,'V') |
---|
| 2214 | 90045 FORMAT (///,20X,'Summary Statistics Follow') |
---|
| 2215 | |
---|
| 2216 | STOP |
---|
| 2217 | END PROGRAM DEMOSTIFF |
---|