share/odepack/fortran/datp.f

   1 *DECK DATP
   2       SUBROUTINE DATP (NEQ, Y, SAVF, P, WGHT, HL0, WK, F, W)
   3       EXTERNAL F
   4       INTEGER NEQ
   5       DOUBLE PRECISION Y, SAVF, P, WGHT, HL0, WK, W
   6       DIMENSION NEQ(*), Y(*), SAVF(*), P(*), WGHT(*), WK(*), W(*)
   7       INTEGER IOWND, IOWNS,
   8      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
   9      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  10      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  11       DOUBLE PRECISION ROWNS,
  12      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  13       COMMON /DLS001/ ROWNS(209),
  14      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  15      2   IOWND(6), IOWNS(6),
  16      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  17      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  18      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  19 C-----------------------------------------------------------------------
  20 C This routine computes the product
  21 C
  22 C              w = (I - hl0*df/dy)*p
  23 C
  24 C This is computed by a call to F and a difference quotient.
  25 C-----------------------------------------------------------------------
  26 C
  27 C      On entry
  28 C
  29 C          NEQ = problem size, passed to F and PSOL (NEQ(1) = N).
  30 C
  31 C            Y = array containing current dependent variable vector.
  32 C
  33 C         SAVF = array containing current value of f(t,y).
  34 C
  35 C            P = real array of length N.
  36 C
  37 C         WGHT = array of length N containing scale factors.
  38 C                1/WGHT(i) are the diagonal elements of the matrix D.
  39 C
  40 C           WK = work array of length N.
  41 C
  42 C      On return
  43 C
  44 C
  45 C            W = array of length N containing desired
  46 C                matrix-vector product.
  47 C
  48 C In addition, this routine uses the Common variables TN, N, NFE.
  49 C-----------------------------------------------------------------------
  50       INTEGER I
  51       DOUBLE PRECISION FAC, PNRM, RPNRM, DVNORM
  52 C
  53       PNRM = DVNORM (N, P, WGHT)
  54       RPNRM = 1.0D0/PNRM
  55       CALL DCOPY (N, Y, 1, W, 1)
  56       DO 20 I = 1,N
  57  20     Y(I) = W(I) + P(I)*RPNRM
  58       CALL F (NEQ, TN, Y, WK)
  59       NFE = NFE + 1
  60       CALL DCOPY (N, W, 1, Y, 1)
  61       FAC = HL0*PNRM
  62       DO 40 I = 1,N
  63  40     W(I) = P(I) - FAC*(WK(I) - SAVF(I))
  64       RETURN
  65 C----------------------- End of Subroutine DATP ------------------------
  66       END