anchor.diff:          real(sp), intent(inout), dimension(n * n) :: y
           real(sp), intent(in), dimension(2 * m - 1, n) :: blocks
   
  -        integer :: i
  -
  -        i = 1
  -        call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:4, :), ms(2), n)
  -        call rgbmv(x((i + 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(5:5, :), ms(3), n)
  -
  -        do i = 2, n - 1
  -            call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
  -            call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:4, :), ms(2), n)
  -            call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(5:5, :), ms(3), n)
  -        end do
  -
  -        i = n
  -        call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(5:5, :), ms(1), n)
  -        call rgbmv(x((i - 1) * n + 1:(i - 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:4, :), ms(2), n)
  -    end subroutine rbbmv_o3
  -
  -    subroutine rbbmv_o5(x, y, sign, blocks, ms, m, n)
  -        implicit none
  -
  -        integer, intent(in) :: m, n
  -        integer, intent(in), dimension(m) :: ms
  -        real(sp), intent(in) :: sign
  -        real(sp), intent(in), dimension(n * n) :: x
  -        real(sp), intent(inout), dimension(n * n) :: y
  -        real(sp), intent(in), dimension(2 * m - 1, n) :: blocks
  -
  -        integer :: i
  -
  -        i = 1
  -        call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:5, :), ms(3), n)
  -        call rgbmv(x((i + 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(6:6, :), ms(4), n)
  -        call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(7:7, :), ms(5), n)
  -
  -        i = 2
  -        call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
  -        call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:5, :), ms(3), n)
  -        call rgbmv(x((i + 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(6:6, :), ms(4), n)
  -        call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(7:7, :), ms(5), n)
  -
  -        do i = 3, n - 2
  -            call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
  -            call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
  -            call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:5, :), ms(3), n)
  -            call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(6:6, :), ms(4), n)
  -            call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(7:7, :), ms(5), n)
  -        end do
  -
  -        i = n - 1
  -        call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
  -        call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
  -        call rgbmv(x((i - 1) * n + 1:(i - 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:5, :), ms(3), n)
  -        call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(6:6, :), ms(4), n)
  -
  -        i = n
  -        call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
  -        call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
  -        call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:5, :), ms(3), n)
  -    end subroutine rbbmv_o5
  -
  -    subroutine rbbmv(x, y, sign, blocks, ms, m, n)
  -        implicit none
  -
  -        integer, intent(in) :: m, n
  -        integer, intent(in), dimension(m) :: ms
  -        real(sp), intent(in) :: sign
  -        real(sp), intent(in), dimension(n * n) :: x
  -        real(sp), intent(inout), dimension(n * n) :: y
  -        real(sp), intent(in), dimension(2 * m - 1, n) :: blocks
  -
           if (m == 3) then
  -            call rbbmv_o3(x, y, sign, blocks, ms, m, n)
  +            call rbbmv_o3(x, y, sign, blocks, ms, n)
           else if (m == 5) then
  -            call rbbmv_o5(x, y, sign, blocks, ms, m, n)
  +            call rbbmv_o5(x, y, sign, blocks, ms, n)
           else if (m == 7) then
  -            ! TODO: implement this
  +            call rbbmv_o7(x, y, sign, blocks, ms, n)
           end if
       end subroutine rbbmv
   
  
positive.code:     public :: divide_derivative_on_radius
      public :: make_laplacian_o3, make_laplacian_o5, make_laplacian_o7, make_laplacian
      public :: make_laplacian_2d_o3, make_laplacian_2d_o5, make_laplacian_2d_o7, make_laplacian_2d
      public :: rbbmv_o3, rbbmv_o5, rbbmv, rgbmv
      public :: revervoir, revervoir_2d
      public :: hamiltonian, hamiltonian_2d
      public :: runge_kutta, runge_kutta_2d
  
