pf_apply_mat_backward Subroutine

  subroutine pf_apply_mat_backward(dst, a, mat, src, zero, flags)
    !! Apply a matrix (tmat or rmat) to src and add to dst.
    class(pf_encap_t), intent(inout) :: dst(:)       !!  destination vector
    real(pfdp),        intent(in)    :: a            !!  scalar
    real(pfdp),        intent(in)    :: mat(:, :)    !!  matrix
    class(pf_encap_t), intent(in)    :: src(:)       !!  src vector
    logical,           intent(in), optional :: zero   !! If false, don't zero out the the dst variable before computing 
    integer,           intent(in), optional :: flags  !! Used for choosing which variable to operate on 

    
    !!  Local variables
    logical :: lzero   !!  local version of input parameter zero
    integer :: which   !!  local version of flags
    integer :: n, m    !!  size of mat   
    integer :: i, j    !!  loop variables

    lzero = .true.; if (present(zero)) lzero = zero    
    which = 2;      if(present(flags)) which = flags
    
    if( which /= 2 ) &
      stop "pf_apply_mat_backward can only be used for restricting the backward integrals with which==2"

    n = size(mat, dim=1)
    m = size(mat, dim=2)
        
    do i = 1, n
      if (lzero) call dst(n+1-i)%setval(0.0_pfdp, 2)
      do j = 1, m
        if (a*mat(i, j) /= 0.0_pfdp)  call dst(n+1-i)%axpy(a * mat(i, j), src(m+1-j), 2)
      end do
    end do
  end subroutine pf_apply_mat_backward