Add symmetry generalization in application bindings

include/fast_matrix_market/app/Blaze.hpp

@@ -34,8 +34,7 @@ namespace fast_matrix_market {
         std::vector<IT> cols(storage_nnz);
         std::vector<VT> vals(storage_nnz);
 
-        auto handler = triplet_parse_handler(rows.begin(), cols.begin(), vals.begin());
-        read_matrix_market_body(instream, header, handler, default_pattern_value, options);
+        read_matrix_market_body_triplet(instream, header, rows, cols, vals, default_pattern_value, options);
 
         // Set the values into the matrix from the triplets.
         // The matrix needs to be constructed row-by-row (if row-major) or column-by-column (if col-major), in order.

include/fast_matrix_market/app/CXSparse.hpp

@@ -14,11 +14,17 @@ namespace fast_matrix_market {
     void read_matrix_market_cxsparse(std::istream &instream,
                                     CS **cs,
                                     ALLOC spalloc,
-                                    const read_options& options = {}) {
+                                    read_options options = {}) {
 
         matrix_market_header header;
         read_header(instream, header);
 
+        // Sanitize symmetry generalization settings
+        if (options.generalize_symmetry && options.generalize_symmetry_app) {
+            // cs_compress drops zero elements
+            options.generalize_coordinate_diagnonal_values = read_options::ExtraZeroElement;
+        }
+
         // allocate
         *cs = spalloc(header.nrows, header.ncols, get_storage_nnz(header, options),
                       header.field == pattern ? 0 : 1, // pattern field means do not allocate values

include/fast_matrix_market/app/Eigen.hpp

@@ -23,7 +23,7 @@ namespace fast_matrix_market {
     void read_matrix_market_eigen(std::istream &instream,
                                   matrix_market_header &header,
                                   SparseType& mat,
-                                  const read_options& options = {},
+                                  read_options options = {},
                                   typename SparseType::Scalar default_pattern_value = 1) {
 
         typedef typename SparseType::Scalar Scalar;
@@ -33,6 +33,12 @@ namespace fast_matrix_market {
         read_header(instream, header);
         mat.resize(header.nrows, header.ncols);
 
+        // Sanitize symmetry generalization settings
+        if (options.generalize_symmetry && options.generalize_symmetry_app) {
+            // setFromTriplets() drops zero elements
+            options.generalize_coordinate_diagnonal_values = read_options::ExtraZeroElement;
+        }
+
         // read into tuples
         std::vector<Triplet> elements;
         elements.resize(get_storage_nnz(header, options));

include/fast_matrix_market/app/GraphBLAS.hpp

@@ -550,13 +550,19 @@ namespace fast_matrix_market {
     void read_body_graphblas_coordinate(std::istream &instream,
                                         matrix_market_header &header,
                                         GrB_Matrix mat,
-                                        const read_options& options) {
+                                        read_options options) {
         size_t storage_nnz = get_storage_nnz(header, options);
 
         // Read into triplets
         std::vector<GrB_Index> rows(storage_nnz);
         std::vector<GrB_Index> cols(storage_nnz);
 
+        // Sanitize symmetry generalization settings
+        if (options.generalize_symmetry && options.generalize_symmetry_app) {
+            // Matrix construction drops zero elements
+            options.generalize_coordinate_diagnonal_values = read_options::ExtraZeroElement;
+        }
+
 #if FMM_GXB_BUILD_SCALAR
         if (header.field == pattern) {
             // read the indices

include/fast_matrix_market/app/triplet.hpp

@@ -36,6 +36,74 @@ namespace fast_matrix_market {
 #define triplet_write_vector typename
 #endif
 
+    /**
+     * Generalize symmetry of triplet.
+     *
+     * Does not duplicate diagonal elements.
+     */
+    template <typename IVEC, typename VVEC>
+    void generalize_symmetry_triplet(IVEC& rows, IVEC& cols, VVEC& values, const symmetry_type& symmetry) {
+        if (symmetry == general) {
+            return;
+        }
+
+        std::size_t num_diagonal_elements = 0;
+
+        // count how many diagonal elements there are (these do not get duplicated)
+        for (std::size_t i = 0; i < rows.size(); ++i) {
+            if (rows[i] == cols[i]) {
+                ++num_diagonal_elements;
+            }
+        }
+
+        // resize vectors
+        auto orig_size = rows.size();
+        auto new_size = 2*orig_size - num_diagonal_elements;
+        rows.resize(new_size);
+        cols.resize(new_size);
+        values.resize(new_size);
+
+        // fill in the new values
+        auto row_iter = rows.begin() + orig_size;
+        auto col_iter = cols.begin() + orig_size;
+        auto val_iter = values.begin() + orig_size;
+        for (std::size_t i = 0; i < orig_size; ++i) {
+            if (rows[i] == cols[i]) {
+                continue;
+            }
+
+            *row_iter = cols[i];
+            *col_iter = rows[i];
+            *val_iter = get_symmetric_value<typename VVEC::value_type>(values[i], symmetry);
+
+            ++row_iter; ++col_iter; ++val_iter;
+        }
+    }
+
+    template <triplet_read_vector IVEC, triplet_read_vector VVEC, typename T>
+    void read_matrix_market_body_triplet(std::istream &instream,
+                                         const matrix_market_header& header,
+                                         IVEC& rows, IVEC& cols, VVEC& values,
+                                         T pattern_value,
+                                         read_options options = {}) {
+        bool app_generalize = false;
+        if (options.generalize_symmetry && options.generalize_symmetry_app) {
+            app_generalize = true;
+            options.generalize_symmetry = false;
+        }
+
+        rows.resize(get_storage_nnz(header, options));
+        cols.resize(get_storage_nnz(header, options));
+        values.resize(get_storage_nnz(header, options));
+
+        auto handler = triplet_parse_handler(rows.begin(), cols.begin(), values.begin());
+        read_matrix_market_body(instream, header, handler, pattern_value, options);
+
+        if (app_generalize) {
+            generalize_symmetry_triplet(rows, cols, values, header.symmetry);
+        }
+    }
+
     /**
      * Read a Matrix Market file into a triplet (i.e. row, column, value vectors).
      */
@@ -44,16 +112,10 @@ namespace fast_matrix_market {
                                     matrix_market_header& header,
                                     IVEC& rows, IVEC& cols, VVEC& values,
                                     const read_options& options = {}) {
-        using VT = typename std::iterator_traits<decltype(values.begin())>::value_type;
-
         read_header(instream, header);
 
-        rows.resize(get_storage_nnz(header, options));
-        cols.resize(get_storage_nnz(header, options));
-        values.resize(get_storage_nnz(header, options));
-
-        auto handler = triplet_parse_handler(rows.begin(), cols.begin(), values.begin());
-        read_matrix_market_body(instream, header, handler, pattern_default_value((const VT*)nullptr), options);
+        using VT = typename std::iterator_traits<decltype(values.begin())>::value_type;
+        read_matrix_market_body_triplet(instream, header, rows, cols, values, pattern_default_value((const VT*)nullptr), options);
     }
 
     /**

include/fast_matrix_market/fast_matrix_market.hpp

@@ -119,6 +119,10 @@ namespace fast_matrix_market {
         return !o;
     }
 
+    inline bool negate(const std::vector<bool>::reference o) {
+        return !o;
+    }
+
     template <typename T>
     T negate(const T& o) {
         return std::negate<T>()(o);

include/fast_matrix_market/read_body.hpp

@@ -104,6 +104,24 @@ namespace fast_matrix_market {
     // Chunks
     ///////////////////////////////////////////////////////////////////
 
+    template <typename RET, typename T>
+    RET get_symmetric_value(const T& v, const symmetry_type& symmetry) {
+        switch (symmetry) {
+            case symmetric:
+                return v;
+            case skew_symmetric:
+                if constexpr (std::is_unsigned_v<T>) {
+                    throw invalid_argument("Cannot load skew-symmetric matrix into unsigned value type.");
+                } else {
+                    return negate(v);
+                }
+            case hermitian:
+                return complex_conjugate(v);
+            case general:
+                return v;
+        }
+    }
+
     template<typename HANDLER, typename IT, typename VT>
     void generalize_symmetry_coordinate(HANDLER& handler,
                                         const matrix_market_header &header,
@@ -445,11 +463,9 @@ namespace fast_matrix_market {
         }
 #endif
 
-        // Verify generalize symmetry is compatible with this file.
-        if (header.symmetry != general && options.generalize_symmetry) {
-            if (header.object != matrix) {
-                throw invalid_mm("Invalid Symmetry: vectors cannot have symmetry. Set generalize_symmetry=false to disregard this symmetry.");
-            }
+        // Sanity check input
+        if (header.object == vector && header.symmetry != general) {
+            throw invalid_mm("Vectors cannot have symmetry.");
         }
 
         if (header.format == array && header.field == pattern) {

include/fast_matrix_market/types.hpp

@@ -88,6 +88,12 @@ namespace fast_matrix_market {
          */
         bool generalize_symmetry = true;
 
+        /**
+         * If true, perform symmetry generalization in the application binding as a post-processing step.
+         * If supported this method can avoid extra diagonal elements.
+         */
+        bool generalize_symmetry_app = true;
+
         /**
          * Generalize Symmetry:
          * How to handle a value on the diagonal of a symmetric coordinate matrix.

tests/armadillo_test.cpp

@@ -117,3 +117,18 @@ TYPED_TEST(ArmadilloTest, SmallMatrices) {
         }
     }
 }
+
+TEST(ArmadilloTest, Symmetric) {
+    arma::SpMat<double> sym, expected;
+
+    {
+        std::ifstream f(kTestMatrixDir + "symmetry/coordinate_symmetric_row.mtx");
+        fast_matrix_market::read_matrix_market_arma(f, sym);
+    }
+    {
+        std::ifstream f(kTestMatrixDir + "symmetry/coordinate_symmetric_row_general.mtx");
+        fast_matrix_market::read_matrix_market_arma(f, expected);
+    }
+    EXPECT_EQ(expected.n_nonzero, sym.n_nonzero);
+    EXPECT_TRUE(arma::approx_equal(expected, sym, "absdiff", 1e-6));
+}

tests/basic_test.cpp

@@ -516,6 +516,7 @@ struct symmetry_problem {
     std::string symmetric;
     std::string general;
     std::string general_dup;
+    std::string general_zero;
 
     bool operator<(const symmetry_problem& rhs) const {
         return symmetric < rhs.symmetric;
@@ -547,6 +548,7 @@ class SymmetrySuite : public testing::TestWithParam<symmetry_problem> {
             p.symmetric = kSymmetrySubdir + std::regex_replace(filename, std::regex("_general"), "");
             p.general = kSymmetrySubdir + filename;
             p.general_dup = kSymmetrySubdir + std::regex_replace(filename, std::regex("_general"), "_general_dup");
+            p.general_zero = kSymmetrySubdir + std::regex_replace(filename, std::regex("_general"), "_general_zero");
             ret.push_back(p);
         }
 
@@ -578,32 +580,41 @@ class SymmetrySuite : public testing::TestWithParam<symmetry_problem> {
 
 using SymMat = triplet_matrix<int64_t, std::complex<double>>;
 TEST_P(SymmetrySuite, SmallTripletCoordinate) {
-    SymMat symmetric, sym_zero, sym_dup, general_zero, general_dup;
+    SymMat symmetric, sym_zero, sym_dup, sym_app, general_zero, general_dup, general_app;
 
     fast_matrix_market::read_options ro_no_gen{};
     ro_no_gen.generalize_symmetry = false;
 
     fast_matrix_market::read_options ro_gen_zero{};
     ro_gen_zero.generalize_symmetry = true;
+    ro_gen_zero.generalize_symmetry_app = false;
     ro_gen_zero.generalize_coordinate_diagnonal_values = fast_matrix_market::read_options::ExtraZeroElement;
 
     fast_matrix_market::read_options ro_gen_dup{};
     ro_gen_dup.generalize_symmetry = true;
+    ro_gen_dup.generalize_symmetry_app = false;
     ro_gen_dup.generalize_coordinate_diagnonal_values = fast_matrix_market::read_options::DuplicateElement;
 
+    fast_matrix_market::read_options ro_gen_app{};
+    ro_gen_app.generalize_symmetry = true;
+    ro_gen_app.generalize_symmetry_app = true;
+
     symmetry_problem p = GetParam();
     read_triplet_file(p.symmetric, symmetric, ro_no_gen);
     read_triplet_file(p.symmetric, sym_zero, ro_gen_zero);
     read_triplet_file(p.symmetric, sym_dup, ro_gen_dup);
-    read_triplet_file(p.general, general_zero, ro_no_gen);
+    read_triplet_file(p.symmetric, sym_app, ro_gen_app);
+    read_triplet_file(p.general_zero, general_zero, ro_no_gen);
     read_triplet_file(p.general_dup, general_dup, ro_no_gen);
+    read_triplet_file(p.general, general_app, ro_no_gen);
 
     EXPECT_EQ(symmetric.nrows, sym_zero.nrows);
     EXPECT_EQ(symmetric.ncols, sym_zero.ncols);
     EXPECT_EQ(symmetric.vals.size() * 2, sym_zero.vals.size());
     EXPECT_EQ(sym_dup.vals.size(), sym_zero.vals.size());
     EXPECT_EQ(sym_zero, general_zero);
     EXPECT_EQ(sym_dup, general_dup);
+    EXPECT_EQ(sym_app, general_app);
 }
 
 class SymmetryTripletArraySuite : public SymmetrySuite {};

tests/blaze_test.cpp

@@ -187,6 +187,22 @@ TYPED_TEST(BlazeMatrixTest, SmallMatrices) {
     }
 }
 
+
+TEST(BlazeMatrixTest, Symmetric) {
+    blaze::CompressedMatrix<double, blaze::rowMajor> sym, expected;
+
+    {
+        std::ifstream f(kTestMatrixDir + "symmetry/coordinate_symmetric_row.mtx");
+        fast_matrix_market::read_matrix_market_blaze(f, sym);
+    }
+    {
+        std::ifstream f(kTestMatrixDir + "symmetry/coordinate_symmetric_row_general.mtx");
+        fast_matrix_market::read_matrix_market_blaze(f, expected);
+    }
+    EXPECT_TRUE(is_equal(sym, expected));
+}
+
+
 /////////////////////////////////////////////////////////
 ///// Vectors
 /////////////////////////////////////////////////////////

tests/cxsparse_test.cpp

@@ -80,12 +80,19 @@ TEST_P(CSCTest, ReadSmall) {
 }
 
 INSTANTIATE_TEST_SUITE_P(
-        CSCTest,
+        CXSparse,
         CSCTest,
         ::testing::Values("eye3.mtx"
                           , "row_3by4.mtx"
                           , "kepner_gilbert_graph.mtx"
                           , "vector_array.mtx"
                           , "vector_coordinate.mtx"
                           , "eye3_pattern.mtx"
-        ));
+        ));
+
+TEST(CXSparse, Symmetric) {
+    cs_dl *sym = read_mtx(kTestMatrixDir + "symmetry/coordinate_symmetric_row.mtx");
+    cs_dl *expected = read_mtx(kTestMatrixDir + "symmetry/coordinate_symmetric_row_general.mtx");
+    EXPECT_EQ(6, sym->nzmax);
+    EXPECT_EQ(5, expected->nzmax);
+}

tests/eigen_test.cpp

@@ -150,4 +150,19 @@ INSTANTIATE_TEST_SUITE_P(
                           , Param{Param::Load_FMM_Vec, "vector_array.mtx"}
                           , Param{Param::Load_FMM_Vec, "vector_coordinate.mtx"}
                           , Param{Param::Load_FMM, "eye3_pattern.mtx"}
-                          ));
+                          ));
+
+
+TEST(EigenTest, Symmetric) {
+    SpRowMajor sym, expected;
+
+    {
+        std::ifstream f(kTestMatrixDir + "symmetry/coordinate_symmetric_row.mtx");
+        fast_matrix_market::read_matrix_market_eigen(f, sym);
+    }
+    {
+        std::ifstream f(kTestMatrixDir + "symmetry/coordinate_symmetric_row_general.mtx");
+        fast_matrix_market::read_matrix_market_eigen(f, expected);
+    }
+    EXPECT_TRUE(expected.isApprox(sym, 1e-6));
+}

tests/graphblas_test.cpp

@@ -534,6 +534,20 @@ TEST(GraphBLASMatrixTest, TypeComment) {
    }
 }
 
+TEST(GraphBLASMatrixTest, Symmetric) {
+    GrB_Matrix sym, expected;
+
+    {
+        std::ifstream f(kTestMatrixDir + "symmetry/coordinate_symmetric_row.mtx");
+        fast_matrix_market::read_matrix_market_graphblas(f, &sym);
+    }
+    {
+        std::ifstream f(kTestMatrixDir + "symmetry/coordinate_symmetric_row_general.mtx");
+        fast_matrix_market::read_matrix_market_graphblas(f, &expected);
+    }
+    EXPECT_TRUE(is_equal(expected, sym));
+}
+
 
 /////////////////////////////////////////////////////////
 ///// Vectors