dprepj Subroutine

subroutine dprepj(Neq, Y, Yh, Nyh, Ewt, Ftem, Savf, Wm, Iwm, f, jac)

NAME

dprepj(3f) - [M_odepack] Compute and process Newton iteration matrix.

DESCRIPTION

DPREPJ is called by DSTODE to compute and process the matrix P = I - hel(1)J , where J is an approximation to the Jacobian.

Here J is computed by the user-supplied routine JAC if MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5.

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

J is stored in WM and replaced by P. If MITER .ne. 3, P is then subjected to LU decomposition in preparation for later solution of linear systems with P as coefficient matrix. This is done by DGEFA if MITER = 1 or 2, and by DGBFA if MITER = 4 or 5.

In addition to variables described in DSTODE and DLSODE prologues, communication with DPREPJ 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.

WM

real work space for matrices. On output it contains the inverse diagonal matrix if MITER = 3 and the LU decomposition of P if MITER is 1, 2 , 4, or 5. Storage of matrix elements starts at WM(3). WM also contains the following matrix-related data: WM(1) = SQRT(UROUND), used in numerical Jacobian increments. WM(2) = H*EL0, saved for later use if MITER = 3.

IWM

integer work space containing pivot information, starting at IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.

EL0

EL(1) (input).

IERPJ

output error flag, = 0 if no trouble, .gt. 0 if P matrix found to be singular.

JCUR

output flag = 1 to indicate that the Jacobian matrix (or approximation) is now current.

This routine also uses the COMMON variables EL0, H, TN, UROUND, MITER, N, NFE, and NJE.

Arguments

Type IntentOptional Attributes 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) :: Ftem(*)
real(kind=dp) :: Savf(*)
real(kind=dp), intent(inout) :: Wm(*)
integer :: Iwm(*)
real :: f
integer :: jac

Calls

proc~~dprepj~~CallsGraph proc~dprepj dprepj.inc::dprepj dgefa dgefa proc~dprepj->dgefa dvnorm dvnorm proc~dprepj->dvnorm none~wrapup dprepj::wrapup proc~dprepj->none~wrapup dgbfa dgbfa none~wrapup->dgbfa

Variables

Type Visibility Attributes Name Initial
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 :: i1
integer, public :: i2
integer, public :: ier
integer, public :: ii
integer, public :: j
integer, public :: j1
integer, public :: jj
integer, public :: lenp
integer, public :: mba
integer, public :: mband
integer, public :: meb1
integer, public :: meband
integer, public :: ml
integer, public :: ml3
integer, public :: mu
integer, public :: np1
real(kind=dp), public :: r
real(kind=dp), public :: r0
real(kind=dp), public :: srur
real(kind=dp), public :: yi
real(kind=dp), public :: yj
real(kind=dp), public :: yjj

Subroutines

subroutine wrapup()

Arguments

None

Source Code

