Handle symmetric with binsparse read/write

BenBrock · BenBrock · commit 3cbb9be40db1 · 2023-11-27T21:29:38.000-08:00
diff --git a/examples/convert_binsparse.cpp b/examples/convert_binsparse.cpp
@@ -14,7 +14,8 @@ void convert_to_binsparse(std::string input_file, std::string output_file,
         T, I, binsparse::__detail::csr_matrix_owning<T, I>>(input_file);
     binsparse::csr_matrix<T, I> matrix{
         x.values().data(),      x.colind().data(),      x.rowptr().data(),
-        std::get<0>(x.shape()), std::get<1>(x.shape()), I(x.size())};
+        std::get<0>(x.shape()), std::get<1>(x.shape()), I(x.size()),
+        x.structure()};
     binsparse::write_csr_matrix(output_file, matrix, user_keys);
     std::cout << "Writing to binsparse file " << output_file << " using "
               << format << " format...\n";
@@ -23,7 +24,8 @@ void convert_to_binsparse(std::string input_file, std::string output_file,
         T, I, binsparse::__detail::coo_matrix_owning<T, I>>(input_file);
     binsparse::coo_matrix<T, I> matrix{
         x.values().data(),      x.rowind().data(),      x.colind().data(),
-        std::get<0>(x.shape()), std::get<1>(x.shape()), I(x.size())};
+        std::get<0>(x.shape()), std::get<1>(x.shape()), I(x.size()),
+        x.structure()};
     binsparse::write_coo_matrix(output_file, matrix, user_keys);
     std::cout << "Writing to binsparse file " << output_file << " using "
               << format << " format...\n";
