dprjs – odepack

subroutine dprjs(Neq, Y, Yh, Nyh, Ewt, Ftem, Savf, Wk, Iwk, f, jac)

DPRJS is called to compute and process the matrix P = I - HEL(1)J, where J is an approximation to the Jacobian. J is computed by columns, either by the user-supplied routine JAC if MITER = 1, or by finite differencing if MITER = 2.

Alternatively, if MITER = 3, a diagonal approximation to J is used.

if MITER = 1 or 2, and if the existing value of the Jacobian (as contained in P) is considered acceptable, then a new value of P is reconstructed from the old value.

In any case, when MITER is 1 or 2, the P matrix is subjected to LU decomposition in CDRV.

P and its LU decomposition are stored (separately) in WK.

In addition to variables described previously, communication with DPRJS uses the following:

Y: array containing predicted values on entry.
FTEM: work array of length N (ACOR in DSTODE).
SAVF: array containing f evaluated at predicted y.
WK: real work space for matrices. On output it contains the inverse diagonal matrix if MITER = 3, and P and its sparse LU decomposition if MITER is 1 or 2. Storage of matrix elements starts at WK(3). WK also contains the following matrix-related data: WK(1) = SQRT(UROUND), used in numerical Jacobian increments. WK(2) = H*EL0, saved for later use if MITER = 3.
IWK: integer work space for matrix-related data, assumed to be equivalenced to WK. In addition, WK(IPRSP) and IWK(IPISP) are assumed to have identical locations.
EL0: EL(1) (input).
IERPJ: output error flag (in Common). = 0 if no error. = 1 if zero pivot found in CDRV. = 2 if a singular matrix arose with MITER = 3. = -1 if insufficient storage for CDRV (should not occur here). = -2 if other error found in CDRV (should not occur here).
JCUR: output flag showing status of (approximate) Jacobian matrix: = 1 to indicate that the Jacobian is now current, or = 0 to indicate that a saved value was used.

This routine also uses other variables in Common.

Arguments

Type	Intent		Name
integer		::	Neq(*)
real(kind=dp),	intent(inout)	::	Y(*)
real(kind=dp),	intent(in)	::	Yh(Nyh,*)
integer,	intent(in)	::	Nyh
real(kind=dp)		::	Ewt(*)
real(kind=dp),	intent(inout)	::	Ftem(*)
real(kind=dp)		::	Savf(*)
real(kind=dp),	intent(inout)	::	Wk(*)
integer		::	Iwk(*)
real		::	f
integer		::	jac

Calls

Help

Variables

Type	Visibility		Name
real(kind=dp),	public	::	con
real(kind=dp),	public	::	di
real(kind=dp),	public	::	fac
real(kind=dp),	public	::	hl0
integer,	public	::	i
integer,	public	::	imul
integer,	public	::	j
integer,	public	::	jj
integer,	public	::	jmax
integer,	public	::	jmin
integer,	public	::	jok
integer,	public	::	k
integer,	public	::	kmax
integer,	public	::	kmin
integer,	public	::	ng
real(kind=dp),	public	::	pij
real(kind=dp),	public	::	r
real(kind=dp),	public	::	r0
real(kind=dp),	public	::	rcon
real(kind=dp),	public	::	rcont
real(kind=dp),	public	::	srur

Subroutines

subroutine wrapup()

Arguments

None

Source Code

subroutine dprjs(Neq,Y,Yh,Nyh,Ewt,Ftem,Savf,Wk,Iwk,f,jac)

integer                       :: Neq(*)
real(kind=dp),intent(inout)   :: Y(*)
integer,intent(in)            :: Nyh
real(kind=dp),intent(in)      :: Yh(Nyh,*)
real(kind=dp)                 :: Ewt(*)
real(kind=dp),intent(inout)   :: Ftem(*)
real(kind=dp)                 :: Savf(*)
real(kind=dp),intent(inout)   :: Wk(*)
integer                       :: Iwk(*)
external                      :: f
external                      :: jac

