datv.f

Go to the documentation of this file.

C Work performed under the auspices of the U.S. Department of Energy
C by Lawrence Livermore National Laboratory under contract number 
C W-7405-Eng-48.
C
      SUBROUTINE datv (NEQ, Y, TN, YPRIME, SAVR, V, WGHT, YPTEM, RES,
     *   ires, psol, z, vtem, wp, iwp, cj, eplin, ier, nre, npsl,
     *   rpar,ipar)
C
C***BEGIN PROLOGUE  DATV
C***DATE WRITTEN   890101   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This routine computes the product
C
C   Z = (D-inverse)*(P-inverse)*(dF/dY)*(D*V),
C
C where F(Y) = G(T, Y, CJ*(Y-A)), CJ is a scalar proportional to 1/H,
C and A involves the past history of Y.  The quantity CJ*(Y-A) is
C an approximation to the first derivative of Y and is stored
C in the array YPRIME.  Note that dF/dY = dG/dY + CJ*dG/dYPRIME.
C
C D is a diagonal scaling matrix, and P is the left preconditioning
C matrix.  V is assumed to have L2 norm equal to 1.
C The product is stored in Z and is computed by means of a
C difference quotient, a call to RES, and one call to PSOL.
C
C      On entry
C
C          NEQ = Problem size, passed to RES and PSOL.
C
C            Y = Array containing current dependent variable vector.
C
C       YPRIME = Array containing current first derivative of y.
C
C         SAVR = Array containing current value of G(T,Y,YPRIME).
C
C            V = Real array of length NEQ (can be the same array as Z).
C
C         WGHT = Array of length NEQ containing scale factors.
C                1/WGHT(I) are the diagonal elements of the matrix D.
C
C        YPTEM = Work array of length NEQ.
C
C         VTEM = Work array of length NEQ used to store the
C                unscaled version of V.
C
C         WP = Real work array used by preconditioner PSOL.
C
C         IWP = Integer work array used by preconditioner PSOL.
C
C           CJ = Scalar proportional to current value of 
C                1/(step size H).
C
C
C      On return
C
C            Z = Array of length NEQ containing desired scaled
C                matrix-vector product.
C
C         IRES = Error flag from RES.
C
C          IER = Error flag from PSOL.
C
C         NRE  = The number of calls to RES.
C
C         NPSL = The number of calls to PSOL.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   RES, PSOL
C
C***END PROLOGUE  DATV
C
      INTEGER neq, ires, iwp, ier, nre, npsl, ipar
      DOUBLE PRECISION y, tn, yprime, savr, v, wght, yptem, z, vtem,
     1   wp, cj, rpar
      dimension y(*), yprime(*), savr(*), v(*), wght(*), yptem(*),
     1   z(*), vtem(*), wp(*), iwp(*), rpar(*), ipar(*)
      INTEGER i
      DOUBLE PRECISION eplin
      EXTERNAL  res, psol
C
      ires = 0
C-----------------------------------------------------------------------
C Set VTEM = D * V.
C-----------------------------------------------------------------------
      DO 10 i = 1,neq
 10     vtem(i) = v(i)/wght(i)
      ier = 0
C-----------------------------------------------------------------------
C Store Y in Z and increment Z by VTEM.
C Store YPRIME in YPTEM and increment YPTEM by VTEM*CJ.
C-----------------------------------------------------------------------
      DO 20 i = 1,neq
        yptem(i) = yprime(i) + vtem(i)*cj
 20     z(i) = y(i) + vtem(i)
C-----------------------------------------------------------------------
C Call RES with incremented Y, YPRIME arguments
C stored in Z, YPTEM.  VTEM is overwritten with new residual.
C-----------------------------------------------------------------------
      CONTINUE
      CALL res(tn,z,yptem,cj,vtem,ires,rpar,ipar)
      nre = nre + 1
      IF (ires .LT. 0) RETURN
C-----------------------------------------------------------------------
C Set Z = (dF/dY) * VBAR using difference quotient.
C (VBAR is old value of VTEM before calling RES)
C-----------------------------------------------------------------------
      DO 70 i = 1,neq
 70     z(i) = vtem(i) - savr(i)
C-----------------------------------------------------------------------
C Apply inverse of left preconditioner to Z.
C-----------------------------------------------------------------------
      CALL psol(neq, tn, y, yprime, savr, yptem, cj, wght, wp, iwp,
     1   z, eplin, ier, rpar, ipar)
      npsl = npsl + 1
      IF (ier .NE. 0) RETURN
C-----------------------------------------------------------------------
C Apply D-inverse to Z and return.
C-----------------------------------------------------------------------
      DO 90 i = 1,neq
 90     z(i) = z(i)*wght(i)
      RETURN
C
C------END OF SUBROUTINE DATV-------------------------------------------
      END