GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
sintdy.f
Go to the documentation of this file.
1  SUBROUTINE sintdy (T, K, YH, NYH, DKY, IFLAG)
2 C***BEGIN PROLOGUE SINTDY
3 C***SUBSIDIARY
4 C***PURPOSE Interpolate solution derivatives.
5 C***TYPE SINGLE PRECISION (SINTDY-S, DINTDY-D)
6 C***AUTHOR Hindmarsh, Alan C., (LLNL)
7 C***DESCRIPTION
8 C
9 C SINTDY computes interpolated values of the K-th derivative of the
10 C dependent variable vector y, and stores it in DKY. This routine
11 C is called within the package with K = 0 and T = TOUT, but may
12 C also be called by the user for any K up to the current order.
13 C (See detailed instructions in the usage documentation.)
14 C
15 C The computed values in DKY are gotten by interpolation using the
16 C Nordsieck history array YH. This array corresponds uniquely to a
17 C vector-valued polynomial of degree NQCUR or less, and DKY is set
18 C to the K-th derivative of this polynomial at T.
19 C The formula for DKY is:
20 C q
21 C DKY(i) = sum c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1)
22 C j=K
23 C where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR.
24 C The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are
25 C communicated by COMMON. The above sum is done in reverse order.
26 C IFLAG is returned negative if either K or T is out of bounds.
27 C
28 C***SEE ALSO SLSODE
29 C***ROUTINES CALLED XERRWV
30 C***COMMON BLOCKS SLS001
31 C***REVISION HISTORY (YYMMDD)
32 C 791129 DATE WRITTEN
33 C 890501 Modified prologue to SLATEC/LDOC format. (FNF)
34 C 890503 Minor cosmetic changes. (FNF)
35 C 930809 Renamed to allow single/double precision versions. (ACH)
36 C 010412 Reduced size of Common block /SLS001/. (ACH)
37 C 031105 Restored 'own' variables to Common block /SLS001/, to
38 C enable interrupt/restart feature. (ACH)
39 C 050427 Corrected roundoff decrement in TP. (ACH)
40 C***END PROLOGUE SINTDY
41 C**End
42  INTEGER K, NYH, IFLAG
43  REAL T, YH, DKY
44  dimension yh(nyh,*), dky(*)
45  INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, CNYH,
46  1 IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
47  1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
48  2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
49  3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
50  REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
51  1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
52  COMMON /sls001/ conit, crate, el(13), elco(13,12),
53  1 hold, rmax, tesco(3,12),
54  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
55  2 init, mxstep, mxhnil, nhnil, nslast, cnyh,
56  3 ialth, ipup, lmax, meo, nqnyh, nslp,
57  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
58  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
59  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
60  INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
61  REAL C, R, S, TP
62  CHARACTER*80 MSG
63 C
64 C***FIRST EXECUTABLE STATEMENT SINTDY
65  iflag = 0
66  IF (k .LT. 0 .OR. k .GT. nq) GO TO 80
67  tp = tn - hu - 100.0e0*uround*sign(abs(tn) + abs(hu), hu)
68  IF ((t-tp)*(t-tn) .GT. 0.0e0) GO TO 90
69 C
70  s = (t - tn)/h
71  ic = 1
72  IF (k .EQ. 0) GO TO 15
73  jj1 = l - k
74  DO 10 jj = jj1,nq
75  10 ic = ic*jj
76  15 c = ic
77  DO 20 i = 1,n
78  20 dky(i) = c*yh(i,l)
79  IF (k .EQ. nq) GO TO 55
80  jb2 = nq - k
81  DO 50 jb = 1,jb2
82  j = nq - jb
83  jp1 = j + 1
84  ic = 1
85  IF (k .EQ. 0) GO TO 35
86  jj1 = jp1 - k
87  DO 30 jj = jj1,j
88  30 ic = ic*jj
89  35 c = ic
90  DO 40 i = 1,n
91  40 dky(i) = c*yh(i,jp1) + s*dky(i)
92  50 CONTINUE
93  IF (k .EQ. 0) RETURN
94  55 r = h**(-k)
95  DO 60 i = 1,n
96  60 dky(i) = r*dky(i)
97  RETURN
98 C
99  80 CALL xerrwd('SINTDY- K (=I1) illegal ',
100  1 30, 51, 0, 1, k, 0, 0, 0.0d0, 0.0d0)
101  iflag = -1
102  RETURN
103  90 CALL xerrwd('SINTDY- T (=R1) illegal ',
104  1 30, 52, 0, 0, 0, 0, 1, dble(t), 0.0d0)
105  CALL xerrwd(
106  1 ' T not in interval TCUR - HU (= R1) to TCUR (=R2) ',
107  1 60, 52, 0, 0, 0, 0, 2, dble(tp), dble(tn))
108  iflag = -2
109  RETURN
110 C----------------------- END OF SUBROUTINE SINTDY ----------------------
111  END
static T abs(T x)
Definition: pr-output.cc:1678
subroutine sintdy(T, K, YH, NYH, DKY, IFLAG)
Definition: sintdy.f:2
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
Definition: xerrwd.f:4