real(kind=dp) :: con, di, fac, hl0, pij, r, r0, rcon, rcont, srur
integer :: i, imul, j, jj, jmax, jmin, jok, k, kmax, kmin, ng

   hl0 = dls1%h*dls1%el0
   con = -hl0
   if ( dls1%miter==3 ) then
   !
   !  If MITER = 3, construct a diagonal approximation to J and P. ---------
      dls1%jcur = 1
      dls1%nje = dls1%nje + 1
      Wk(2) = hl0
      dls1%ierpj = 0
      r = dls1%el0*0.1D0
      do i = 1, dls1%n
         Y(i) = Y(i) + r*(dls1%h*Savf(i)-Yh(i,2))
      enddo
      call f(Neq,dls1%tn,Y,Wk(3))
      dls1%nfe = dls1%nfe + 1
      do i = 1, dls1%n
         r0 = dls1%h*Savf(i) - Yh(i,2)
         di = 0.1D0*r0 - dls1%h*(Wk(i+2)-Savf(i))
         Wk(i+2) = 1.0D0
         if ( abs(r0)>=dls1%uround/Ewt(i) ) then
            if ( abs(di)==0.0D0 ) then
               dls1%ierpj = 2
               return
            else
               Wk(i+2) = 0.1D0*r0/di
            endif
         endif
      enddo
      return
   else
   !  See whether J should be reevaluated (JOK = 0) or not (JOK = 1). ------
      jok = 1
      if ( dls1%nst==0 .or. dls1%nst>=dlss%nslj+dlss%msbj ) jok = 0
      if ( dls1%icf==1 .and. abs(dls1%rc-1.0D0)<dlss%ccmxj ) jok = 0
      if ( dls1%icf==2 ) jok = 0
      if ( jok==1 ) then
   !
   !  If JOK = 1, reconstruct new P from old P. ----------------------------
         dls1%jcur = 0
         rcon = con/dlss%con0
         rcont = abs(con)/dlss%conmin
         if ( rcont<=dlss%rbig .or. dlss%iplost/=1 ) then
            kmin = Iwk(dlss%ipian)
            do j = 1, dls1%n
               kmax = Iwk(dlss%ipian+j) - 1
               do k = kmin, kmax
                  i = Iwk(dlss%ibjan+k)
                  pij = Wk(dlss%iba+k)
                  if ( i==j ) then
                     pij = pij - 1.0D0
                     if ( abs(pij)<dlss%psmall ) then
                        dlss%iplost = 1
                        dlss%conmin = min(abs(dlss%con0),dlss%conmin)
                     endif
                  endif
                  pij = pij*rcon
                  if ( i==j ) pij = pij + 1.0D0
                  Wk(dlss%iba+k) = pij
               enddo
               kmin = kmax + 1
            enddo
            call wrapup()
            return
         endif
      endif
   !
   !  MITER = 1 or 2, and the Jacobian is to be reevaluated. ---------------
      dls1%jcur = 1
      dls1%nje = dls1%nje + 1
      dlss%nslj = dls1%nst
      dlss%iplost = 0
      dlss%conmin = abs(con)
      if ( dls1%miter==2 ) then
   !
   !  If MITER = 2, make NGP calls to F to approximate J and P. ------------
         fac = dvnorm(dls1%n,Savf,Ewt)
         r0 = 1000.0D0*abs(dls1%h)*dls1%uround*dls1%n*fac
         if ( r0==0.0D0 ) r0 = 1.0D0
         srur = Wk(1)
         jmin = Iwk(dlss%ipigp)
         do ng = 1, dlss%ngp
            jmax = Iwk(dlss%ipigp+ng) - 1
            do j = jmin, jmax
               jj = Iwk(dlss%ibjgp+j)
               r = max(srur*abs(Y(jj)),r0/Ewt(jj))
               Y(jj) = Y(jj) + r
            enddo
            call f(Neq,dls1%tn,Y,Ftem)
            do j = jmin, jmax
               jj = Iwk(dlss%ibjgp+j)
               Y(jj) = Yh(jj,1)
               r = max(srur*abs(Y(jj)),r0/Ewt(jj))
               fac = -hl0/r
               kmin = Iwk(dlss%ibian+jj)
               kmax = Iwk(dlss%ibian+jj+1) - 1
               do k = kmin, kmax
                  i = Iwk(dlss%ibjan+k)
                  Wk(dlss%iba+k) = (Ftem(i)-Savf(i))*fac
                  if ( i==jj ) Wk(dlss%iba+k) = Wk(dlss%iba+k) + 1.0D0
               enddo
            enddo
            jmin = jmax + 1
         enddo
         dls1%nfe = dls1%nfe + dlss%ngp
      else
   !
   !  If MITER = 1, call JAC, multiply by scalar, and add identity. --------
         kmin = Iwk(dlss%ipian)
         do j = 1, dls1%n
            kmax = Iwk(dlss%ipian+j) - 1
            do i = 1, dls1%n
               Ftem(i) = 0.0D0
            enddo
            call jac(Neq,dls1%tn,Y,j,Iwk(dlss%ipian),Iwk(dlss%ipjan),Ftem)
            do k = kmin, kmax
               i = Iwk(dlss%ibjan+k)
               Wk(dlss%iba+k) = Ftem(i)*con
               if ( i==j ) Wk(dlss%iba+k) = Wk(dlss%iba+k) + 1.0D0
            enddo
            kmin = kmax + 1
         enddo
      endif
   endif

   call wrapup()

contains

subroutine wrapup()
   !
   !  Do numerical factorization of P matrix. ------------------------------
   dlss%nlu = dlss%nlu + 1
   dlss%con0 = con
   dls1%ierpj = 0

   do i = 1, dls1%n
      Ftem(i) = 0.0D0
   enddo

   call cdrv(dls1%n,Iwk(dlss%ipr),Iwk(dlss%ipc),Iwk(dlss%ipic), &
    & Iwk(dlss%ipian),Iwk(dlss%ipjan),Wk(dlss%ipa),Ftem,Ftem,dlss%nsp, &
    & Iwk(dlss%ipisp),Wk(dlss%iprsp),dlss%iesp,2,dlss%iys)

   if ( dlss%iys==0 ) return
   imul = (dlss%iys-1)/dls1%n
   dls1%ierpj = -2
   if ( imul==8 ) dls1%ierpj = 1
   if ( imul==10 ) dls1%ierpj = -1

end subroutine wrapup

end subroutine dprjs

dprjs Subroutine

Contents

Variables

Subroutines

Source Code

subroutine dprjs(Neq, Y, Yh, Nyh, Ewt, Ftem, Savf, Wk, Iwk, f, jac)

Arguments

Calls

Variables

Subroutines

subroutine wrapup()

Arguments

Source Code