dvjac Subroutine

private subroutine dvjac(me, y, yh, ldyh, ewt, ftem, savf, wm, iwm, ierpj)

dvjac is called by dvnlsd to compute and process the matrix p = i - h*rl1*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. if jsv = -1, j is computed from scratch in all cases. if jsv = 1 and miter = 1, 2, 4, or 5, and if the saved value of j is considered acceptable, then p is constructed from the saved j. j is stored in wm and replaced by p. if miter /= 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.

Type Bound

dvode_t

Arguments

Type IntentOptional Attributes Name
class(dvode_t), intent(inout) :: me
real(kind=wp), intent(inout) :: y(*)

vector containing predicted values on entry.
real(kind=wp), intent(in) :: yh(ldyh,*)

the nordsieck array, an ldyh by lmax array.
integer, intent(in) :: ldyh

a constant >= n, the first dimension of yh.
real(kind=wp), intent(in) :: ewt(*)

an error weight vector of length n.
real(kind=wp), intent(in) :: ftem(*)
real(kind=wp), intent(in) :: savf(*)

array containing f evaluated at predicted y.
real(kind=wp), intent(inout) :: wm(*)

real work space for matrices. in the 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). storage of the saved jacobian starts at wm(locjs). wm also contains the following matrix-related data: wm(1) = sqrt(uround), used in numerical jacobian step. wm(2) = h*rl1, saved for later use if miter = 3.
integer, intent(inout) :: iwm(*)

integer work space containing pivot information, starting at iwm(31), 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.
integer, intent(out) :: ierpj

output error flag, = 0 if no trouble, 1 if the p matrix is found to be singular.

Calls

proc~~dvjac~~CallsGraph proc~dvjac dvode_module::dvode_t%dvjac proc~dacopy dvode_module::dacopy proc~dvjac->proc~dacopy proc~dcopy dvode_blas_module::dcopy proc~dvjac->proc~dcopy proc~dgbfa dvode_linpack_module::dgbfa proc~dvjac->proc~dgbfa proc~dgefa dvode_linpack_module::dgefa proc~dvjac->proc~dgefa proc~dscal dvode_blas_module::dscal proc~dvjac->proc~dscal proc~dacopy->proc~dcopy proc~dgbfa->proc~dscal proc~daxpy dvode_blas_module::daxpy proc~dgbfa->proc~daxpy proc~idamax dvode_blas_module::idamax proc~dgbfa->proc~idamax proc~dgefa->proc~dscal proc~dgefa->proc~daxpy proc~dgefa->proc~idamax

Called by

proc~~dvjac~~CalledByGraph proc~dvjac dvode_module::dvode_t%dvjac proc~dvnlsd dvode_module::dvode_t%dvnlsd proc~dvnlsd->proc~dvjac proc~dvstep dvode_module::dvode_t%dvstep proc~dvstep->proc~dvnlsd proc~dvode dvode_module::dvode_t%dvode proc~dvode->proc~dvstep

