sintdy.f

Go to the documentation of this file.

      SUBROUTINE sintdy (T, K, YH, NYH, DKY, IFLAG)
C***BEGIN PROLOGUE  SINTDY
C***SUBSIDIARY
C***PURPOSE  Interpolate solution derivatives.
C***TYPE      SINGLE PRECISION (SINTDY-S, DINTDY-D)
C***AUTHOR  Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C
C  SINTDY computes interpolated values of the K-th derivative of the
C  dependent variable vector y, and stores it in DKY.  This routine
C  is called within the package with K = 0 and T = TOUT, but may
C  also be called by the user for any K up to the current order.
C  (See detailed instructions in the usage documentation.)
C
C  The computed values in DKY are gotten by interpolation using the
C  Nordsieck history array YH.  This array corresponds uniquely to a
C  vector-valued polynomial of degree NQCUR or less, and DKY is set
C  to the K-th derivative of this polynomial at T.
C  The formula for DKY is:
C               q
C   DKY(i)  =  sum  c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1)
C              j=K
C  where  c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR.
C  The quantities  nq = NQCUR, l = nq+1, N = NEQ, tn, and h are
C  communicated by COMMON.  The above sum is done in reverse order.
C  IFLAG is returned negative if either K or T is out of bounds.
C
C***SEE ALSO  SLSODE
C***ROUTINES CALLED  XERRWV
C***COMMON BLOCKS    SLS001
C***REVISION HISTORY  (YYMMDD)
C   791129  DATE WRITTEN
C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
C   890503  Minor cosmetic changes.  (FNF)
C   930809  Renamed to allow single/double precision versions. (ACH)
C   010412  Reduced size of Common block /SLS001/. (ACH)
C   031105  Restored 'own' variables to Common block /SLS001/, to
C           enable interrupt/restart feature. (ACH)
C   050427  Corrected roundoff decrement in TP. (ACH)
C***END PROLOGUE  SINTDY
C**End
      INTEGER K, NYH, IFLAG
      REAL T, YH, DKY
      dimension yh(nyh,*), dky(*)
      INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, CNYH,
     1   IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
     1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
     2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
     3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
      REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
     1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
      COMMON /sls001/ conit, crate, el(13), elco(13,12),
     1   hold, rmax, tesco(3,12),
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     2   init, mxstep, mxhnil, nhnil, nslast, cnyh,
     3   ialth, ipup, lmax, meo, nqnyh, nslp,
     3   icf, ierpj, iersl, jcur, jstart, kflag, l,
     4   lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
      REAL C, R, S, TP
      CHARACTER*80 MSG
C
C***FIRST EXECUTABLE STATEMENT  SINTDY
      iflag = 0
      IF (k .LT. 0 .OR. k .GT. nq) GO TO 80
      tp = tn - hu -  100.0e0*uround*sign(abs(tn) + abs(hu), hu)
      IF ((t-tp)*(t-tn) .GT. 0.0e0) GO TO 90
C
      s = (t - tn)/h
      ic = 1
      IF (k .EQ. 0) GO TO 15
      jj1 = l - k
      DO 10 jj = jj1,nq
 10     ic = ic*jj
 15   c = ic
      DO 20 i = 1,n
 20     dky(i) = c*yh(i,l)
      IF (k .EQ. nq) GO TO 55
      jb2 = nq - k
      DO 50 jb = 1,jb2
        j = nq - jb
        jp1 = j + 1
        ic = 1
        IF (k .EQ. 0) GO TO 35
        jj1 = jp1 - k
        DO 30 jj = jj1,j
 30       ic = ic*jj
 35     c = ic
        DO 40 i = 1,n
 40       dky(i) = c*yh(i,jp1) + s*dky(i)
 50     CONTINUE
      IF (k .EQ. 0) RETURN
 55   r = h**(-k)
      DO 60 i = 1,n
 60     dky(i) = r*dky(i)
      RETURN
C
 80   CALL xerrwd('SINTDY-  K (=I1) illegal      ',
     1     30, 51, 0, 1, k, 0, 0, 0.0e0, 0.0e0)
      iflag = -1
      RETURN
 90   CALL xerrwd('SINTDY-  T (=R1) illegal      ',
     1     30, 52, 0, 0, 0, 0, 1, t, 0.0e0)
      CALL xerrwd(
     1   '      T not in interval TCUR - HU (= R1) to TCUR (=R2)      ',
     1    60, 52, 0, 0, 0, 0, 2, tp, tn)
      iflag = -2
      RETURN
C----------------------- END OF SUBROUTINE SINTDY ----------------------
      END