datv Subroutine

subroutine datv(Neq, Y, Savf, V, Wght, Ftem, f, psol, Z, Vtem, Wp, Iwp, Hl0, Jpre, Ier, Npsl)

This routine computes the product

  (D-inverse)*(P1-inverse)*(I - hl0*df/dy)*(P2-inverse)*(D*v)

where D is a diagonal scaling matrix, and P1 and P2 are the left and right preconditioning matrices, respectively.

v is assumed to have WRMS norm equal to 1. The product is stored in z. This is computed by a difference quotient, a call to F, and two calls to PSOL.

On entry

NEQ

problem size, passed to F and PSOL (NEQ(1) = N).

Y

array containing current dependent variable vector.

SAVF

array containing current value of f(t,y).

V

real array of length N (can be the same array as Z).

WGHT

array of length N containing scale factors. 1/WGHT(i) are the diagonal elements of the matrix D.

FTEM

work array of length N.

VTEM

work array of length N used to store the unscaled version of V.

WP

real work array used by preconditioner PSOL.

IWP

integer work array used by preconditioner PSOL.

HL0

current value of (step size h) * (coefficient l0).

JPRE

preconditioner type flag.

On return

Z

array of length N containing desired scaled matrix-vector product.

IER

error flag from PSOL.

NPSL

the number of calls to PSOL.

In addition, this routine uses the Common variables TN, N, NFE.

Arguments

Type IntentOptional Attributes Name
integer :: Neq(*)
real(kind=dp) :: Y(*)
real(kind=dp) :: Savf(*)
real(kind=dp), intent(in) :: V(*)
real(kind=dp), intent(in) :: Wght(*)
real(kind=dp) :: Ftem(*)
real :: f
real :: psol
real(kind=dp), intent(inout) :: Z(*)
real(kind=dp), intent(inout) :: Vtem(*)
real(kind=dp), intent(inout) :: Wp(*)
integer :: Iwp(*)
real(kind=dp) :: Hl0
integer, intent(inout) :: Jpre
integer, intent(inout) :: Ier
integer, intent(inout) :: Npsl

Calls

proc~~datv~~CallsGraph proc~datv datv.inc::datv dnrm2 dnrm2 proc~datv->dnrm2

Variables

Type Visibility Attributes Name Initial
real(kind=dp), public :: fac
integer, public :: i
real(kind=dp), public :: rnorm
real(kind=dp), public :: tempn

Source Code

subroutine datv(Neq,Y,Savf,V,Wght,Ftem,f,psol,Z,Vtem,Wp,Iwp,Hl0,Jpre,Ier,Npsl)

integer                     :: Neq(*)
real(kind=dp)               :: Y(*)
real(kind=dp)               :: Savf(*)
real(kind=dp),intent(in)    :: V(*)
real(kind=dp),intent(in)    :: Wght(*)
real(kind=dp)               :: Ftem(*)
external                    :: f
external                    :: psol
real(kind=dp),intent(inout) :: Z(*)
real(kind=dp),intent(inout) :: Vtem(*)
real(kind=dp),intent(inout) :: Wp(*)
integer                     :: Iwp(*)
real(kind=dp)               :: Hl0
integer, intent(inout)      :: Jpre
integer, intent(inout)      :: Ier
integer, intent(inout)      :: Npsl

real(kind=dp) :: fac, rnorm, tempn
integer :: i

   !  Set VTEM = D * V.
   do i = 1, dls1%n
      Vtem(i) = V(i)/Wght(i)
   enddo

   Ier = 0

   if ( Jpre>=2 ) then

      !  JPRE = 2 or 3.  Apply inverse of right preconditioner to VTEM.
      call psol(Neq,dls1%tn,Y,Savf,Ftem,Hl0,Wp,Iwp,Vtem,2,Ier)
      Npsl = Npsl + 1

      if ( Ier/=0 ) return

      !  Calculate L-2 norm of (D-inverse) * VTEM.
      do i = 1, dls1%n
         Z(i) = Vtem(i)*Wght(i)
      enddo

      tempn = dnrm2(dls1%n,Z,1)
      rnorm = 1.0D0/tempn

      !  Save Y in Z and increment Y by VTEM/norm.
      !X!call dcopy(dls1%n,Y,1,Z,1)
      Z(1:dls1%n)=Y(1:dls1%n)!X!

      do i = 1, dls1%n
         Y(i) = Z(i) + Vtem(i)*rnorm
      enddo

      fac = Hl0*tempn
   else

      !  JPRE = 0 or 1.  Save Y in Z and increment Y by VTEM.
      !X!call dcopy(dls1%n,Y,1,Z,1)
      Z(1:dls1%n)=Y(1:dls1%n)!X!

      do i = 1, dls1%n
         Y(i) = Z(i) + Vtem(i)
      enddo

      fac = Hl0
   endif

   !  For all JPRE, call F with incremented Y argument, and restore Y.
   call f(Neq,dls1%tn,Y,Ftem)
   dls1%nfe = dls1%nfe + 1
   !X!call dcopy(dls1%n,Z,1,Y,1)
   Y(1:dls1%n)=Z(1:dls1%n)!X!

   !  Set Z = (identity - hl0*Jacobian) * VTEM, using difference quotient.
   do i = 1, dls1%n
      Z(i) = Ftem(i) - Savf(i)
   enddo
   do i = 1, dls1%n
      Z(i) = Vtem(i) - fac*Z(i)
   enddo

   !  Apply inverse of left preconditioner to Z, if nontrivial.
   if ( Jpre/=0 .and. Jpre/=2 ) then
      call psol(Neq,dls1%tn,Y,Savf,Ftem,Hl0,Wp,Iwp,Z,1,Ier)
      Npsl = Npsl + 1
      if ( Ier/=0 ) return
   endif

   !  Apply D-inverse to Z and return.
   do i = 1, dls1%n
      Z(i) = Z(i)*Wght(i)
   enddo

end subroutine datv