subroutine imk_evaluate(this, lev, t, m, flags, step)
use pf_mod_dtype
class(pf_imk_t), intent(inout) :: this
real(pfdp), intent(in ) :: t
integer, intent(in ) :: m
class(pf_level_t), intent(inout) :: lev
integer, optional, intent(in) :: flags, step
integer :: i
real(pfdp) :: dt
! Propagate to get y=exp(Om)
!prop needs e^{Q (omega)} and apply to Y
if (this%debug) &
print*, 'level', lev%index, 'in evaluate ', m, '------------------'
if (this%rk) then
! 't' in f_evals are meaningless since i have a time-independent matrix, A
dt = this%dt
do i = 1, 5
call lev%Q(i)%setval(0.0_pfdp)
end do
call lev%Q(1)%copy(lev%q0, flags=1)
call this%f_eval(lev%Q(1), t, lev%index, this%A(1))
! commutator_p flags=1 hack copies Q(1)%y -> Q(1)%array
! all subsequent RK stages are done on Q(m)%array
call this%commutator_p(this%A(1), lev%Q(1), lev%F(1,1), flags=1)
call lev%Q(2)%axpy(1.0_pfdp, lev%Q(1))
call lev%Q(2)%axpy(0.5_pfdp*dt, lev%F(1,1))
call this%f_eval(lev%Q(2), t, lev%index, this%A(2))
call this%commutator_p(this%A(2), lev%Q(2), lev%F(2,1))
call lev%Q(3)%axpy(1.0_pfdp, lev%Q(1))
call lev%Q(3)%axpy(0.5_pfdp*dt, lev%F(2,1))
call this%f_eval(lev%Q(3), t, lev%index, this%A(3))
call this%commutator_p(this%A(3), lev%Q(3), lev%F(3,1))
call lev%Q(4)%axpy(1.0_pfdp, lev%Q(1))
call lev%Q(4)%axpy(dt, lev%F(3,1))
call this%f_eval(lev%Q(4), t, lev%index, this%A(4))
call this%commutator_p(this%A(4), lev%Q(4), lev%F(4,1))
call lev%Q(5)%axpy(1.0_pfdp, lev%Q(1))
call lev%Q(5)%axpy(dt/6.0_pfdp, lev%F(1,1))
call lev%Q(5)%axpy(dt/3.0_pfdp, lev%F(2,1))
call lev%Q(5)%axpy(dt/3.0_pfdp, lev%F(3,1))
call lev%Q(5)%axpy(dt/6.0_pfdp, lev%F(4,1))
! commutator_p flags=2 hack copies Q(m)%array -> Q(m)%y
! only the Q argument in this case matters
! in order to get solution back to qend
call this%commutator_p(this%A(1), lev%Q(1), lev%F(1,1), flags=2)
call this%commutator_p(this%A(2), lev%Q(2), lev%F(1,1), flags=2)
call this%commutator_p(this%A(3), lev%Q(3), lev%F(1,1), flags=2)
call this%commutator_p(this%A(4), lev%Q(4), lev%F(1,1), flags=2)
call this%commutator_p(this%A(5), lev%Q(5), lev%F(1,1), flags=2)
else
if (this%debug) then
endif
if (m > 1) then
if (this%debug) print*, 'propagating'
call this%propagate(lev%q0, lev%Q(m))
else
if (this%debug) print*, 'copying'
call lev%Q(1)%copy(lev%q0, flags=1)
end if
if (this%debug) print*, 'Q'
if (this%debug) call lev%Q(m)%eprint()
! Compute A(y,t)
call this%f_eval(lev%Q(m), t, lev%index, this%A(m))
if (this%debug) print*, 'A'
if (this%debug) call this%A(m)%eprint()
! Compute the series expansion for dexpinv
if (m > 1) then
call this%dexpinv(this%A(m), lev%Q(m), lev%F(m,1))
else
call lev%F(1,1)%copy(this%A(1))
endif
if (this%debug) print*, 'depxinv'
if (this%debug) call lev%F(m,1)%eprint()
endif
if (this%debug) print*, 'end evaluate ------------'
end subroutine imk_evaluate