GNU Octave  4.0.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
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 IOWND, IOWNS,
46  1 icf, ierpj, iersl, jcur, jstart, kflag, l,
47  2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
48  3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
49  REAL ROWNS,
50  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
51  COMMON /sls001/ rowns(209),
52  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
53  2 iownd(6), iowns(6),
54  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
55  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
56  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
57  INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
58  REAL C, R, S, TP
59  CHARACTER*80 MSG
60 C
61 C***FIRST EXECUTABLE STATEMENT SINTDY
62  iflag = 0
63  IF (k .LT. 0 .OR. k .GT. nq) go to 80
64  tp = tn - hu - 100.0e0*uround*sign(abs(tn) + abs(hu), hu)
65  IF ((t-tp)*(t-tn) .GT. 0.0e0) go to 90
66 C
67  s = (t - tn)/h
68  ic = 1
69  IF (k .EQ. 0) go to 15
70  jj1 = l - k
71  DO 10 jj = jj1,nq
72  10 ic = ic*jj
73  15 c = ic
74  DO 20 i = 1,n
75  20 dky(i) = c*yh(i,l)
76  IF (k .EQ. nq) go to 55
77  jb2 = nq - k
78  DO 50 jb = 1,jb2
79  j = nq - jb
80  jp1 = j + 1
81  ic = 1
82  IF (k .EQ. 0) go to 35
83  jj1 = jp1 - k
84  DO 30 jj = jj1,j
85  30 ic = ic*jj
86  35 c = ic
87  DO 40 i = 1,n
88  40 dky(i) = c*yh(i,jp1) + s*dky(i)
89  50 CONTINUE
90  IF (k .EQ. 0) RETURN
91  55 r = h**(-k)
92  DO 60 i = 1,n
93  60 dky(i) = r*dky(i)
94  RETURN
95 C
96  80 CALL xerrwd('SINTDY- K (=I1) illegal ',
97  1 30, 51, 0, 1, k, 0, 0, 0.0e0, 0.0e0)
98  iflag = -1
99  RETURN
100  90 CALL xerrwd('SINTDY- T (=R1) illegal ',
101  1 30, 52, 0, 0, 0, 0, 1, t, 0.0e0)
102  CALL xerrwd(
103  1 ' T not in interval TCUR - HU (= R1) to TCUR (=R2) ',
104  1 60, 52, 0, 0, 0, 0, 2, tp, tn)
105  iflag = -2
106  RETURN
107 C----------------------- END OF SUBROUTINE SINTDY ----------------------
108  END
std::string dimension(void) const
subroutine xerrwd(MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
Definition: xerrwd.f:3
subroutine sintdy(T, K, YH, NYH, DKY, IFLAG)
Definition: sintdy.f:1
T abs(T x)
Definition: pr-output.cc:3062