Source Code

   subroutine dvjac(me,y,yh,ldyh,ewt,ftem,savf,wm,iwm,ierpj)

      class(dvode_t),intent(inout) :: me
      real(wp),intent(inout) :: y(*) !! vector containing predicted values on entry.
      integer,intent(in) :: ldyh !! a constant >= n, the first dimension of `yh`.
      real(wp),intent(in) :: yh(ldyh,*) !! the nordsieck array, an `ldyh` by `lmax` array.
      real(wp),intent(in) :: ewt(*) !! an error weight vector of length `n`.
      real(wp),intent(in) :: ftem(*)
      real(wp),intent(in) :: savf(*) !! array containing f evaluated at predicted y.
      real(wp),intent(inout) :: wm(*) !! real work space for matrices.  in the 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).
                                      !! storage of the saved jacobian starts at wm(locjs).
                                      !! wm also contains the following matrix-related data:
                                      !! wm(1) = sqrt(uround), used in numerical jacobian step.
                                      !! wm(2) = h*rl1, saved for later use if miter = 3.
      integer,intent(inout) :: iwm(*) !! integer work space containing pivot information,
                                      !! starting at iwm(31), 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.
      integer,intent(out) :: ierpj !! output error flag,  = 0 if no trouble, 1 if the p
                                   !! matrix is found to be singular.

      real(wp) :: con , di , fac , hrl1 , r , r0 , &
                  srur , yi , yj , yjj
      integer :: i , i1 , i2 , ier , ii , j , j1 , jj , jok , lenp , mba , &
                 mband , meb1 , meband , ml , ml3 , mu , np1

      real(wp),parameter :: one = 1.0_wp
      real(wp),parameter :: thou = 1000.0_wp
      real(wp),parameter :: pt1 = 0.1_wp

      ierpj = 0
      hrl1 = me%dat%h*me%dat%rl1 ! rl1 = 1/el(2)
      ! see whether j should be evaluated (jok = -1) or not (jok = 1). -------
      jok = me%dat%jsv
      if ( me%dat%jsv==1 ) then
         if ( me%dat%nst==0 .or. me%dat%nst>me%dat%nslj+me%dat%msbj ) jok = -1
         if ( me%dat%icf==1 .and. me%dat%drc<me%dat%ccmxj ) jok = -1
         if ( me%dat%icf==2 ) jok = -1
      endif
      ! end of setting jok. --------------------------------------------------

      if ( jok==-1 .and. me%dat%miter==1 ) then
         ! if jok = -1 and miter = 1, call me%jac to evaluate jacobian. ------------
         me%dat%nje = me%dat%nje + 1
         me%dat%nslj = me%dat%nst
         me%dat%jcur = 1
         lenp = me%dat%n*me%dat%n
         do i = 1 , lenp
            wm(i+2) = zero
         enddo
         call me%jac(me%dat%n,me%dat%tn,y(1:me%dat%n),0,0,wm(3),me%dat%n)
         if ( me%dat%jsv==1 ) call dcopy(lenp,wm(3),1,wm(me%dat%locjs),1)
      endif

      if ( jok==-1 .and. me%dat%miter==2 ) then
         ! if miter = 2, make n calls to f to approximate the jacobian. ---------
         me%dat%nje = me%dat%nje + 1
         me%dat%nslj = me%dat%nst
         me%dat%jcur = 1
         fac = me%dvnorm(me%dat%n,savf(1:me%dat%n),ewt(1:me%dat%n))
         r0 = thou*abs(me%dat%h)*me%dat%uround*real(me%dat%n,wp)*fac
         if ( r0==zero ) r0 = one
         srur = wm(1)
         j1 = 2
         do j = 1 , me%dat%n
            yj = y(j)
            r = max(srur*abs(yj),r0/ewt(j))
            y(j) = y(j) + r
            fac = one/r
            call me%f(me%dat%n,me%dat%tn,y(1:me%dat%n),ftem(1:me%dat%n))
            do i = 1 , me%dat%n
               wm(i+j1) = (ftem(i)-savf(i))*fac
            enddo
            y(j) = yj
            j1 = j1 + me%dat%n
         enddo
         me%dat%nfe = me%dat%nfe + me%dat%n
         lenp = me%dat%n*me%dat%n
         if ( me%dat%jsv==1 ) call dcopy(lenp,wm(3),1,wm(me%dat%locjs),1)
      endif

      if ( jok==1 .and. (me%dat%miter==1 .or. me%dat%miter==2) ) then
         me%dat%jcur = 0
         lenp = me%dat%n*me%dat%n
         call dcopy(lenp,wm(me%dat%locjs),1,wm(3),1)
      endif

      if ( me%dat%miter==1 .or. me%dat%miter==2 ) then
         ! multiply jacobian by scalar, add identity, and do lu decomposition. --
         con = -hrl1
         call dscal(lenp,con,wm(3),1)
         j = 3
         np1 = me%dat%n + 1
         do i = 1 , me%dat%n
            wm(j) = wm(j) + one
            j = j + np1
         enddo
         me%dat%nlu = me%dat%nlu + 1
         call dgefa(wm(3),me%dat%n,me%dat%n,iwm(31),ier)
         if ( ier/=0 ) ierpj = 1
         return
      endif
      ! end of code block for miter = 1 or 2. --------------------------------

      if ( me%dat%miter==3 ) then
         ! if miter = 3, construct a diagonal approximation to j and p. ---------
         me%dat%nje = me%dat%nje + 1
         me%dat%jcur = 1
         wm(2) = hrl1
         r = me%dat%rl1*pt1
         do i = 1 , me%dat%n
            y(i) = y(i) + r*(me%dat%h*savf(i)-yh(i,2))
         enddo
         call me%f(me%dat%n,me%dat%tn,y(1:me%dat%n),wm(3))
         me%dat%nfe = me%dat%nfe + 1
         do i = 1 , me%dat%n
            r0 = me%dat%h*savf(i) - yh(i,2)
            di = pt1*r0 - me%dat%h*(wm(i+2)-savf(i))
            wm(i+2) = one
            if ( abs(r0)>=me%dat%uround/ewt(i) ) then
               if ( abs(di)==zero ) then
                  ierpj = 1
                  return
               end if
               wm(i+2) = pt1*r0/di
            endif
         enddo
         return
      endif
      ! end of code block for miter = 3. -------------------------------------

      ! set constants for miter = 4 or 5. ------------------------------------
      ml = iwm(1)
      mu = iwm(2)
      ml3 = ml + 3
      mband = ml + mu + 1
      meband = mband + ml
      lenp = meband*me%dat%n

      if ( jok==-1 .and. me%dat%miter==4 ) then
         ! if jok = -1 and miter = 4, call me%jac to evaluate jacobian. ------------
         me%dat%nje = me%dat%nje + 1
         me%dat%nslj = me%dat%nst
         me%dat%jcur = 1
         do i = 1 , lenp
            wm(i+2) = zero
         enddo
         call me%jac(me%dat%n,me%dat%tn,y(1:me%dat%n),ml,mu,wm(ml3),meband)
         if ( me%dat%jsv==1 ) call dacopy(mband,me%dat%n,wm(ml3),meband,wm(me%dat%locjs),mband)
      endif

      if ( jok==-1 .and. me%dat%miter==5 ) then
         ! if miter = 5, make ml+mu+1 calls to f to approximate the jacobian. ---
         me%dat%nje = me%dat%nje + 1
         me%dat%nslj = me%dat%nst
         me%dat%jcur = 1
         mba = min(mband,me%dat%n)
         meb1 = meband - 1
         srur = wm(1)
         fac = me%dvnorm(me%dat%n,savf(1:me%dat%n),ewt(1:me%dat%n))
         r0 = thou*abs(me%dat%h)*me%dat%uround*real(me%dat%n,wp)*fac
         if ( r0==zero ) r0 = one
         do j = 1 , mba
            do i = j , me%dat%n , mband
               yi = y(i)
               r = max(srur*abs(yi),r0/ewt(i))
               y(i) = y(i) + r
            enddo
            call me%f(me%dat%n,me%dat%tn,y(1:me%dat%n),ftem(1:me%dat%n))
            do jj = j , me%dat%n , mband
               y(jj) = yh(jj,1)
               yjj = y(jj)
               r = max(srur*abs(yjj),r0/ewt(jj))
               fac = one/r
               i1 = max(jj-mu,1)
               i2 = min(jj+ml,me%dat%n)
               ii = jj*meb1 - ml + 2
               do i = i1 , i2
                  wm(ii+i) = (ftem(i)-savf(i))*fac
               enddo
            enddo
         enddo
         me%dat%nfe = me%dat%nfe + mba
         if ( me%dat%jsv==1 ) call dacopy(mband,me%dat%n,wm(ml3),meband,wm(me%dat%locjs),mband)
      endif

      if ( jok==1 ) then
         me%dat%jcur = 0
         call dacopy(mband,me%dat%n,wm(me%dat%locjs),mband,wm(ml3),meband)
      endif

      ! multiply jacobian by scalar, add identity, and do lu decomposition.
      con = -hrl1
      call dscal(lenp,con,wm(3),1)
      ii = mband + 2
      do i = 1 , me%dat%n
         wm(ii) = wm(ii) + one
         ii = ii + meband
      enddo
      me%dat%nlu = me%dat%nlu + 1
      call dgbfa(wm(3),meband,me%dat%n,ml,mu,iwm(31),ier)
      if ( ier/=0 ) ierpj = 1
      ! end of code block for miter = 4 or 5. --------------------------------

   end subroutine dvjac