| 1 |
PROGRAM DEMO |
| 2 |
C |
| 3 |
C*********************************************************************** |
| 4 |
C |
| 5 |
C SDASSL routines |
| 6 |
C by Linda Petzold, Andreas Kroener, Wolfgang Marquardt |
| 7 |
C Created: 15/03/83 |
| 8 |
C Version: 1.1 Rev: 1989/12/11 |
| 9 |
C Date last modified: 1994/09/02 |
| 10 |
C |
| 11 |
C This file is part of the SDASSL differential/algebraic system solver. |
| 12 |
C |
| 13 |
C Copyright (C) 1983, 1989, 1994 Linda Petzold, Andreas Kroener, |
| 14 |
C Wolfgang Marquardt |
| 15 |
C |
| 16 |
C The SDASSL differential/algebraic system solver is free software; |
| 17 |
C you can redistribute it and/or modify it under the terms of the GNU |
| 18 |
C General Public License as published by the Free Software Foundation; |
| 19 |
C either version 2 of the License, or (at your option) any later version. |
| 20 |
C |
| 21 |
C The SDASSL system solver is distributed in hope that it will be |
| 22 |
C useful, but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 23 |
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 24 |
C General Public License for more details. |
| 25 |
C |
| 26 |
C You should have received a copy of the GNU General Public License along |
| 27 |
C with the program; if not, write to the Free Software Foundation, Inc., |
| 28 |
C 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING. |
| 29 |
C |
| 30 |
C*********************************************************************** |
| 31 |
C |
| 32 |
C Author: A. Kroener 25 - June - 1990 |
| 33 |
C Institut f"ur Systemdynamik und Regelungstechnik |
| 34 |
C Universit"at Stuttgart |
| 35 |
C Pfaffenwaldring 9 |
| 36 |
C D - 7000 Stuttgart 80 |
| 37 |
C |
| 38 |
C |
| 39 |
C Datum Bearbeiter Aenderungen |
| 40 |
C |
| 41 |
C*********************************************************************** |
| 42 |
C |
| 43 |
C Aufgabe des Programms: |
| 44 |
C ---------------------- |
| 45 |
C |
| 46 |
C Driver routine for the sparse DA-integrator SDASSL. |
| 47 |
C For further details see the updated description |
| 48 |
C SDASSL_INFO.TXT. |
| 49 |
C The values of the control arrau INFO(1:11) are input from the |
| 50 |
C data file INFO.DAT. |
| 51 |
C In the data file CONTRL.DAT the user provides the required scalar |
| 52 |
C relative (RTOL) and absolute (ATOL) error tolerances and the |
| 53 |
C number of examples to be calculated. Two simple examples are |
| 54 |
C included. |
| 55 |
C |
| 56 |
C Each of the examples consists of three routines: |
| 57 |
C INIT# reads from the file INIT#.DAT the system dimension, |
| 58 |
C initial and final time and an proposed initial time |
| 59 |
C step. |
| 60 |
C Initial conditions for the states Y and their |
| 61 |
C derivatives YP are also provided in this file. |
| 62 |
C |
| 63 |
C PTN# provides the sparsity pattern of the merged jacobian |
| 64 |
C dF/dY and dF/dYP in sparse colum vector notation. |
| 65 |
C |
| 66 |
C FUN# evaluates the system equations 0=F(Y,YP,T) |
| 67 |
C |
| 68 |
C*********************************************************************** |
| 69 |
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
| 70 |
IMPLICIT INTEGER (I-N) |
| 71 |
C |
| 72 |
PARAMETER (NGL=20,NII2 = 1) |
| 73 |
PARAMETER (NJAC = NGL*NGL) |
| 74 |
PARAMETER (NLUD = 2*NJAC) |
| 75 |
PARAMETER (MLUD = NLUD) |
| 76 |
PARAMETER (LIW = 20 + NGL*(2+13/NII2) + (NJAC+NLUD+MLUD)/NII2) |
| 77 |
PARAMETER (LRW = 40 + 10*NGL + NJAC + NLUD) |
| 78 |
C |
| 79 |
EXTERNAL PTN1, FUN1, INIT1 |
| 80 |
EXTERNAL PTN2, FUN2, INIT2 |
| 81 |
C |
| 82 |
DIMENSION Y(NGL), YP(NGL) |
| 83 |
DIMENSION RWS(LRW) |
| 84 |
DIMENSION IWS(LIW) |
| 85 |
C |
| 86 |
COMMON / ANVE / IDUM1,IDUM2,IDUM3,LENB |
| 87 |
COMMON / JACV / LENFX,LENBX,IDUM4,IDUM5 |
| 88 |
C |
| 89 |
NG = NGL |
| 90 |
NOUT = 6 |
| 91 |
IU1 = 11 |
| 92 |
IU2 = 12 |
| 93 |
C |
| 94 |
OPEN (IU1, FILE='CONTRL.DAT', STATUS='OLD', FORM='FORMATTED', |
| 95 |
$ ACCESS='SEQUENTIAL') |
| 96 |
C |
| 97 |
C ---- read tolerances for integration |
| 98 |
C |
| 99 |
READ (IU1,*) RTOL, ATOL |
| 100 |
C |
| 101 |
C ---- loop over several examples |
| 102 |
C |
| 103 |
100 CONTINUE |
| 104 |
C |
| 105 |
C ---- read no. of example |
| 106 |
C |
| 107 |
READ (IU1,*,END=200) JBSP |
| 108 |
WRITE (NOUT,*) 'Example No.', JBSP |
| 109 |
C |
| 110 |
C |
| 111 |
N = NGL |
| 112 |
ISTEP = 0 |
| 113 |
C |
| 114 |
IF (JBSP .EQ. 1) THEN |
| 115 |
C |
| 116 |
C ---- initialization and initial values |
| 117 |
C |
| 118 |
CALL INIT1 (N, Y, YP, T0, TE, H0) |
| 119 |
C |
| 120 |
C ---- integration |
| 121 |
C |
| 122 |
CALL INTG (N, Y, YP, IWS, IWS(16), (LIW-15), RWS, LRW, T0, TE, |
| 123 |
1 H0, RTOL, ATOL, PTN1, FUN1) |
| 124 |
ELSE IF (JBSP .EQ. 2) THEN |
| 125 |
CALL INIT2 (N, Y, YP, T0, TE, H0) |
| 126 |
CALL INTG (N, Y, YP, IWS, IWS(16), (LIW-15), RWS, LRW, T0, TE, |
| 127 |
1 H0, RTOL, ATOL, PTN2, FUN2) |
| 128 |
ELSE |
| 129 |
STOP ' Error in no. of example' |
| 130 |
END IF |
| 131 |
C |
| 132 |
GOTO 100 |
| 133 |
C |
| 134 |
200 CLOSE (IU1) |
| 135 |
WRITE (NOUT,'(/5X,A/)') 'DEMO successfully completed' |
| 136 |
STOP |
| 137 |
END |
| 138 |
|
| 139 |
SUBROUTINE INTG (N, Y, YP, INFO, IWORK, LIW, RWORK, LRW, T0, |
| 140 |
1 TFINAL, H0, RTOL, ATOL, PTN, FUN) |
| 141 |
C |
| 142 |
C*********************************************************************** |
| 143 |
C |
| 144 |
C supply parameters for integrator DDASSL |
| 145 |
C |
| 146 |
C*********************************************************************** |
| 147 |
C |
| 148 |
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
| 149 |
C |
| 150 |
DIMENSION Y(N), YP(N) |
| 151 |
DIMENSION INFO(1), IWORK(1), RWORK(1) |
| 152 |
C |
| 153 |
C COMMON / IO / NOUT |
| 154 |
C COMMON / HLP / NN, JBSP |
| 155 |
C |
| 156 |
EXTERNAL PTN, FUN |
| 157 |
C |
| 158 |
OPEN (UNIT=20,FILE='INFO.DAT',STATUS='OLD') |
| 159 |
READ(20,21)(INFO(I),I=1,15) |
| 160 |
21 FORMAT(I3) |
| 161 |
WRITE(6,25)(info(i),I = 1,15) |
| 162 |
25 FORMAT(1x,I3) |
| 163 |
WRITE(6,*) 'info complete' |
| 164 |
CLOSE(20) |
| 165 |
|
| 166 |
T = T0 |
| 167 |
DELTAT = .1 |
| 168 |
TOUT = T0 + DELTAT ! next output time |
| 169 |
RWORK(1) = T0 + DELTAT ! TSTOP |
| 170 |
RWORK(2) = 1.0 ! HMAX |
| 171 |
RWORK(3) = H0 |
| 172 |
C |
| 173 |
IWORK(1) = 4 |
| 174 |
IWORK(2) = 4 |
| 175 |
IWORK(3) = 1 |
| 176 |
ICHAN = 0 |
| 177 |
|
| 178 |
C |
| 179 |
C ---- integration |
| 180 |
C |
| 181 |
10 CONTINUE |
| 182 |
CALL OUT (T, N, Y, YP) |
| 183 |
DO 40 WHILE (T .LE. TFINAL) |
| 184 |
DO 45 WHILE (T .LE. TOUT) |
| 185 |
CALL DDASSL (FUN, N, T, Y, YP, TOUT, INFO, RTOL, ATOL, IDID, |
| 186 |
$ RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC,PTN, |
| 187 |
1 ICHAN) |
| 188 |
C |
| 189 |
WRITE(6,*)idid |
| 190 |
CALL OUT (T, N, Y, YP) |
| 191 |
IF (idid .LT. 0 ) GO TO 60 |
| 192 |
IF (IDID .EQ. -1) THEN |
| 193 |
INFO(1) = 1 |
| 194 |
ICHAN = 0 |
| 195 |
ELSEIF (idid .EQ. 1 .OR. idid .EQ. 3) THEN |
| 196 |
INFO(1) = 1 |
| 197 |
ICHAN = 0 |
| 198 |
ENDIF |
| 199 |
IF ((IDID .EQ. 3) .OR. (IDID .EQ. 2) ) GO TO 55 |
| 200 |
45 END DO |
| 201 |
55 TOUT = t + DELTAT |
| 202 |
RWORK(1) = t + DELTAT |
| 203 |
40 END DO |
| 204 |
60 CONTINUE |
| 205 |
|
| 206 |
RETURN |
| 207 |
END |
| 208 |
|
| 209 |
SUBROUTINE OUT (T, N, Y, YP) |
| 210 |
C |
| 211 |
C*********************************************************************** |
| 212 |
C |
| 213 |
C solution output |
| 214 |
C |
| 215 |
C*********************************************************************** |
| 216 |
C |
| 217 |
IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
| 218 |
IMPLICIT INTEGER (I-N) |
| 219 |
DIMENSION Y(N), YP(N) |
| 220 |
C |
| 221 |
SAVE /ISTEP/ |
| 222 |
DATA ISTEP /1/ |
| 223 |
SAVE /TT/ |
| 224 |
DATA TT /-1./ |
| 225 |
|
| 226 |
IF (TT .GT. T ) ISTEP = 1 |
| 227 |
C |
| 228 |
WRITE (6,*) ISTEP, T |
| 229 |
WRITE (6,*)(Y(I),I=1,N) |
| 230 |
WRITE (6,*)(YP(I),I=1,N) |
| 231 |
ISTEP = ISTEP + 1 |
| 232 |
TT = T |
| 233 |
C |
| 234 |
RETURN |
| 235 |
END |
| 236 |
SUBROUTINE FUN1 (T, Y, YP, DELTA, IRES, RPAR, IPAR) |
| 237 |
C |
| 238 |
DOUBLE PRECISION T, Y(2), YP(2), DELTA(2), RPAR(1), RL |
| 239 |
C |
| 240 |
INTEGER IRES, IPAR(1) |
| 241 |
C |
| 242 |
RL = 8.0D-3 |
| 243 |
DELTA(1) = -YP(1) - RL / (Y(1)*Y(1)) * DEXP (-1.5D0*Y(2)) |
| 244 |
DELTA(2) = - 5.D0 * Y(1) + DEXP (-0.5D0*Y(2)) |
| 245 |
IRES = 0 |
| 246 |
C |
| 247 |
RETURN |
| 248 |
END |
| 249 |
SUBROUTINE FUN2 (T, Y, YP, DELTA, IRES, RPAR, IPAR) |
| 250 |
C |
| 251 |
DOUBLE PRECISION T, Y(3), YP(3), DELTA(3), RPAR(1) |
| 252 |
DOUBLE PRECISION R_1,SK_1,SK_2,SK_3 |
| 253 |
C |
| 254 |
INTEGER IRES, IPAR(1) |
| 255 |
C |
| 256 |
r_1 = .01 |
| 257 |
SK_1 = 100 |
| 258 |
SK_2 = 1 |
| 259 |
SK_3 = 1000 |
| 260 |
C |
| 261 |
DELTA(1) = YP(1) - sk_1 * (Y(2)) |
| 262 |
DELTA(2) = YP(2) - sk_2 * (-Y(2) - Y(3) + R_1 ) |
| 263 |
DELTA(3) = YP(3) - sk_3 * ( Y(1) - Y(3)) |
| 264 |
C |
| 265 |
IRES = 0 |
| 266 |
C |
| 267 |
RETURN |
| 268 |
END |
| 269 |
C**************************************************************************** |
| 270 |
C SUBROUTINE WHICH RETURNS THE MERGED SPARSITY PATTERNS OF THE JACOBIAN |
| 271 |
C********************************************************************* |
| 272 |
|
| 273 |
SUBROUTINE PTN1(IP,M,IRN,LIRN,LDIM,RDUM,IDUM) |
| 274 |
INTEGER M,LIRN,LDIM |
| 275 |
INTEGER IP(M),IRN(LIRN), IDUM(1) |
| 276 |
DOUBLEPRECISION RDUM(1) |
| 277 |
C |
| 278 |
|
| 279 |
IP(1) = 2 |
| 280 |
IP(2) = 2 |
| 281 |
C |
| 282 |
IRN(1) = 1 |
| 283 |
IRN(2) = 2 |
| 284 |
IRN(3) = 1 |
| 285 |
IRN(4) = 2 |
| 286 |
|
| 287 |
LDIM = 4 |
| 288 |
|
| 289 |
END |
| 290 |
C**************************************************************************** |
| 291 |
C SUBROUTINE WHICH RETURNS THE MERGED SPARSITY PATTERNS OF THE JACOBIAN |
| 292 |
C********************************************************************* |
| 293 |
|
| 294 |
SUBROUTINE PTN2(IP,M,IRN,LIRN,LDIM,RDUM,IDUM) |
| 295 |
INTEGER M,LIRN,LDIM |
| 296 |
INTEGER IP(M),IRN(LIRN), IDUM(1) |
| 297 |
DOUBLEPRECISION RDUM(1) |
| 298 |
C |
| 299 |
|
| 300 |
IP(1) = 2 |
| 301 |
IP(2) = 2 |
| 302 |
IP(3) = 2 |
| 303 |
C |
| 304 |
IRN(1) = 1 |
| 305 |
IRN(2) = 3 |
| 306 |
IRN(3) = 1 |
| 307 |
IRN(4) = 2 |
| 308 |
IRN(5) = 2 |
| 309 |
IRN(6) = 3 |
| 310 |
|
| 311 |
LDIM = 6 |
| 312 |
|
| 313 |
END |
| 314 |
C |
| 315 |
C*********************************************************************** |
| 316 |
C |
| 317 |
C initialization |
| 318 |
C |
| 319 |
C*********************************************************************** |
| 320 |
C |
| 321 |
C N : system order (on input at least 100) |
| 322 |
C Y0 : initial values |
| 323 |
C YP0 : initial values of derivatives |
| 324 |
C T0 : initial time |
| 325 |
C TE : end time |
| 326 |
C H0 : initial stepsize (suggestion) |
| 327 |
C |
| 328 |
C************************************************************************ |
| 329 |
C |
| 330 |
SUBROUTINE INIT1 (N, Y0, YP0, T0, TE, H0) |
| 331 |
|
| 332 |
INTEGER N, IUNIT, I |
| 333 |
INTEGER IDUM1,IDUM2,IDUM3,LENB |
| 334 |
INTEGER LENFX,LENBX,IDUM4,IDUM5 |
| 335 |
C |
| 336 |
DOUBLE PRECISION Y0(N), YP0(N) |
| 337 |
DOUBLE PRECISION T0, TE, H0 |
| 338 |
C |
| 339 |
COMMON / ANVE / IDUM1,IDUM2,IDUM3,LENB |
| 340 |
COMMON / JACV / LENFX,LENBX,IDUM4,IDUM5 |
| 341 |
|
| 342 |
LENB = 1 |
| 343 |
LENFX = 4 |
| 344 |
LENBX = 0 |
| 345 |
C |
| 346 |
IUNIT = 21 |
| 347 |
C |
| 348 |
OPEN (IUNIT, FILE='INIT1.DAT', STATUS='OLD', FORM ='FORMATTED', |
| 349 |
1 ACCESS='SEQUENTIAL') |
| 350 |
C |
| 351 |
READ (IUNIT,*) N, T0, TE, H0 |
| 352 |
C |
| 353 |
READ (IUNIT,*) (Y0(I), I=1,N) |
| 354 |
READ (IUNIT,*) (YP0(I), I=1,N) |
| 355 |
|
| 356 |
CLOSE (IUNIT) |
| 357 |
RETURN |
| 358 |
END |
| 359 |
C |
| 360 |
C*********************************************************************** |
| 361 |
C |
| 362 |
C initialization |
| 363 |
C |
| 364 |
C*********************************************************************** |
| 365 |
C |
| 366 |
C N : system order (on input at least 100) |
| 367 |
C Y0 : initial values |
| 368 |
C YP0 : initial values of derivatives |
| 369 |
C T0 : initial time |
| 370 |
C TE : end time |
| 371 |
C H0 : initial stepsize (suggestion) |
| 372 |
C |
| 373 |
C************************************************************************ |
| 374 |
C |
| 375 |
SUBROUTINE INIT2 (N, Y0, YP0, T0, TE, H0) |
| 376 |
|
| 377 |
INTEGER N, IUNIT, I |
| 378 |
INTEGER IDUM1,IDUM2,IDUM3,LENB |
| 379 |
INTEGER LENFX,LENBX,IDUM4,IDUM5 |
| 380 |
C |
| 381 |
DOUBLE PRECISION Y0(N), YP0(N) |
| 382 |
DOUBLE PRECISION T0, TE, H0 |
| 383 |
C |
| 384 |
COMMON / ANVE / IDUM1,IDUM2,IDUM3,LENB |
| 385 |
COMMON / JACV / LENFX,LENBX,IDUM4,IDUM5 |
| 386 |
|
| 387 |
LENB = 3 |
| 388 |
LENFX = 5 |
| 389 |
LENBX = 0 |
| 390 |
C |
| 391 |
IUNIT = 21 |
| 392 |
C |
| 393 |
OPEN (IUNIT, FILE='INIT2.DAT', STATUS='OLD', FORM ='FORMATTED', |
| 394 |
1 ACCESS='SEQUENTIAL') |
| 395 |
C |
| 396 |
READ (IUNIT,*) N, T0, TE, H0 |
| 397 |
C |
| 398 |
READ (IUNIT,*) (Y0(I), I=1,N) |
| 399 |
READ (IUNIT,*) (YP0(I), I=1,N) |
| 400 |
|
| 401 |
CLOSE (IUNIT) |
| 402 |
RETURN |
| 403 |
END |
Parent Directory
| 
Revision Log