subroutine dprepj(Neq,Y,Yh,Nyh,Ewt,Ftem,Savf,Wm,Iwm,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)               :: Ftem(*)
real(kind=dp)               :: Savf(*)
real(kind=dp),intent(inout) :: Wm(*)
integer                     :: Iwm(*)
external                    :: f
external                    :: jac
!
real(kind=dp) :: con, di, fac, hl0, r, r0, srur, yi, yj, yjj
integer :: i, i1, i2, ier, ii, j, j1, jj, lenp, mba, mband, meb1, meband, ml, ml3, mu, np1
!
   dls1%nje = dls1%nje + 1
   dls1%ierpj = 0
   dls1%jcur = 1
   hl0 = dls1%h*dls1%el0
   select case (dls1%miter)
   case (2)
   !  If MITER = 2, make N calls to F to approximate J. --------------------
      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 = Wm(1)
      j1 = 2
      do j = 1 , dls1%n
         yj = Y(j)
         r = max(srur*abs(yj),r0/Ewt(j))
         Y(j) = Y(j) + r
         fac = -hl0/r
         call f(Neq,dls1%tn,Y,Ftem)
         do i = 1 , dls1%n
            Wm(i+j1) = (Ftem(i)-Savf(i))*fac
         enddo
         Y(j) = yj
         j1 = j1 + dls1%n
      enddo
      dls1%nfe = dls1%nfe + dls1%n
   case (3)
   !  If MITER = 3, construct a diagonal approximation to J and P. ---------
      Wm(2) = hl0
      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,Wm(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*(Wm(i+2)-Savf(i))
         Wm(i+2) = 1.0D0
         if ( abs(r0)>=dls1%uround/Ewt(i) ) then
            if ( abs(di)==0.0D0 ) then
               dls1%ierpj = 1
               return
            else
               Wm(i+2) = 0.1D0*r0/di
            endif
         endif
      enddo
      return
   case (4)
   !  If MITER = 4, call JAC and multiply by scalar. -----------------------
      ml = Iwm(1)
      mu = Iwm(2)
      ml3 = ml + 3
      mband = ml + mu + 1
      meband = mband + ml
      lenp = meband*dls1%n
      do i = 1 , lenp
         Wm(i+2) = 0.0D0
      enddo
      call jac(Neq,dls1%tn,Y,ml,mu,Wm(ml3),meband)
      con = -hl0
      do i = 1 , lenp
         Wm(i+2) = Wm(i+2)*con
      enddo
      call wrapup()
      return
   case (5)
   !  If MITER = 5, make MBAND calls to F to approximate J. ----------------
      ml = Iwm(1)
      mu = Iwm(2)
      mband = ml + mu + 1
      mba = min(mband,dls1%n)
      meband = mband + ml
      meb1 = meband - 1
      srur = Wm(1)
      fac = dvnorm(dls1%n,Savf,Ewt)
      r0 = 1000.0D0*abs(dls1%h)*dls1%uround*dls1%n*fac
      if ( r0==0.0D0 ) r0 = 1.0D0
      do j = 1 , mba
         do i = j , dls1%n , mband
            yi = Y(i)
            r = max(srur*abs(yi),r0/Ewt(i))
            Y(i) = Y(i) + r
         enddo
         call f(Neq,dls1%tn,Y,Ftem)
         do jj = j , dls1%n , mband
            Y(jj) = Yh(jj,1)
            yjj = Y(jj)
            r = max(srur*abs(yjj),r0/Ewt(jj))
            fac = -hl0/r
            i1 = max(jj-mu,1)
            i2 = min(jj+ml,dls1%n)
            ii = jj*meb1 - ml + 2
            do i = i1 , i2
               Wm(ii+i) = (Ftem(i)-Savf(i))*fac
            enddo
         enddo
      enddo
      dls1%nfe = dls1%nfe + mba
      call wrapup()
      return
   case default
   !  If MITER = 1, call JAC and multiply by scalar. -----------------------
      lenp = dls1%n*dls1%n
      do i = 1 , lenp
         Wm(i+2) = 0.0D0
      enddo
      call jac(Neq,dls1%tn,Y,0,0,Wm(3),dls1%n)
      con = -hl0
      do i = 1 , lenp
         Wm(i+2) = Wm(i+2)*con
      enddo
   endselect
   !  Add identity matrix. -------------------------------------------------
   j = 3
   np1 = dls1%n + 1
   do i = 1 , dls1%n
      Wm(j) = Wm(j) + 1.0D0
      j = j + np1
   enddo
   !  Do LU decomposition on P. --------------------------------------------
   call dgefa(Wm(3),dls1%n,dls1%n,Iwm(21),ier)
   if ( ier/=0 ) dls1%ierpj = 1

contains

subroutine wrapup()
   !  Add identity matrix. -------------------------------------------------
   ii = mband + 2
   do i = 1 , dls1%n
      Wm(ii) = Wm(ii) + 1.0D0
      ii = ii + meband
   enddo
   !  Do LU decomposition of P. --------------------------------------------
   call dgbfa(Wm(3),meband,dls1%n,ml,mu,Iwm(21),ier)
   if ( ier/=0 ) dls1%ierpj = 1
end subroutine wrapup

end subroutine dprepj