sintdy.f

Go to the documentation of this file.
00001       SUBROUTINE SINTDY (T, K, YH, NYH, DKY, IFLAG)
00002 C***BEGIN PROLOGUE  SINTDY
00003 C***SUBSIDIARY
00004 C***PURPOSE  Interpolate solution derivatives.
00005 C***TYPE      SINGLE PRECISION (SINTDY-S, DINTDY-D)
00006 C***AUTHOR  Hindmarsh, Alan C., (LLNL)
00007 C***DESCRIPTION
00008 C
00009 C  SINTDY computes interpolated values of the K-th derivative of the
00010 C  dependent variable vector y, and stores it in DKY.  This routine
00011 C  is called within the package with K = 0 and T = TOUT, but may
00012 C  also be called by the user for any K up to the current order.
00013 C  (See detailed instructions in the usage documentation.)
00014 C
00015 C  The computed values in DKY are gotten by interpolation using the
00016 C  Nordsieck history array YH.  This array corresponds uniquely to a
00017 C  vector-valued polynomial of degree NQCUR or less, and DKY is set
00018 C  to the K-th derivative of this polynomial at T.
00019 C  The formula for DKY is:
00020 C               q
00021 C   DKY(i)  =  sum  c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1)
00022 C              j=K
00023 C  where  c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR.
00024 C  The quantities  nq = NQCUR, l = nq+1, N = NEQ, tn, and h are
00025 C  communicated by COMMON.  The above sum is done in reverse order.
00026 C  IFLAG is returned negative if either K or T is out of bounds.
00027 C
00028 C***SEE ALSO  SLSODE
00029 C***ROUTINES CALLED  XERRWV
00030 C***COMMON BLOCKS    SLS001
00031 C***REVISION HISTORY  (YYMMDD)
00032 C   791129  DATE WRITTEN
00033 C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
00034 C   890503  Minor cosmetic changes.  (FNF)
00035 C   930809  Renamed to allow single/double precision versions. (ACH)
00036 C   010412  Reduced size of Common block /SLS001/. (ACH)
00037 C   031105  Restored 'own' variables to Common block /SLS001/, to
00038 C           enable interrupt/restart feature. (ACH)
00039 C   050427  Corrected roundoff decrement in TP. (ACH)
00040 C***END PROLOGUE  SINTDY
00041 C**End
00042       INTEGER K, NYH, IFLAG
00043       REAL T, YH, DKY
00044       DIMENSION YH(NYH,*), DKY(*)
00045       INTEGER IOWND, IOWNS,
00046      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
00047      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
00048      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
00049       REAL ROWNS,
00050      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
00051       COMMON /SLS001/ ROWNS(209),
00052      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
00053      2   IOWND(6), IOWNS(6),
00054      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
00055      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
00056      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
00057       INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
00058       REAL C, R, S, TP
00059       CHARACTER*80 MSG
00060 C
00061 C***FIRST EXECUTABLE STATEMENT  SINTDY
00062       IFLAG = 0
00063       IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80
00064       TP = TN - HU -  100.0E0*UROUND*SIGN(ABS(TN) + ABS(HU), HU)
00065       IF ((T-TP)*(T-TN) .GT. 0.0E0) GO TO 90
00066 C
00067       S = (T - TN)/H
00068       IC = 1
00069       IF (K .EQ. 0) GO TO 15
00070       JJ1 = L - K
00071       DO 10 JJ = JJ1,NQ
00072  10     IC = IC*JJ
00073  15   C = IC
00074       DO 20 I = 1,N
00075  20     DKY(I) = C*YH(I,L)
00076       IF (K .EQ. NQ) GO TO 55
00077       JB2 = NQ - K
00078       DO 50 JB = 1,JB2
00079         J = NQ - JB
00080         JP1 = J + 1
00081         IC = 1
00082         IF (K .EQ. 0) GO TO 35
00083         JJ1 = JP1 - K
00084         DO 30 JJ = JJ1,J
00085  30       IC = IC*JJ
00086  35     C = IC
00087         DO 40 I = 1,N
00088  40       DKY(I) = C*YH(I,JP1) + S*DKY(I)
00089  50     CONTINUE
00090       IF (K .EQ. 0) RETURN
00091  55   R = H**(-K)
00092       DO 60 I = 1,N
00093  60     DKY(I) = R*DKY(I)
00094       RETURN
00095 C
00096  80   CALL XERRWD('SINTDY-  K (=I1) illegal      ', 
00097      1     30, 51, 0, 1, K, 0, 0, 0.0E0, 0.0E0)
00098       IFLAG = -1
00099       RETURN
00100  90   CALL XERRWD('SINTDY-  T (=R1) illegal      ', 
00101      1     30, 52, 0, 0, 0, 0, 1, T, 0.0E0)
00102       CALL XERRWD(
00103      1   '      T not in interval TCUR - HU (= R1) to TCUR (=R2)      ',
00104      1    60, 52, 0, 0, 0, 0, 2, TP, TN)
00105       IFLAG = -2
00106       RETURN
00107 C----------------------- END OF SUBROUTINE SINTDY ----------------------
00108       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines