diff --git a/raptor/CMakeLists.txt b/raptor/CMakeLists.txt
index 353ca2aa..47724495 100644
--- a/raptor/CMakeLists.txt
+++ b/raptor/CMakeLists.txt
@@ -19,7 +19,7 @@ endif()
 
 add_library(raptor ${core_SOURCES} ${core_HEADERS} 
 	${gallery_SOURCES} ${gallery_HEADERS} ${ext_gallery_HEADERS} 
-	${util_SOURCES} ${util_HEADERS}
+	${util_SOURCES} ${util_HEADERS} ${root_util_HEADERS}
 	${par_SOURCES} strength.cpp 
 	${ruge_stuben_SOURCES} ${ruge_stuben_HEADERS}
         ${aggregation_SOURCES} ${aggregation_HEADERS}
@@ -46,6 +46,7 @@ install(FILES ${core_HEADERS} DESTINATION "include/raptor/core")
 install(FILES ${gallery_HEADERS}  DESTINATION "include/raptor/gallery")
 install(FILES ${ext_gallery_HEADERS} DESTINATION "include/raptor/gallery/external")
 install(FILES ${util_HEADERS} DESTINATION "include/raptor/util/linalg")
+install(FILES ${root_util_HEADERS} DESTINATION "include/raptor/util")
 install(FILES ${ext_util_HEADERS} DESTINATION "include/raptor/util/linalg/external")
 install(FILES ${ruge_stuben_HEADERS} DESTINATION "include/raptor/ruge_stuben")
 install(FILES ${aggregation_HEADERS} DESTINATION "include/raptor/aggregation")
diff --git a/raptor/core/matrix_traits.hpp b/raptor/core/matrix_traits.hpp
index b9604bd6..5537c238 100644
--- a/raptor/core/matrix_traits.hpp
+++ b/raptor/core/matrix_traits.hpp
@@ -11,6 +11,7 @@ using is_bsr_or_csr = std::enable_if_t<std::is_same_v<T, ParCSRMatrix> ||
 
 template <class T> struct is_bsr : std::false_type {};
 template <> struct is_bsr<ParBSRMatrix> : std::true_type {};
+template <> struct is_bsr<BSRMatrix> : std::true_type {};
 template <class T> inline constexpr bool is_bsr_v = is_bsr<T>::value;
 
 inline BSRMatrix & bsr_cast(Matrix &mat) { return dynamic_cast<BSRMatrix &>(mat); }
diff --git a/raptor/util/CMakeLists.txt b/raptor/util/CMakeLists.txt
index 8dd041cb..cc5c772f 100644
--- a/raptor/util/CMakeLists.txt
+++ b/raptor/util/CMakeLists.txt
@@ -5,12 +5,16 @@ add_subdirectory(linalg)
 
 add_subdirectory(linalg/external)
 
+set(root_util_HEADERS
+  util/writer.hpp
+  PARENT_SCOPE)
 set(util_HEADERS
     ${linalg_HEADERS}
     PARENT_SCOPE)
 set(util_SOURCES
-    ${linalg_SOURCES}
-    PARENT_SCOPE)
+  util/writer.cpp
+  ${linalg_SOURCES}
+  PARENT_SCOPE)
 
 set(ext_util_HEADERS
     ${ext_linalg_HEADERS}
diff --git a/raptor/util/writer.cpp b/raptor/util/writer.cpp
new file mode 100644
index 00000000..2af2a429
--- /dev/null
+++ b/raptor/util/writer.cpp
@@ -0,0 +1,119 @@
+#include <array>
+#include <sstream>
+#include <fstream>
+
+#include "raptor/core/matrix_traits.hpp"
+#include "raptor/util/writer.hpp"
+
+namespace raptor {
+
+namespace impl {
+void write_header(const char *fname, const raptor::ParCSRMatrix & mat) {
+	std::ofstream ofile(std::string(fname) + ".hdr", std::ios_base::binary);
+
+	int nprocs; MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+	std::array<int, 4> buf{
+		0, // csr
+		mat.global_num_rows,
+		mat.global_num_cols,
+		nprocs};
+
+	ofile.write(reinterpret_cast<const char*>(buf.data()), sizeof(int)*buf.size());
+}
+
+
+void write_header(const char *fname, const raptor::ParBSRMatrix & mat) {
+	std::ofstream ofile(std::string(fname) + ".hdr", std::ios_base::binary);
+
+	int nprocs; MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+	auto & diag = dynamic_cast<raptor::BSRMatrix&>(*mat.on_proc);
+	std::array<int, 6> buf{
+		1, // bsr
+		mat.global_num_rows,
+		mat.global_num_cols,
+		nprocs,
+		diag.b_rows,
+		diag.b_cols};
+
+	ofile.write(reinterpret_cast<const char*>(buf.data()), sizeof(int)*buf.size());
+}
+
+
+void write_rowptr(std::ostream & out,
+                  const raptor::Matrix & diag,
+                  const raptor::Matrix & offd) {
+	for (int i = 0; i < diag.n_rows + 1; ++i) {
+		int ptr = diag.idx1[i] + offd.idx1[i];
+		out.write(reinterpret_cast<const char *>(&ptr), sizeof(int));
+	}
+}
+
+
+template<class T>
+void write_rows(std::ostream & out,
+                const T & diag, const T & offd,
+                const std::vector<int> & diag_colmap, const std::vector<int> & offd_colmap) {
+	out.write(reinterpret_cast<const char*>(&diag.n_rows), sizeof(diag.n_rows));
+	write_rowptr(out, diag, offd);
+
+	for (int i = 0; i < diag.n_rows; ++i) {
+		auto write_colinds = [&](const T & mat, const std::vector<int> & colmap) {
+			for (int j = mat.idx1[i]; j < mat.idx1[i + 1]; ++j) {
+				int gcol = colmap[mat.idx2[j]];
+				out.write(reinterpret_cast<const char *>(&gcol), sizeof(int));
+			}
+		};
+
+		write_colinds(diag, diag_colmap);
+		write_colinds(offd, offd_colmap);
+	}
+	for (int i = 0; i < diag.n_rows; ++i) {
+		auto write_values = [&](const T & mat) {
+			if constexpr (is_bsr_v<T>) {
+				for (int j = mat.idx1[i]; j < mat.idx1[i + 1]; ++j) {
+					out.write(reinterpret_cast<const char *>(mat.block_vals[j]),
+					          mat.b_size * sizeof(double));
+				}
+			} else {
+				out.write(reinterpret_cast<const char *>(&mat.vals[mat.idx1[i]]),
+				          (mat.idx1[i + 1] - mat.idx1[i]) * sizeof(double));
+			}
+		};
+
+		write_values(diag);
+		write_values(offd);
+	}
+}
+
+
+std::ofstream get_rank_file(const char * fname) {
+	int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+	std::ostringstream rank_file;
+	rank_file << fname << '.' << rank;
+	return std::ofstream(rank_file.str(), std::ios_base::binary);
+}
+
+
+template <class T, is_bsr_or_csr<T> = true>
+void write(const char * fname, const T & mat) {
+	write_header(fname, mat);
+
+	auto cast = [](const auto & m) -> const auto & {
+		return dynamic_cast<const sequential_matrix_t<T> &>(m);
+	};
+
+	auto ofile = get_rank_file(fname);
+	write_rows(ofile, cast(*mat.on_proc), cast(*mat.off_proc),
+	           mat.on_proc_column_map, mat.off_proc_column_map);
+}
+} // namespace impl
+
+void write(const char *fname, const ParCSRMatrix &mat) {
+	impl::write(fname, mat);
+}
+
+void write(const char *fname, const ParBSRMatrix & mat) {
+	impl::write(fname, mat);
+}
+
+} // namespace raptor
diff --git a/raptor/util/writer.hpp b/raptor/util/writer.hpp
new file mode 100644
index 00000000..09ad68e9
--- /dev/null
+++ b/raptor/util/writer.hpp
@@ -0,0 +1,13 @@
+#ifndef RAPTOR_UTIL_WRITER_HPP
+#define RAPTOR_UTIL_WRITER_HPP
+
+#include "raptor/core/par_matrix.hpp"
+
+namespace raptor
+{
+
+void write(const char *fname, const ParCSRMatrix &mat);
+void write(const char *fname, const ParBSRMatrix &mat);
+
+}
+#endif