anchor.diff:      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 :: rbbmv_o3, rbbmv_o5, rbbmv_o7, rbbmv, rgbmv
       public :: revervoir, revervoir_2d
       public :: hamiltonian, hamiltonian_2d
       public :: runge_kutta, runge_kutta_2d
  
positive.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, m, n)
          implicit none
  
          integer, intent(in) :: m, n
  
negative.code:         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, n)
          else if (m == 5) then
              call rbbmv_o5(x, y, sign, blocks, ms, n)
          else if (m == 7) then
              call rbbmv_o7(x, y, sign, blocks, ms, n)
          end if
      end subroutine rbbmv
  
  