negative.code:     !   \brief Perform one of the matrix-vector operations   y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y in
      !   case of block band matrix.
      !   TODO: support high-order approximations.
      subroutine rbbmv_o3(x, y, sign, blocks, ms, n)
          implicit none
  
          integer, parameter :: m = 3
          integer, intent(in) :: n
          integer, intent(in), dimension(m) :: ms
          real(sp), intent(in) :: sign
          real(sp), intent(in), dimension(n * n) :: x
          real(sp), intent(inout), dimension(n * n) :: y
          real(sp), intent(in), dimension(2 * m - 1, n) :: blocks
  
          integer :: i
  
          i = 1
          call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:4, :), ms(2), n)
          call rgbmv(x((i + 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(5:5, :), ms(3), n)
  
          do i = 2, n - 1
              call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
              call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:4, :), ms(2), n)
              call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(5:5, :), ms(3), n)
          end do
  
          i = n
          call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(5:5, :), ms(1), n)
          call rgbmv(x((i - 1) * n + 1:(i - 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:4, :), ms(2), n)
      end subroutine rbbmv_o3
  
      subroutine rbbmv_o5(x, y, sign, blocks, ms, n)
          implicit none
  
          integer, parameter :: m = 5
          integer, intent(in) :: n
          integer, intent(in), dimension(m) :: ms
          real(sp), intent(in) :: sign
          real(sp), intent(in), dimension(n * n) :: x
          real(sp), intent(inout), dimension(n * n) :: y
          real(sp), intent(in), dimension(2 * m - 1, n) :: blocks
  
          integer :: i
  
          i = 1
          call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:7, :), ms(3), n)
          call rgbmv(x((i + 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(8:8, :), ms(4), n)
          call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(9:9, :), ms(5), n)
  
          i = 2
          call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
          call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:7, :), ms(3), n)
          call rgbmv(x((i + 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(8:8, :), ms(4), n)
          call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(9:9, :), ms(5), n)
  
          do i = 3, n - 2
              call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
              call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
              call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:7, :), ms(3), n)
              call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(8:8, :), ms(4), n)
              call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(9:9, :), ms(5), n)
          end do
  
          i = n - 1
          call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
          call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
          call rgbmv(x((i - 1) * n + 1:(i - 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:7, :), ms(3), n)
          call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(8:8, :), ms(4), n)
  
          i = n
          call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
          call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
          call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:7, :), ms(3), n)
      end subroutine rbbmv_o5
  
      subroutine rbbmv_o7(x, y, sign, blocks, ms, n)
          implicit none
  
          integer, parameter :: m = 7
          integer, intent(in) :: n
          integer, intent(in), dimension(m) :: ms
          real(sp), intent(in) :: sign
          real(sp), intent(in), dimension(n * n) :: x
          real(sp), intent(inout), dimension(n * n) :: y
          real(sp), intent(in), dimension(2 * m - 1, n) :: blocks
  
          integer :: i
  
          i = 1
          call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(4:10, :), ms(4), n)
          call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(11:11, :), ms(5), n)
          call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(12:12, :), ms(6), n)
          call rgbmv(x((i + 2) * n + 1:(i + 3) * n), y((i - 1) * n + 1:i * n), sign, blocks(13:13, :), ms(7), n)
  
          i = 2
          call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:3, :), ms(3), n)
          call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(4:10, :), ms(4), n)
          call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(11:11, :), ms(5), n)
          call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(12:12, :), ms(6), n)
          call rgbmv(x((i + 2) * n + 1:(i + 3) * n), y((i - 1) * n + 1:i * n), sign, blocks(13:13, :), ms(7), n)
  
          i = 3
          call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
          call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:3, :), ms(3), n)
          call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(4:10, :), ms(4), n)
          call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(11:11, :), ms(5), n)
          call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(12:12, :), ms(6), n)
          call rgbmv(x((i + 2) * n + 1:(i + 3) * n), y((i - 1) * n + 1:i * n), sign, blocks(13:13, :), ms(7), n)
  
          do i = 4, n - 3
              call rgbmv(x((i - 4) * n + 1:(i - 3) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
              call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
              call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:3, :), ms(3), n)
              call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(4:10, :), ms(4), n)
              call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(11:11, :), ms(5), n)
              call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(12:12, :), ms(6), n)
              call rgbmv(x((i + 2) * n + 1:(i + 3) * n), y((i - 1) * n + 1:i * n), sign, blocks(13:13, :), ms(7), n)
          end do
  
          i = n - 2
          call rgbmv(x((i - 4) * n + 1:(i - 3) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
          call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
          call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:3, :), ms(3), n)
          call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(4:10, :), ms(4), n)
          call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(11:11, :), ms(5), n)
          call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(12:12, :), ms(6), n)
  
          i = n - 1
          call rgbmv(x((i - 4) * n + 1:(i - 3) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
          call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
          call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:3, :), ms(3), n)
          call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(4:10, :), ms(4), n)
          call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(11:11, :), ms(5), n)
  
          i = n
          call rgbmv(x((i - 4) * n + 1:(i - 3) * n), y((i - 1) * n + 1:i * n), sign, blocks(1:1, :), ms(1), n)
          call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(2:2, :), ms(2), n)
          call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(3:3, :), ms(3), n)
          call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(4:10, :), ms(4), n)
      end subroutine rbbmv_o7
  
      subroutine rbbmv(x, y, sign, blocks, ms, m, n)
          implicit none
  
          integer, intent(in) :: m, n
  