Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <vector>
       8             : 
       9             : #include "Utilities/Gsl.hpp"
      10             : 
      11             : /// Struct used to fill SimpleSparsematrix
      12           1 : struct SparseMatrixElement {
      13           0 :   SparseMatrixElement(const size_t row, const size_t column, const double val)
      14             :       : row_dest(row), column_src(column), value(val) {}
      15           0 :   size_t row_dest, column_src;
      16           0 :   double value;
      17             : };
      18             : 
      19             : /*!
      20             :  * \brief A simple fast sparse matrix.
      21             :  *
      22             :  * Supports only a limited number of operations to keep it simple:
      23             :  * filling the matrix, accessing elements, and multiplication by
      24             :  * a vector on the right.
      25             :  *
      26             :  * We use SimpleSparseMatrix because filling a Blaze sparse matrix is
      27             :  * enormously slow. For example, here are timings for creating a
      28             :  * rank-one TensorYlmFilter sparse matrix using SimpleSparseMatrix
      29             :  * (column labeled "this") and using blaze::CompressedMatrix (column
      30             :  * labeled "Blaze").  The column "Blaze-wo-blaze" is the same as the
      31             :  * "Blaze" column (which includes timing of all the TensorYlmFilter
      32             :  * code that doesn't explicitly depend on blaze) minus the timings of
      33             :  * the blaze calls that 1. resize the matrix, 2. reserve elements, and
      34             :  * 3. fill the blaze matrix.  The timings for those individual blaze
      35             :  * calls are in subsequent columns.  Note that "Blaze-wo-blaze" is
      36             :  * roughly equivalent to "this" in speed.  Also note that the "this"
      37             :  * timings were done in SpEC, which uses SimpleSparseMatrix [which was
      38             :  * ported from SpEC to SpECTRE and simplified here].
      39             :  *
      40             :  * l_max this Blaze Blaze-wo-blaze blaze.resize blaze.reserve blaze.fill
      41             :  * ----------------------------------------------------------------------------
      42             :  *  8    0ms     5ms     1ms              1ms           2ms          0ms
      43             :  * 18    1ms   101ms     2ms             39ms          51ms          8ms
      44             :  * 28    3ms   479ms     3ms            181ms         246ms         48ms
      45             :  * 38    4ms  1555ms     4ms            571ms         822ms        157ms
      46             :  * 48    5ms  3668ms     6ms           1395ms        1873ms        393ms
      47             :  * 58    6ms  7604ms     6ms           2908ms        3855ms        832ms
      48             :  * 68    7ms 14176ms     8ms           5416ms        7214ms       1536ms
      49             :  * 78    8ms 24234ms     8ms           9288ms       12304ms       2632ms
      50             :  *
      51             :  */
      52           1 : class SimpleSparseMatrix {
      53             :  public:
      54           0 :   SimpleSparseMatrix() = default;
      55             : 
      56             :   /// Fills with data.
      57             :   /// Normally called only by SparseMatrixFiller.
      58             :   /// ASSUMES that data has already been sorted by row, then column.
      59           1 :   void fill(const std::vector<SparseMatrixElement>& data);
      60             : 
      61             :   /// Does a += m*b where m is the sparse matrix.
      62             :   /// T is a container of doubles, like a std::vector<double>.
      63             :   /// Assumes b points into contiguous storage that allows one
      64             :   /// to access b[i] for b_offset <= i <= b_offset+max(column_indices).
      65             :   /// Assumes a points into contiguous storage that allows one
      66             :   /// to access a[i] for a_offset <= a <= a_offset+max(row_indices).
      67             :   template <typename T>
      68           1 :   void increment_multiply_on_right(gsl::not_null<T*> a, size_t a_offset,
      69             :                                    const T& b, size_t b_offset) const;
      70             : 
      71             :   /// Number of nonzero elements.
      72           1 :   size_t size() const { return matrix_elements_.size(); }
      73             : 
      74             :   /// Obtains element at (destindex,srcindex).
      75             :   /// Slow, should not use very often in critical situations.
      76           1 :   double operator()(size_t row_dest_index, size_t column_src_index) const;
      77             : 
      78             :  private:
      79             :   // Holds all nonzero matrix elements.
      80           0 :   std::vector<double> matrix_elements_;
      81             :   // Holds indices of all rows of the matrix that have nonzero matrix elements.
      82           0 :   std::vector<size_t> row_dest_indices_;
      83             :   // Holds the number of columns in every row that has nonzero matrix elements.
      84             :   // row_dest_indices_ and num_columns_per_row_ have the same size (which is the
      85             :   // number of rows with nonzero matrix elements).
      86             :   // The sum of all the num_columns_per_row_ is the size of matrix_elements_.
      87           0 :   std::vector<size_t> num_columns_per_row_;
      88             :   // Holds the indices of all columns that have nonzero matrix elements.
      89             :   // column_indices_ has the same size as matrix_elements_, so many
      90             :   // of the column_src_indices_ might be repeated.
      91           0 :   std::vector<size_t> column_src_indices_;
      92             : };