diff --git a/include/binsparse/binsparse.hpp b/include/binsparse/binsparse.hpp
@@ -34,6 +34,11 @@ void write_dense_matrix(std::string fname, dense_matrix<T, I, Order> m,
   j["binsparse"]["nnz"] = m.m * m.n;
   j["binsparse"]["data_types"]["values"] = type_info<T>::label();
 
+  if (m.structure != general) {
+    j["binsparse"]["structure"] =
+        __detail::get_structure_name(m.structure).value();
+  }
+
   for (auto&& v : user_keys.items()) {
     j[v.key()] = v.value();
   }
@@ -67,7 +72,13 @@ auto read_dense_matrix(std::string fname, Allocator&& alloc = Allocator{}) {
 
   auto values = hdf5_tools::read_dataset<T>(f, "values", alloc);
 
-  return dense_matrix<T, I, Order>{values.data(), nrows, ncols};
+  structure_t structure = general;
+
+  if (binsparse_metadata.contains("structure")) {
+    structure = __detail::parse_structure(binsparse_metadata["structure"]);
+  }
+
+  return dense_matrix<T, I, Order>{values.data(), nrows, ncols, structure};
 }
 
 // CSR Format
@@ -96,6 +107,11 @@ void write_csr_matrix(std::string fname, csr_matrix<T, I> m,
   j["binsparse"]["data_types"]["indices_1"] = type_info<I>::label();
   j["binsparse"]["data_types"]["values"] = type_info<T>::label();
 
+  if (m.structure != general) {
+    j["binsparse"]["structure"] =
+        __detail::get_structure_name(m.structure).value();
+  }
+
   for (auto&& v : user_keys.items()) {
     j[v.key()] = v.value();
   }
@@ -130,8 +146,14 @@ csr_matrix<T, I> read_csr_matrix(std::string fname, Allocator&& alloc) {
   auto colind = hdf5_tools::read_dataset<I>(f, "indices_1", i_alloc);
   auto row_ptr = hdf5_tools::read_dataset<I>(f, "pointers_to_1", i_alloc);
 
-  return csr_matrix<T, I>{values.data(), colind.data(), row_ptr.data(),
-                          nrows,         ncols,         nnz};
+  structure_t structure = general;
+
+  if (binsparse_metadata.contains("structure")) {
+    structure = __detail::parse_structure(binsparse_metadata["structure"]);
+  }
+
+  return csr_matrix<T, I>{values.data(), colind.data(), row_ptr.data(), nrows,
+                          ncols,         nnz,           structure};
 }
 
 template <typename T, typename I>
@@ -165,6 +187,11 @@ void write_csc_matrix(std::string fname, csc_matrix<T, I> m,
   j["binsparse"]["data_types"]["indices_1"] = type_info<I>::label();
   j["binsparse"]["data_types"]["values"] = type_info<T>::label();
 
+  if (m.structure != general) {
+    j["binsparse"]["structure"] =
+        __detail::get_structure_name(m.structure).value();
+  }
+
   for (auto&& v : user_keys.items()) {
     j[v.key()] = v.value();
   }
@@ -199,8 +226,14 @@ csc_matrix<T, I> read_csc_matrix(std::string fname, Allocator&& alloc) {
   auto rowind = hdf5_tools::read_dataset<I>(f, "indices_1", i_alloc);
   auto col_ptr = hdf5_tools::read_dataset<I>(f, "pointers_to_1", i_alloc);
 
-  return csc_matrix<T, I>{values.data(), rowind.data(), col_ptr.data(),
-                          nrows,         ncols,         nnz};
+  structure_t structure = general;
+
+  if (binsparse_metadata.contains("structure")) {
+    structure = __detail::parse_structure(binsparse_metadata["structure"]);
+  }
+
+  return csc_matrix<T, I>{values.data(), rowind.data(), col_ptr.data(), nrows,
+                          ncols,         nnz,           structure};
 }
 
 template <typename T, typename I>
@@ -234,6 +267,11 @@ void write_coo_matrix(std::string fname, coo_matrix<T, I> m,
   j["binsparse"]["data_types"]["indices_1"] = type_info<I>::label();
   j["binsparse"]["data_types"]["values"] = type_info<T>::label();
 
+  if (m.structure != general) {
+    j["binsparse"]["structure"] =
+        __detail::get_structure_name(m.structure).value();
+  }
+
   for (auto&& v : user_keys.items()) {
     j[v.key()] = v.value();
   }
@@ -270,8 +308,14 @@ coo_matrix<T, I> read_coo_matrix(std::string fname, Allocator&& alloc) {
   auto rows = hdf5_tools::read_dataset<I>(f, "indices_0", i_alloc);
   auto cols = hdf5_tools::read_dataset<I>(f, "indices_1", i_alloc);
 
-  return coo_matrix<T, I>{values.data(), rows.data(), cols.data(),
-                          nrows,         ncols,       nnz};
+  structure_t structure = general;
+
+  if (binsparse_metadata.contains("structure")) {
+    structure = __detail::parse_structure(binsparse_metadata["structure"]);
+  }
+
+  return coo_matrix<T, I>{values.data(), rows.data(), cols.data(), nrows,
+                          ncols,         nnz,         structure};
 }
 
 template <typename T, typename I>
diff --git a/include/binsparse/detail.hpp b/include/binsparse/detail.hpp
@@ -27,6 +27,32 @@ inline std::string unalias_format(const std::string& format) {
   }
 }
 
+inline std::optional<std::string> get_structure_name(structure_t structure) {
+  if (structure == general) {
+    return {};
+  } else if (structure == symmetric) {
+    return "symmetric_lower";
+  } else if (structure == skew_symmetric) {
+    return "skew_symmetric_lower";
+  } else if (structure == hermitian) {
+    return "hermitian";
+  } else {
+    throw std::runtime_error("get_structure_name: unknown structure");
+  }
+}
+
+inline structure_t parse_structure(const std::string& structure) {
+  if (structure == "symmetric_lower") {
+    return symmetric;
+  } else if (structure == "skew_symmetric_lower") {
+    return skew_symmetric;
+  } else if (structure == "hermitian") {
+    return hermitian;
+  } else {
+    throw std::runtime_error("parse_structure: unsupported structure");
+  }
+}
+
 } // namespace __detail
 
 } // namespace binsparse
diff --git a/include/binsparse/matrix_market/matrix_market_read.hpp b/include/binsparse/matrix_market/matrix_market_read.hpp
@@ -11,7 +11,8 @@ namespace __detail {
 template <typename T, typename I>
 class csr_matrix_owning {
 public:
-  csr_matrix_owning(std::tuple<I, I> shape) : shape_(shape) {}
+  csr_matrix_owning(std::tuple<I, I> shape, structure_t structure = general)
+      : shape_(shape), structure_(structure) {}
 
   auto values() {
     return std::ranges::views::all(values_);
@@ -80,17 +81,23 @@ class csr_matrix_owning {
     return values_.size();
   }
 
+  auto structure() const {
+    return structure_;
+  }
+
 private:
   std::tuple<I, I> shape_;
   std::vector<T> values_;
   std::vector<I> rowptr_;
   std::vector<I> colind_;
+  structure_t structure_;
 };
 
 template <typename T, typename I>
 class coo_matrix_owning {
 public:
-  coo_matrix_owning(std::tuple<I, I> shape) : shape_(shape) {}
+  coo_matrix_owning(std::tuple<I, I> shape, structure_t structure = general)
+      : shape_(shape), structure_(structure) {}
 
   auto values() {
     return std::ranges::views::all(values_);
@@ -144,11 +151,16 @@ class coo_matrix_owning {
     return values_.size();
   }
 
+  auto structure() const {
+    return structure_;
+  }
+
 private:
   std::tuple<I, I> shape_;
   std::vector<T> values_;
   std::vector<I> rowind_;
   std::vector<I> colind_;
+  structure_t structure_;
 };
 
 /// Read in the Matrix Market file at location `file_path` and
@@ -170,9 +182,8 @@ inline MatrixType mmread(std::string file_path, bool one_indexed = true) {
   std::string buf;
 
   // Make sure the file is matrix market matrix, coordinate, and check whether
-  // it is symmetric. If the matrix is symmetric, non-diagonal elements will
-  // be inserted in both (i, j) and (j, i).  Error out if skew-symmetric or
-  // Hermitian.
+  // it is symmetric. If the matrix is symmetric.
+  // Error out if skew-symmetric or Hermitian.
   std::getline(f, buf);
   std::istringstream ss(buf);
   std::string item;
@@ -200,11 +211,11 @@ inline MatrixType mmread(std::string file_path, bool one_indexed = true) {
   }
   // TODO: do something with real vs. integer vs. pattern?
   ss >> item;
-  bool symmetric;
+  structure_t structure;
   if (item == "general") {
-    symmetric = false;
+    structure = general;
   } else if (item == "symmetric") {
-    symmetric = true;
+    structure = symmetric;
   } else {
     throw std::runtime_error(file_path + " has an unsupported matrix type");
   }
@@ -224,18 +235,11 @@ inline MatrixType mmread(std::string file_path, bool one_indexed = true) {
   ss.str(buf);
   ss >> m >> n >> nnz;
 
-  // NOTE for symmetric matrices: `nnz` holds the number of stored values in
-  // the matrix market file, while `matrix.nnz_` will hold the total number of
-  // stored values (including "mirrored" symmetric values).
-  MatrixType m_out({m, n});
+  MatrixType m_out({m, n}, structure);
 
   using coo_type = std::vector<std::tuple<std::tuple<I, I>, T>>;
   coo_type matrix;
-  if (symmetric) {
-    matrix.reserve(2 * nnz);
-  } else {
-    matrix.reserve(nnz);
-  }
+  matrix.reserve(nnz);
 
   size_type c = 0;
   while (std::getline(f, buf)) {
@@ -260,10 +264,6 @@ inline MatrixType mmread(std::string file_path, bool one_indexed = true) {
 
     matrix.push_back({{i, j}, v});
 
-    if (symmetric && i != j) {
-      matrix.push_back({{j, i}, v});
-    }
-
     c++;
     if (c > nnz) {
       throw std::runtime_error("read_MatrixMarket: error reading Matrix Market "
