apply_matrix Subroutine

   subroutine apply_matrix(mesh, diag, coeff, nb, x, ax)
      type(mesh_t), intent(in) :: mesh
      real(rk), intent(in) :: diag(:)
      real(rk), intent(in) :: coeff(:,:)
      integer, intent(in) :: nb(:,:)
      real(rk), intent(in) :: x(:)
      real(rk), intent(out) :: ax(:)

      integer :: c, lf, neighbor

      do c = 1, mesh%ncells
         ax(c) = diag(c) * x(c)
         do lf = 1, mesh%ncell_faces(c)
            neighbor = nb(lf, c)
            if (neighbor > 0) then
               ax(c) = ax(c) - coeff(lf, c) * x(neighbor)
            end if
         end do
      end do
   end subroutine apply_matrix