Merge pull request #61 from HaoZeke/updFMMlib

MAINT: Update `fast_matrix_market` to 1.7.4

inst/include/fast_matrix_market/app/array.hpp

@@ -105,7 +105,9 @@ namespace fast_matrix_market {
         header.nnz = values.size();
 
         header.object = matrix;
-        header.field = get_field_type((const VT*)nullptr);
+        if (options.fill_header_field_type) {
+            header.field = get_field_type((const VT *) nullptr);
+        }
         header.format = array;
         header.symmetry = general;
 
@@ -121,4 +123,4 @@ namespace fast_matrix_market {
 #undef array_read_vector
 #undef array_write_vector
 #endif
-}
+}

inst/include/fast_matrix_market/app/doublet.hpp

@@ -89,7 +89,7 @@ namespace fast_matrix_market {
         header.object = vector;
         if (header.nnz > 0 && (values.cbegin() == values.cend())) {
             header.field = pattern;
-        } else if (header.field != pattern) {
+        } else if (header.field != pattern && options.fill_header_field_type) {
             header.field = get_field_type((const VT *) nullptr);
         }
         header.format = coordinate;
@@ -109,4 +109,4 @@ namespace fast_matrix_market {
 #undef doublet_read_vector
 #undef doublet_write_vector
 #endif
-}
+}

inst/include/fast_matrix_market/app/triplet.hpp

@@ -36,6 +36,75 @@ 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;
+        }
+
+        auto nnz = get_storage_nnz(header, options);
+        rows.resize(nnz);
+        cols.resize(nnz);
+        values.resize(nnz);
+
+        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 +113,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);
     }
 
     /**
@@ -88,7 +151,7 @@ namespace fast_matrix_market {
         header.object = matrix;
         if (header.nnz > 0 && (values.cbegin() == values.cend())) {
             header.field = pattern;
-        } else if (header.field != pattern) {
+        } else if (header.field != pattern && options.fill_header_field_type) {
             header.field = get_field_type((const VT *) nullptr);
         }
         header.format = coordinate;
@@ -122,7 +185,7 @@ namespace fast_matrix_market {
         header.object = matrix;
         if (header.nnz > 0 && (values.cbegin() == values.cend())) {
             header.field = pattern;
-        } else if (header.field != pattern) {
+        } else if (header.field != pattern && options.fill_header_field_type) {
             header.field = get_field_type((const VT *) nullptr);
         }
         header.format = coordinate;
@@ -143,4 +206,4 @@ namespace fast_matrix_market {
 #undef triplet_read_vector
 #undef triplet_write_vector
 #endif
-}
+}

inst/include/fast_matrix_market/app/user_type_string.hpp

@@ -0,0 +1,86 @@
+// Copyright (C) 2022-2023 Adam Lugowski. All rights reserved.
+// Use of this source code is governed by the BSD 2-clause license found in the LICENSE.txt file.
+// SPDX-License-Identifier: BSD-2-Clause
+
+/*
+ * Support std::string as a value type.
+ *
+ * Supports such things as a triplet where the value type is std::vector<std::string>, and such structure can read
+ * any Matrix Market file regardless of field type.
+ *
+ * IMPORTANT!
+ *
+ * If using this as a template for your own user type, note the weird #include ordering requirements!
+ *
+ * ```
+ * #include <fast_matrix_market/types.hpp>
+ *
+ * namespace fast_matrix_market {
+ *   // declare your type stuff here
+ * }
+ *
+ * // Now include the main header.
+ * #include <fast_matrix_market/fast_matrix_market.hpp>
+ * ```
+ */
+
+namespace fast_matrix_market {
+    /**
+     * (Needed for read) Parse a value.
+     *
+     * This method must find the end of the value, too. Likely the next newline or end of string is the end.
+     *
+     * @param pos starting character
+     * @param end end of string. Do not dereference >end.
+     * @param out out parameter of parsed value
+     * @return a pointer to the next character following the value. Likely the newline.
+     */
+    inline const char *read_value(const char *pos, const char *end, std::string &out, [[maybe_unused]] const read_options& options = {}) {
+        const char *field_start = pos;
+        // find the end of the line
+        while (pos != end && *pos != '\n') {
+            ++pos;
+        }
+        out = std::string(field_start, (pos - field_start));
+
+        return pos;
+    }
+
+    /**
+     * (Optional, used for read) Declare that read_value(std::string) declared above can read complex values too.
+     */
+    template<> struct can_read_complex<std::string> : std::true_type {};
+
+    /**
+     * (Needed for read) Used to handle skew-symmetric files. Must be defined regardless.
+     */
+    inline std::string negate(const std::string& o) {
+        return "-" + o;
+    }
+
+    /**
+     * (Needed for read) Default value to set for patterns (if option selected). Must be defined regardless.
+     */
+    inline std::string pattern_default_value([[maybe_unused]] const std::string* type) {
+        return "";
+    }
+
+    // If using default dense array loader, type must also work with std::plus<>.
+
+    /**
+     * (Needed for write) Used to determine what field type to set in the MatrixMarket header.
+     *
+     * IMPORTANT! This is the default type. If it does not match your actual data, set the appropriate type
+     * in the header, then use write_options::fill_header_field_type = false.
+     */
+    inline field_type get_field_type([[maybe_unused]] const std::string* type) {
+        return real;
+    }
+
+    /**
+     * (Needed for write) Write a value.
+     */
+    inline std::string value_to_string(const std::string& value, [[maybe_unused]] int precision) {
+        return value;
+    }
+}

inst/include/fast_matrix_market/chunking.hpp

@@ -57,7 +57,7 @@ namespace fast_matrix_market {
 
     template <typename ITER>
     bool is_all_spaces(ITER begin, ITER end) {
-        return std::all_of(begin, end, [](char c) { return c == ' ' || c == '\t'; });
+        return std::all_of(begin, end, [](char c) { return c == ' ' || c == '\t' || c == '\r'; });
     }
 
     /**

inst/include/fast_matrix_market/fast_matrix_market.hpp

@@ -8,30 +8,30 @@
 #pragma once
 
 #include <algorithm>
-#include <complex>
 #include <future>
 #include <iostream>
 #include <map>
 #include <string>
 #include <sstream>
 #include <utility>
+#include <vector>
 
 #include "types.hpp"
 
+// Support std::string as a user type
+#include "app/user_type_string.hpp"
+
 namespace fast_matrix_market {
 
     // Version macros.
     // Keep in sync with python/pyproject.toml
 #define FAST_MATRIX_MARKET_VERSION_MAJOR 1
 #define FAST_MATRIX_MARKET_VERSION_MINOR 7
-#define FAST_MATRIX_MARKET_VERSION_PATCH 2
+#define FAST_MATRIX_MARKET_VERSION_PATCH 4
 
     constexpr std::string_view kSpace = " ";
     constexpr std::string_view kNewline = "\n";
 
-    template<class T> struct is_complex : std::false_type {};
-    template<class T> struct is_complex<std::complex<T>> : std::true_type {};
-
     /**
      *
      */
@@ -120,6 +120,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);