Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <blaze/math/CompressedMatrix.h>
       7             : #include <cstddef>
       8             : #include <unordered_map>
       9             : #include <utility>
      10             : #include <vector>
      11             : 
      12             : #include "DataStructures/SimpleSparseMatrix.hpp"
      13             : #include "Utilities/ErrorHandling/Assert.hpp"
      14             : #include "Utilities/Gsl.hpp"
      15             : 
      16             : namespace sparse_matrix_filler_detail {
      17             : /// Hash on a pair of size_ts
      18             : struct PairHasher {
      19             :   size_t operator()(const std::pair<size_t, size_t>& key) const {
      20             :     ASSERT(key.first < 4294967295, "Key.first too large: " << key.first);
      21             :     ASSERT(key.second < 4294967295, "Key.second too large: " << key.second);
      22             :     return key.first bitor (key.second << 32);
      23             :   }
      24             : };
      25             : }  // namespace sparse_matrix_filler_detail
      26             : 
      27             : /*!
      28             :  * \ingroup DataStructuresGroup
      29             :  * \brief A data structure used to fill sparse matrices.
      30             :  *
      31             :  * We normally use blaze::CompressedMatrix to hold sparse matrices and to
      32             :  * operate on them.
      33             :  *
      34             :  * Unfortunately to fill a blaze::CompressedMatrix efficiently you
      35             :  * need to fill nonzero matrix elements in a well-defined order
      36             :  * (i.e. fill all nonzero elements in the top row in the order of
      37             :  * nonzero columns, and repeat for each other row, in order) using
      38             :  * CompressedMatrix's reserve(), append(), and finalize() functions.
      39             :  *
      40             :  * However, we are interested in a use case where:
      41             :  *
      42             :  *   1. We are filling the matrix elements using a complicated
      43             :  *      algorithm that computes (as opposed to systematically loops
      44             :  *      over) the matrix element indices to be filled, so we are
      45             :  *      effectively filling matrix elements in a random order.
      46             :  *
      47             :  *   2. The algorithm can traverse the same matrix element hundreds
      48             :  *      of times, meaning that the first time an element is traversed
      49             :  *      it should be filled, but subsequent times it is traversed it
      50             :  *      should be incremented.
      51             :  *
      52             :  * SparseMatrixFiller provides a function to fill matrix elements, keeping
      53             :  * in mind the above two issues.  And it provides a function to copy its
      54             :  * sparse matrix into blaze::CompressedMatrix efficiently.
      55             :  *
      56             :  * We consider two ways to keep track of unique elements:
      57             :  *
      58             :  * 1. Store matrix elements in a full vector that is indexed by
      59             :  *    matrix_elements[src_index+max_src_index*dest_index].  Set this
      60             :  *    vector to zero initially, and increment it each time 'add' is called.
      61             :  *    This is simple, but has a major issue: memory.
      62             :  *    In SpEC, I (Mark Scheel) have measured that computing TensorYlm
      63             :  *    filtering matrices using this method, after l_max=40 , the
      64             :  *    slope of CPU-h vs l_max for this method has a kink, and CPU-h
      65             :  *    increases like l_max^2 instead of l_max.  Presumably this is
      66             :  *    because of cache misses when accessing the very large
      67             :  *    matrix_elements vector at random locations.
      68             :  *
      69             :  * 2. Have a std::unordered_map<std::pair<size_t,size_t>, size_t>>
      70             :  *    called `element_index` that keeps track of the unique elements.
      71             :  *    Also have three vectors `dest_indices`, `src_indices`,
      72             :  *    and `matrix_elements` that keep the nonzero elements.
      73             :  *    element_index is indexed by
      74             :  *    element_index[{dest_index,src_index}] = index-into-matrix-elements
      75             :  *    So each time we add a new element, we check element_index to see
      76             :  *    if that element has not already been added, and if not, we push_back
      77             :  *    the three vectors.  If it has already been added, we simply do
      78             :  *    matrix_elements[index-into-matrix-elements] += new_element.
      79             :  *
      80             :  * I (Mark Scheel) have found in SpEC for computing TensorYlm filtering
      81             :  * matrices that the CPU cost of algorithm 2. above scales linearly
      82             :  * with l_max at least up to l_max=140, but for l_max < 40 or so,
      83             :  * algorithm 2. is about a factor of 3 slower than algorithm 1.
      84             :  * So SpEC uses both algorithms, depending on the value of l_max.
      85             :  *
      86             :  * For SpECTRE, we plan to do the same: for small l_max we will use
      87             :  * algorithm 1 and for large l_max we will use algorithm 2.  The
      88             :  * cutoff for l_max will depend on other factors (like the rank of the
      89             :  * Tensor being filtered) and will be determined by whoever calls
      90             :  * SparseMatrixFiller.
      91             :  */
      92           1 : class SparseMatrixFiller {
      93             :  public:
      94             :   /// Assume num_rows and num_cols are equal.
      95             :   /// If use_map_method is true, uses algorithm 2.
      96             :   /// Scale is the overall scale of the matrix; scale is used to remove
      97             :   /// elements that are smaller than roundoff.
      98           1 :   SparseMatrixFiller(size_t num_cols, bool use_map_method, double scale);
      99             : 
     100             :   /// Add (or increment) an element such that for v = M x, where M is
     101             :   /// the matrix and v and x are vectors, we compute the contribution
     102             :   /// to v[dest_index] from x[src_index].
     103           1 :   void add(double element, size_t dest_index, size_t src_index);
     104             : 
     105             :   /// Fills the matrix.
     106             :   /// Note that the order in which blaze matrices must be filled depends
     107             :   /// on whether the matrix is rowMajor or columnMajor, so changing
     108             :   /// the template parameter here will fail without changing the code.
     109           1 :   void fill(gsl::not_null<blaze::CompressedMatrix<double, blaze::rowMajor>*>
     110             :                 matrix) const;
     111             : 
     112             :   /// Fills a SimpleSparseMatrix.
     113           1 :   void fill(gsl::not_null<SimpleSparseMatrix*> matrix) const;
     114             : 
     115             :  private:
     116           0 :   void fill_sparse_matrix_elements(
     117             :       gsl::not_null<std::vector<SparseMatrixElement>*> data) const;
     118           0 :   size_t num_rows_{0}, num_cols_{0};
     119           0 :   bool use_map_method_{true};
     120           0 :   double scale_{1.0};
     121           0 :   std::vector<double> matrix_elements_{};
     122             :   std::unordered_map<std::pair<size_t, size_t>, size_t,
     123             :                      sparse_matrix_filler_detail::PairHasher>
     124           0 :       element_index_{};
     125           0 :   std::vector<size_t> dest_indices_{};
     126           0 :   std::vector<size_t> src_indices_{};
     127             : };