dvjust Subroutine

private subroutine dvjust(me, yh, ldyh, iord)

this subroutine adjusts the yh array on reduction of order, and also when the order is increased for the stiff option (meth = 2).

Type Bound

dvode_t

Arguments

Type IntentOptional Attributes Name
class(dvode_t), intent(inout) :: me
real(kind=wp), intent(inout) :: yh(ldyh,*)
integer, intent(in) :: ldyh

leading dimension of yh
integer, intent(in) :: iord

an integer flag used when meth = 2 to indicate an order increase (iord = +1) or an order decrease (iord = -1).

Calls

proc~~dvjust~~CallsGraph proc~dvjust dvode_module::dvode_t%dvjust proc~daxpy dvode_blas_module::daxpy proc~dvjust->proc~daxpy

Called by

proc~~dvjust~~CalledByGraph proc~dvjust dvode_module::dvode_t%dvjust proc~dvstep dvode_module::dvode_t%dvstep proc~dvstep->proc~dvjust proc~dvode dvode_module::dvode_t%dvode proc~dvode->proc~dvstep

Source Code

   subroutine dvjust(me,yh,ldyh,iord)

      class(dvode_t),intent(inout) :: me
      integer,intent(in) :: ldyh !! leading dimension of `yh`
      integer,intent(in) :: iord !! an integer flag used when meth = 2 to indicate an order
                                 !! increase (iord = +1) or an order decrease (iord = -1).
      real(wp),intent(inout) :: yh(ldyh,*)

      real(wp) :: alph0 , alph1 , hsum , prod , t1 , xi , xiold
      integer :: i , iback , j , jp1 , lp1 , nqm1 , nqm2 , nqp1

      real(wp),parameter :: one = 1.0_wp

      if ( (me%dat%nq==2) .and. (iord/=1) ) return
      nqm1 = me%dat%nq - 1
      nqm2 = me%dat%nq - 2
      if ( me%dat%meth==2 ) then
         !-----------------------------------------------------------------------
         ! stiff option...
         ! check to see if the order is being increased or decreased.
         !-----------------------------------------------------------------------
         if ( iord==1 ) then
            ! order increase. ------------------------------------------------------
            do j = 1 , me%dat%lmax
               me%dat%el(j) = zero
            enddo
            me%dat%el(3) = one
            alph0 = -one
            alph1 = one
            prod = one
            xiold = one
            hsum = me%dat%hscal
            if ( me%dat%nq/=1 ) then
               do j = 1 , nqm1
                  ! construct coefficients of x*x*(x+xi(1))*...*(x+xi(j)). ---------------
                  jp1 = j + 1
                  hsum = hsum + me%dat%tau(jp1)
                  xi = hsum/me%dat%hscal
                  prod = prod*xi
                  alph0 = alph0 - one/real(jp1,wp)
                  alph1 = alph1 + one/xi
                  do iback = 1 , jp1
                     i = (j+4) - iback
                     me%dat%el(i) = me%dat%el(i)*xiold + me%dat%el(i-1)
                  enddo
                  xiold = xi
               enddo
            endif
            t1 = (-alph0-alph1)/prod
            ! load column l + 1 in yh array. ---------------------------------------
            lp1 = me%dat%l + 1
            do i = 1 , me%dat%n
               yh(i,lp1) = t1*yh(i,me%dat%lmax)
            enddo
            ! add correction terms to yh array. ------------------------------------
            nqp1 = me%dat%nq + 1
            do j = 3 , nqp1
               call daxpy(me%dat%n,me%dat%el(j),yh(1,lp1),1,yh(1,j),1)
            enddo
         else
            ! order decrease. ------------------------------------------------------
            do j = 1 , me%dat%lmax
               me%dat%el(j) = zero
            enddo
            me%dat%el(3) = one
            hsum = zero
            do j = 1 , nqm2
               ! construct coefficients of x*x*(x+xi(1))*...*(x+xi(j)). ---------------
               hsum = hsum + me%dat%tau(j)
               xi = hsum/me%dat%hscal
               jp1 = j + 1
               do iback = 1 , jp1
                  i = (j+4) - iback
                  me%dat%el(i) = me%dat%el(i)*xi + me%dat%el(i-1)
               enddo
            enddo
            ! subtract correction terms from yh array. -----------------------------
            do j = 3 , me%dat%nq
               do i = 1 , me%dat%n
                  yh(i,j) = yh(i,j) - yh(i,me%dat%l)*me%dat%el(j)
               enddo
            enddo
            return
         endif
      !-----------------------------------------------------------------------
      ! nonstiff option...
      ! check to see if the order is being increased or decreased.
      !-----------------------------------------------------------------------
      else if ( iord==1 ) then
         ! order increase. ------------------------------------------------------
         ! zero out next column in yh array. ------------------------------------
         lp1 = me%dat%l + 1
         do i = 1 , me%dat%n
            yh(i,lp1) = zero
         enddo
         return
      else
         ! order decrease. ------------------------------------------------------
         do j = 1 , me%dat%lmax
            me%dat%el(j) = zero
         enddo
         me%dat%el(2) = one
         hsum = zero
         do j = 1 , nqm2
            ! construct coefficients of x*(x+xi(1))*...*(x+xi(j)). -----------------
            hsum = hsum + me%dat%tau(j)
            xi = hsum/me%dat%hscal
            jp1 = j + 1
            do iback = 1 , jp1
               i = (j+3) - iback
               me%dat%el(i) = me%dat%el(i)*xi + me%dat%el(i-1)
            enddo
         enddo
         ! construct coefficients of integrated polynomial. ---------------------
         do j = 2 , nqm1
            me%dat%el(j+1) = real(me%dat%nq,wp)*me%dat%el(j)/real(j,wp)
         enddo
         ! subtract correction terms from yh array. -----------------------------
         do j = 3 , me%dat%nq
            do i = 1 , me%dat%n
               yh(i,j) = yh(i,j) - yh(i,me%dat%l)*me%dat%el(j)
            enddo
         enddo
         return
      endif

   end subroutine dvjust