diff --git a/include/hdf5_utils.hpp b/include/hdf5_utils.hpp
index 181e984e..d89c71c5 100644
--- a/include/hdf5_utils.hpp
+++ b/include/hdf5_utils.hpp
@@ -7,6 +7,7 @@
 
 #include "mfem.hpp"
 #include "hdf5.h"
+#include "linalg_utils.hpp"
 
 using namespace mfem;
 
@@ -15,6 +16,7 @@ namespace hdf5_utils
 
 inline hid_t GetType(int) { return (H5T_NATIVE_INT); }
 inline hid_t GetType(double) { return (H5T_NATIVE_DOUBLE); }
+inline hid_t GetType(bool) { return (H5T_NATIVE_HBOOL); }
 
 hid_t GetNativeType(hid_t type);
 
@@ -169,6 +171,9 @@ void WriteDataset(hid_t &source, std::string dataset, const DenseMatrix &value);
 void ReadDataset(hid_t &source, std::string dataset, DenseTensor &value);
 void WriteDataset(hid_t &source, std::string dataset, const DenseTensor &value);
 
+void ReadDataset(hid_t &source, std::string dataset, MatrixBlocks &value);
+void WriteDataset(hid_t &source, std::string dataset, const MatrixBlocks &value);
+
 inline bool pathExists(hid_t id, const std::string& path)
 {
   return H5Lexists(id, path.c_str(), H5P_DEFAULT) > 0;
diff --git a/include/linalg_utils.hpp b/include/linalg_utils.hpp
index e2aa5d7d..a033a11b 100644
--- a/include/linalg_utils.hpp
+++ b/include/linalg_utils.hpp
@@ -61,6 +61,17 @@ struct MatrixBlocks
 
    ~MatrixBlocks() { DeletePointers(blocks); }
 
+   void SetSize(const int w, const int h)
+   {
+      for (int i = 0; i < nrows; i++)
+         for (int j = 0; j < ncols; j++)
+            if (blocks(i, j)) delete blocks(i, j);
+
+      nrows = w; ncols = h;
+      blocks.SetSize(w, h);
+      blocks = NULL;
+   }
+
    SparseMatrix* operator()(const int i, const int j) const
    {
       assert((i >= 0) && (i < nrows));
diff --git a/src/hdf5_utils.cpp b/src/hdf5_utils.cpp
index c3a4b3b3..5c7b3361 100644
--- a/src/hdf5_utils.cpp
+++ b/src/hdf5_utils.cpp
@@ -241,7 +241,7 @@ BlockMatrix* ReadBlockMatrix(hid_t &source, std::string matrix_name,
    assert(grp_id >= 0);
 
    Array<int> size, row_offsets, col_offsets;
-   Array2D<int> zero_blocks;
+   Array2D<bool> zero_blocks;
    ReadDataset(grp_id, "size", size);
    ReadDataset(grp_id, "row_offsets", row_offsets);
    ReadDataset(grp_id, "col_offsets", col_offsets);
@@ -305,7 +305,7 @@ void WriteBlockMatrix(hid_t &source, std::string matrix_name, BlockMatrix* mat)
    WriteDataset(grp_id, "row_offsets", row_offsets);
    WriteDataset(grp_id, "col_offsets", col_offsets);
 
-   Array2D<int> zero_blocks(size[0], size[1]);
+   Array2D<bool> zero_blocks(size[0], size[1]);
    for (int i = 0; i < size[0]; i++)
       for (int j = 0; j < size[1]; j++)
          zero_blocks(i, j) = (mat->IsZeroBlock(i, j) || (mat->GetBlock(i,j).NumNonZeroElems() == 0));
@@ -418,4 +418,65 @@ void WriteDataset(hid_t &source, std::string dataset, const DenseTensor &value)
    assert(errf >= 0);
 }
 
+void ReadDataset(hid_t &source, std::string dataset, MatrixBlocks &value)
+{
+   herr_t errf = 0;
+   hid_t grp_id;
+   grp_id = H5Gopen2(source, dataset.c_str(), H5P_DEFAULT);
+   assert(grp_id >= 0);
+
+   int nrows, ncols;
+   ReadAttribute(grp_id, "nrows", nrows);
+   ReadAttribute(grp_id, "ncols", ncols);
+   value.SetSize(nrows, ncols);
+
+   Array2D<bool> zero_blocks;
+   ReadDataset(grp_id, "zero_blocks", zero_blocks);
+   assert(zero_blocks.NumRows() == nrows);
+   assert(zero_blocks.NumCols() == ncols);
+
+   std::string block_name;
+   for (int i = 0; i < nrows; i++)
+      for (int j = 0; j < ncols; j++)
+      {
+         if (zero_blocks(i, j)) continue;
+
+         block_name = "block_" + std::to_string(i) + "_" + std::to_string(j);  
+         value.blocks(i, j) = ReadSparseMatrix(grp_id, block_name);
+      }
+
+   errf = H5Gclose(grp_id);
+   assert(errf >= 0);
+}
+
+void WriteDataset(hid_t &source, std::string dataset, const MatrixBlocks &value)
+{
+   herr_t errf = 0;
+   hid_t grp_id;
+   grp_id = H5Gcreate(source, dataset.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+   assert(grp_id >= 0);
+
+   WriteAttribute(grp_id, "nrows", value.nrows);
+   WriteAttribute(grp_id, "ncols", value.ncols);
+
+   Array2D<bool> zero_blocks(value.nrows, value.ncols);
+   for (int i = 0; i < value.nrows; i++)
+      for (int j = 0; j < value.ncols; j++)
+         zero_blocks(i, j) = (value.blocks(i, j) || (value.blocks(i, j)->NumNonZeroElems() == 0));
+   WriteDataset(grp_id, "zero_blocks", zero_blocks);
+
+   std::string block_name;
+   for (int i = 0; i < value.nrows; i++)
+      for (int j = 0; j < value.ncols; j++)
+      {
+         block_name = "block_" + std::to_string(i) + "_" + std::to_string(j);
+         if (zero_blocks(i, j)) continue;
+
+         WriteSparseMatrix(grp_id, block_name, value.blocks(i, j));
+      }
+
+   errf = H5Gclose(grp_id);
+   assert(errf >= 0);
+}
+
 }