diff --git a/resources/example.mtx b/resources/example.mtx
new file mode 100644
index 0000000..428295f
--- /dev/null
+++ b/resources/example.mtx
@@ -0,0 +1,51 @@
+%%MatrixMarket matrix coordinate real general
+25 25 49
+1 1 1.000e+00
+2 2 2.000e+00
+3 3 3.000e+00
+4 4 4.000e+00
+5 5 5.000e+00
+6 6 6.000e+00
+7 7 7.000e+00
+8 8 8.000e+00
+9 9 9.000e+00
+10 10 1.000e+01
+11 11 2.000e+01
+12 12 3.000e+01
+13 13 4.000e+01
+14 14 5.000e+01
+15 15 6.000e+01
+16 16 7.000e+01
+17 17 8.000e+01
+18 18 8.000e+01
+19 19 9.000e+01
+20 20 1.000e+02
+21 21 2.000e+02
+22 22 2.000e+02
+23 23 3.000e+02
+24 24 4.000e+02
+25 25 5.000e+02
+1 2 1.000e+00
+2 3 2.000e+00
+3 4 3.000e+00
+4 5 4.000e+00
+5 6 5.000e+00
+6 7 6.000e+00
+7 8 7.000e+00
+8 9 8.000e+00
+9 10 9.000e+00
+10 11 1.000e+01
+11 12 2.000e+01
+12 13 3.000e+01
+13 14 4.000e+01
+14 15 5.000e+01
+15 16 6.000e+01
+16 17 7.000e+01
+17 18 8.000e+01
+18 19 9.000e+01
+19 20 1.000e+01
+20 21 2.000e+01
+21 22 3.000e+01
+22 23 4.000e+01
+23 24 5.000e+01
+24 25 6.000e+01
diff --git a/scripts/build_run.sh b/scripts/build_run.sh
index 65400c9..ae8895e 100755
--- a/scripts/build_run.sh
+++ b/scripts/build_run.sh
@@ -17,3 +17,4 @@ mpirun -n 3 ./build/src/example4
 mpirun -n 3 ./build/src/example5
 mpirun -n 3 ./build/src/example6
 mpirun -n 3 ./build/src/example7
+mpirun -n 3 ./build/src/example8
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 7e03f09..fcda0fa 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -46,3 +46,8 @@ add_executable(example7 example7.cpp)
 
 target_compile_definitions(example7 INTERFACE DR_FORMAT)
 target_link_libraries(example7 DR::mpi fmt::fmt)
+
+add_executable(example8 example8.cpp)
+
+target_compile_definitions(example8 INTERFACE DR_FORMAT)
+target_link_libraries(example8 DR::mpi fmt::fmt)
diff --git a/src/example7.cpp b/src/example7.cpp
index 7330319..5a2fbb1 100644
--- a/src/example7.cpp
+++ b/src/example7.cpp
@@ -10,13 +10,12 @@ int main() {
   dr::mp::init(sycl::default_selector_v);
   using I = long;
   using V = double;
+
   dr::views::csr_matrix_view<V, I> local_data;
   auto root = 0;
-  auto n = 10;
-  auto up = 1;    // number of diagonals above main diagonal
-  auto down = up; // number of diagonals below main diagonal
   if (root == dr::mp::rank()) {
-    local_data = dr::generate_band_csr<V, I>(n, up, down);
+    auto source = "./resources/example.mtx";
+    local_data = dr::read_csr<double, long>(source);
   }
 
   dr::mp::distributed_sparse_matrix<
@@ -38,7 +37,7 @@ int main() {
   gemv(root, res, matrix, broadcasted_b);
 
   if (root == dr::mp::rank()) {
-    fmt::print("Band matrix {} x {} with bandwitch {}\n", n, n, up * 2);
+    fmt::print("Matrix imported from {}\n", "./resources/example.mtx");
     fmt::print("Input: ");
     for (auto x : b) {
       fmt::print("{} ", x);
diff --git a/src/example8.cpp b/src/example8.cpp
new file mode 100644
index 0000000..391cdaa
--- /dev/null
+++ b/src/example8.cpp
@@ -0,0 +1,74 @@
+// SPDX-FileCopyrightText: Intel Corporation
+//
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include <dr/mp.hpp>
+#include <fmt/core.h>
+#include <random>
+
+/* Sparse band matrix vector multiplication */
+int main() {
+  dr::mp::init(sycl::default_selector_v);
+  using I = long;
+  using V = double;
+  dr::views::csr_matrix_view<V, I> local_data;
+  auto root = 0;
+  if (root == dr::mp::rank()) {
+    auto size = 10;
+    auto nnz = 20;
+    auto colInd = new I[nnz];
+    auto rowInd = new I[size + 1];
+    auto values = new V[nnz];
+    std::uniform_real_distribution<double> unif(0, 1);
+    std::default_random_engine re;
+    for (auto i = 0; i <= size; i++) {
+      rowInd[i] = i * 2; // two elements per row
+    }
+    for (auto i = 0; i < nnz; i++) {
+      colInd[i] = (i % 2) * (std::max(i / 2, 1)); // column on 0 and diagonal (with additional entry in first row)
+      values[i] = unif(re);
+    }
+
+    local_data = dr::views::csr_matrix_view<V, I>(values, rowInd, colInd, {size, size}, nnz, root);
+  }
+
+  dr::mp::distributed_sparse_matrix<
+      V, I, dr::mp::MpiBackend,
+      dr::mp::csr_eq_distribution<V, I, dr::mp::MpiBackend>>
+      matrix(local_data, root);
+
+  std::vector<double> b;
+  b.reserve(matrix.shape().second);
+  std::vector<double> res(matrix.shape().first);
+  for (auto i = 0; i < matrix.shape().second; i++) {
+    b.push_back(i);
+  }
+
+  dr::mp::broadcasted_vector<double> broadcasted_b;
+  broadcasted_b.broadcast_data(matrix.shape().second, 0, b,
+                               dr::mp::default_comm());
+
+  gemv(root, res, matrix, broadcasted_b);
+
+  if (root == dr::mp::rank()) {
+    fmt::print("Matrix with {} x {} and number of non-zero entries equal to {} and entries:\n", matrix.shape().first, matrix.shape().second, matrix.size());
+    for (auto [i, v]: matrix) {
+      auto [n, m] = i;
+      fmt::print("Matrix entry <{}, {}, {}>\n", n, m, v);
+    }
+    fmt::print("Input: ");
+    for (auto x : b) {
+      fmt::print("{} ", x);
+    }
+    fmt::print("\n");
+    fmt::print("Matrix vector multiplication res: ");
+    for (auto x : res) {
+      fmt::print("{} ", x);
+    }
+    fmt::print("\n");
+  }
+
+  dr::mp::finalize();
+
+  return 0;
+}