Mercurial > octave-nkf

annotate libcruft/odepack/sintdy.f @ 7793:96ba591be50f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add some more support for single precision to libcruft functions
author David Bateman <dbateman@free.fr>
date Sun, 11 May 2008 22:51:50 +0200
parents
children 2b458dfe31ae
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7793
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
1 SUBROUTINE SINTDY (T, K, YH, NYH, DKY, IFLAG)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
2 C***BEGIN PROLOGUE SINTDY
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
3 C***SUBSIDIARY
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
4 C***PURPOSE Interpolate solution derivatives.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
5 C***TYPE SINGLE PRECISION (SINTDY-S, DINTDY-D)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
6 C***AUTHOR Hindmarsh, Alan C., (LLNL)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
7 C***DESCRIPTION
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
8 C
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
9 C SINTDY computes interpolated values of the K-th derivative of the
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
10 C dependent variable vector y, and stores it in DKY. This routine
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
11 C is called within the package with K = 0 and T = TOUT, but may
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
12 C also be called by the user for any K up to the current order.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
13 C (See detailed instructions in the usage documentation.)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
14 C
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
15 C The computed values in DKY are gotten by interpolation using the
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
16 C Nordsieck history array YH. This array corresponds uniquely to a
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
17 C vector-valued polynomial of degree NQCUR or less, and DKY is set
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
18 C to the K-th derivative of this polynomial at T.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
19 C The formula for DKY is:
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
20 C q
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
21 C DKY(i) = sum c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
22 C j=K
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
23 C where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
24 C The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
25 C communicated by COMMON. The above sum is done in reverse order.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
26 C IFLAG is returned negative if either K or T is out of bounds.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
27 C
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
28 C***SEE ALSO SLSODE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
29 C***ROUTINES CALLED XERRWV
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
30 C***COMMON BLOCKS SLS001
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
31 C***REVISION HISTORY (YYMMDD)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
32 C 791129 DATE WRITTEN
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
33 C 890501 Modified prologue to SLATEC/LDOC format. (FNF)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
34 C 890503 Minor cosmetic changes. (FNF)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
35 C 930809 Renamed to allow single/double precision versions. (ACH)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
36 C 010412 Reduced size of Common block /SLS001/. (ACH)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
37 C 031105 Restored 'own' variables to Common block /SLS001/, to
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
38 C enable interrupt/restart feature. (ACH)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
39 C 050427 Corrected roundoff decrement in TP. (ACH)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
40 C***END PROLOGUE SINTDY
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
41 C**End
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
42 INTEGER K, NYH, IFLAG
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
43 REAL T, YH, DKY
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
44 DIMENSION YH(NYH,*), DKY(*)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
45 INTEGER IOWND, IOWNS,
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
46 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
47 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
48 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
49 REAL ROWNS,
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
50 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
51 COMMON /SLS001/ ROWNS(209),
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
52 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
53 2 IOWND(6), IOWNS(6),
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
54 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
55 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
56 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
57 INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
58 REAL C, R, S, TP
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
59 CHARACTER*80 MSG
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
60 C
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
61 C***FIRST EXECUTABLE STATEMENT SINTDY
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
62 IFLAG = 0
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
63 IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
64 TP = TN - HU - 100.0E0*UROUND*SIGN(ABS(TN) + ABS(HU), HU)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
65 IF ((T-TP)*(T-TN) .GT. 0.0E0) GO TO 90
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
66 C
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
67 S = (T - TN)/H
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
68 IC = 1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
69 IF (K .EQ. 0) GO TO 15
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
70 JJ1 = L - K
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
71 DO 10 JJ = JJ1,NQ
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
72 10 IC = IC*JJ
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
73 15 C = IC
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
74 DO 20 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
75 20 DKY(I) = C*YH(I,L)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
76 IF (K .EQ. NQ) GO TO 55
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
77 JB2 = NQ - K
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
78 DO 50 JB = 1,JB2
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
79 J = NQ - JB
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
80 JP1 = J + 1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
81 IC = 1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
82 IF (K .EQ. 0) GO TO 35
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
83 JJ1 = JP1 - K
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
84 DO 30 JJ = JJ1,J
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
85 30 IC = IC*JJ
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
86 35 C = IC
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
87 DO 40 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
88 40 DKY(I) = C*YH(I,JP1) + S*DKY(I)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
89 50 CONTINUE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
90 IF (K .EQ. 0) RETURN
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
91 55 R = H**(-K)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
92 DO 60 I = 1,N
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
93 60 DKY(I) = R*DKY(I)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
94 RETURN
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
95 C
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
96 80 MSG = 'SINTDY- K (=I1) illegal '
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
97 CALL XERRWV (MSG, 30, 51, 0, 1, K, 0, 0, 0.0E0, 0.0E0)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
98 IFLAG = -1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
99 RETURN
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
100 90 MSG = 'SINTDY- T (=R1) illegal '
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
101 CALL XERRWV (MSG, 30, 52, 0, 0, 0, 0, 1, T, 0.0E0)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
102 MSG=' T not in interval TCUR - HU (= R1) to TCUR (=R2) '
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
103 CALL XERRWV (MSG, 60, 52, 0, 0, 0, 0, 2, TP, TN)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
104 IFLAG = -2
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
105 RETURN
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
106 C----------------------- END OF SUBROUTINE SINTDY ----------------------
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
107 END