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              :   ///
      68              :   /// Both source and destination are allowed to have a stride.
      69              :   template <typename T>
      70            1 :   void increment_multiply_on_right(gsl::not_null<T*> a, size_t a_offset,
      71              :                                    size_t a_stride, const T& b, size_t b_offset,
      72              :                                    size_t b_stride) const;
      73              : 
      74              :   /// Number of nonzero elements.
      75            1 :   size_t size() const { return matrix_elements_.size(); }
      76              : 
      77              :   /// Obtains element at (destindex,srcindex).
      78              :   /// Slow, should not use very often in critical situations.
      79            1 :   double operator()(size_t row_dest_index, size_t column_src_index) const;
      80              : 
      81              :  private:
      82              :   // Holds all nonzero matrix elements.
      83            0 :   std::vector<double> matrix_elements_;
      84              :   // Holds indices of all rows of the matrix that have nonzero matrix elements.
      85            0 :   std::vector<size_t> row_dest_indices_;
      86              :   // Holds the number of columns in every row that has nonzero matrix elements.
      87              :   // row_dest_indices_ and num_columns_per_row_ have the same size (which is the
      88              :   // number of rows with nonzero matrix elements).
      89              :   // The sum of all the num_columns_per_row_ is the size of matrix_elements_.
      90            0 :   std::vector<size_t> num_columns_per_row_;
      91              :   // Holds the indices of all columns that have nonzero matrix elements.
      92              :   // column_indices_ has the same size as matrix_elements_, so many
      93              :   // of the column_src_indices_ might be repeated.
      94            0 :   std::vector<size_t> column_src_indices_;
      95              : };