MFEM  v4.5.2
Finite element discretization library
sparsemat.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_SPARSEMAT_HPP
13 #define MFEM_SPARSEMAT_HPP
14 
15 // Data types for sparse matrix
16 
17 #include "../general/backends.hpp"
18 #include "../general/mem_alloc.hpp"
19 #include "../general/mem_manager.hpp"
20 #include "../general/device.hpp"
21 #include "../general/table.hpp"
22 #include "../general/globals.hpp"
23 #include "densemat.hpp"
24 
25 #if defined(MFEM_USE_HIP)
26 #if (HIP_VERSION_MAJOR * 100 + HIP_VERSION_MINOR) < 502
27 #include <hipsparse.h>
28 #else
29 #include <hipsparse/hipsparse.h>
30 #endif
31 #endif
32 
33 
34 namespace mfem
35 {
36 
37 class
38 #if defined(__alignas_is_defined)
39  alignas(double)
40 #endif
41  RowNode
42 {
43 public:
44  double Value;
45  RowNode *Prev;
46  int Column;
47 };
48 
49 /// Data type sparse matrix
51 {
52 protected:
53  /// @name Arrays used by the CSR storage format.
54  /** */
55  ///@{
56  /// @brief %Array with size (#height+1) containing the row offsets.
57  /** The data for row r, 0 <= r < height, is at offsets j, I[r] <= j < I[r+1].
58  The offsets, j, are indices in the #J and #A arrays. The first entry in
59  this array is always zero, I[0] = 0, and the last entry, I[height], gives
60  the total number of entries stored (at a minimum, all nonzeros must be
61  represented) in the sparse matrix. */
63  /** @brief %Array with size #I[#height], containing the column indices for
64  all matrix entries, as indexed by the #I array. */
66  /** @brief %Array with size #I[#height], containing the actual entries of the
67  sparse matrix, as indexed by the #I array. */
69  ///@}
70 
71  /** @brief %Array of linked lists, one for every row. This array represents
72  the linked list (LIL) storage format. */
73  RowNode **Rows;
74 
75  mutable int current_row;
76  mutable int* ColPtrJ;
77  mutable RowNode ** ColPtrNode;
78 
79  /// Transpose of A. Owned. Used to perform MultTranspose() on devices.
80  mutable SparseMatrix *At;
81 
82 #ifdef MFEM_USE_MEMALLOC
85 #endif
86 
87  /// Are the columns sorted already.
88  bool isSorted;
89 
90  void Destroy(); // Delete all owned data
91  void SetEmpty(); // Init all entries with empty values
92 
93  bool useGPUSparse = true; // Use cuSPARSE or hipSPARSE if available
94 
95  // Initialize cuSPARSE/hipSPARSE
96  void InitGPUSparse();
97 
98 #ifdef MFEM_USE_CUDA_OR_HIP
99  // common for hipSPARSE and cuSPARSE
100  static int SparseMatrixCount;
101  static size_t bufferSize;
102  static void *dBuffer;
103  mutable bool initBuffers = false;
104 
105 #if defined(MFEM_USE_CUDA)
106  cusparseStatus_t status;
107  static cusparseHandle_t handle;
108  cusparseMatDescr_t descr = 0;
109 
110 #if CUDA_VERSION >= 10010
111  mutable cusparseSpMatDescr_t matA_descr;
112  mutable cusparseDnVecDescr_t vecX_descr;
113  mutable cusparseDnVecDescr_t vecY_descr;
114 #else // CUDA_VERSION >= 10010
115  mutable cusparseMatDescr_t matA_descr;
116 #endif // CUDA_VERSION >= 10010
117 
118 #else // defined(MFEM_USE_CUDA)
119  hipsparseStatus_t status;
120  static hipsparseHandle_t handle;
121  hipsparseMatDescr_t descr = 0;
122 
123  mutable hipsparseSpMatDescr_t matA_descr;
124  mutable hipsparseDnVecDescr_t vecX_descr;
125  mutable hipsparseDnVecDescr_t vecY_descr;
126 #endif // defined(MFEM_USE_CUDA)
127 #endif // MFEM_USE_CUDA_OR_HIP
128 
129 public:
130  /// Create an empty SparseMatrix.
132  {
133  SetEmpty();
134 
135  InitGPUSparse();
136  }
137 
138  /** @brief Create a sparse matrix with flexible sparsity structure using a
139  row-wise linked list (LIL) format. */
140  /** New entries are added as needed by methods like AddSubMatrix(),
141  SetSubMatrix(), etc. Calling Finalize() will convert the SparseMatrix to
142  the more compact compressed sparse row (CSR) format. */
143  explicit SparseMatrix(int nrows, int ncols = -1);
144 
145  /** @brief Create a sparse matrix in CSR format. Ownership of @a i, @a j, and
146  @a data is transferred to the SparseMatrix. */
147  SparseMatrix(int *i, int *j, double *data, int m, int n);
148 
149  /** @brief Create a sparse matrix in CSR format. Ownership of @a i, @a j, and
150  @a data is optionally transferred to the SparseMatrix. */
151  /** If the parameter @a data is NULL, then the internal #A array is allocated
152  by this constructor (initializing it with zeros and taking ownership,
153  regardless of the parameter @a owna). */
154  SparseMatrix(int *i, int *j, double *data, int m, int n, bool ownij,
155  bool owna, bool issorted);
156 
157  /** @brief Create a sparse matrix in CSR format where each row has space
158  allocated for exactly @a rowsize entries. */
159  /** SetRow() can then be called or the #I, #J, #A arrays can be used
160  directly. */
161  SparseMatrix(int nrows, int ncols, int rowsize);
162 
163  /// Copy constructor (deep copy).
164  /** If @a mat is finalized and @a copy_graph is false, the #I and #J arrays
165  will use a shallow copy (copy the pointers only) without transferring
166  ownership.
167  If @a mt is MemoryType::PRESERVE the memory type of the resulting
168  SparseMatrix's #I, #J, and #A arrays will be the same as @a mat,
169  otherwise the type will be @a mt for those arrays that are deep
170  copied. */
171  SparseMatrix(const SparseMatrix &mat, bool copy_graph = true,
173 
174  /// Create a SparseMatrix with diagonal @a v, i.e. A = Diag(v)
175  SparseMatrix(const Vector & v);
176 
177  /// @brief Sets the height and width of the matrix.
178  /** @warning This does not modify in any way the underlying CSR or LIL
179  representation of the matrix.
180 
181  This function should generally be called when manually constructing the
182  CSR #I, #J, and #A arrays in conjunction with the
183  SparseMatrix::SparseMatrix() constructor. */
184  void OverrideSize(int height_, int width_);
185 
186  /** @brief Runtime option to use cuSPARSE or hipSPARSE. Only valid when using
187  a CUDA or HIP backend.
188 
189  @note This option is enabled by default, so typically one would use this
190  method to disable the use of cuSPARSE/hipSPARSE. */
191  void UseGPUSparse(bool useGPUSparse_ = true) { useGPUSparse = useGPUSparse_;}
192  /// Deprecated equivalent of UseGPUSparse().
193  MFEM_DEPRECATED
194  void UseCuSparse(bool useCuSparse_ = true) { UseGPUSparse(useCuSparse_); }
195 
196  /// Assignment operator: deep copy
197  SparseMatrix& operator=(const SparseMatrix &rhs);
198 
199  /** @brief Clear the contents of the SparseMatrix and make it a reference to
200  @a master */
201  /** After this call, the matrix will point to the same data as @a master but
202  it will not own its data. The @a master must be finalized. */
203  void MakeRef(const SparseMatrix &master);
204 
205  /// For backward compatibility, define Size() to be synonym of Height().
206  int Size() const { return Height(); }
207 
208  /// Clear the contents of the SparseMatrix.
209  void Clear() { Destroy(); SetEmpty(); }
210 
211  /** @brief Clear the cuSPARSE/hipSPARSE descriptors.
212  This must be called after releasing the device memory of A. */
213  void ClearGPUSparse();
214  /// Deprecated equivalent of ClearGPUSparse().
215  MFEM_DEPRECATED
217 
218  /// Check if the SparseMatrix is empty.
219  bool Empty() const { return (A == NULL) && (Rows == NULL); }
220 
221  /// Return the array #I.
222  inline int *GetI() { return I; }
223  /// Return the array #I, const version.
224  inline const int *GetI() const { return I; }
225 
226  /// Return the array #J.
227  inline int *GetJ() { return J; }
228  /// Return the array #J, const version.
229  inline const int *GetJ() const { return J; }
230 
231  /// Return the element data, i.e. the array #A.
232  inline double *GetData() { return A; }
233  /// Return the element data, i.e. the array #A, const version.
234  inline const double *GetData() const { return A; }
235 
236  // Memory access methods for the #I array.
237  Memory<int> &GetMemoryI() { return I; }
238  const Memory<int> &GetMemoryI() const { return I; }
239  const int *ReadI(bool on_dev = true) const
240  { return mfem::Read(I, Height()+1, on_dev); }
241  int *WriteI(bool on_dev = true)
242  { return mfem::Write(I, Height()+1, on_dev); }
243  int *ReadWriteI(bool on_dev = true)
244  { return mfem::ReadWrite(I, Height()+1, on_dev); }
245  const int *HostReadI() const
246  { return mfem::Read(I, Height()+1, false); }
247  int *HostWriteI()
248  { return mfem::Write(I, Height()+1, false); }
250  { return mfem::ReadWrite(I, Height()+1, false); }
251 
252  // Memory access methods for the #J array.
253  Memory<int> &GetMemoryJ() { return J; }
254  const Memory<int> &GetMemoryJ() const { return J; }
255  const int *ReadJ(bool on_dev = true) const
256  { return mfem::Read(J, J.Capacity(), on_dev); }
257  int *WriteJ(bool on_dev = true)
258  { return mfem::Write(J, J.Capacity(), on_dev); }
259  int *ReadWriteJ(bool on_dev = true)
260  { return mfem::ReadWrite(J, J.Capacity(), on_dev); }
261  const int *HostReadJ() const
262  { return mfem::Read(J, J.Capacity(), false); }
263  int *HostWriteJ()
264  { return mfem::Write(J, J.Capacity(), false); }
266  { return mfem::ReadWrite(J, J.Capacity(), false); }
267 
268  // Memory access methods for the #A array.
270  const Memory<double> &GetMemoryData() const { return A; }
271  const double *ReadData(bool on_dev = true) const
272  { return mfem::Read(A, A.Capacity(), on_dev); }
273  double *WriteData(bool on_dev = true)
274  { return mfem::Write(A, A.Capacity(), on_dev); }
275  double *ReadWriteData(bool on_dev = true)
276  { return mfem::ReadWrite(A, A.Capacity(), on_dev); }
277  const double *HostReadData() const
278  { return mfem::Read(A, A.Capacity(), false); }
279  double *HostWriteData()
280  { return mfem::Write(A, A.Capacity(), false); }
282  { return mfem::ReadWrite(A, A.Capacity(), false); }
283 
284  /// Returns the number of elements in row @a i.
285  int RowSize(const int i) const;
286 
287  /// Returns the maximum number of elements among all rows.
288  int MaxRowSize() const;
289 
290  /// Return a pointer to the column indices in a row.
291  int *GetRowColumns(const int row);
292  /// Return a pointer to the column indices in a row, const version.
293  const int *GetRowColumns(const int row) const;
294 
295  /// Return a pointer to the entries in a row.
296  double *GetRowEntries(const int row);
297  /// Return a pointer to the entries in a row, const version.
298  const double *GetRowEntries(const int row) const;
299 
300  /// Change the width of a SparseMatrix.
301  /*!
302  * If width_ = -1 (DEFAULT), this routine will set the new width
303  * to the actual Width of the matrix awidth = max(J) + 1.
304  * Values 0 <= width_ < awidth are not allowed (error check in Debug Mode only)
305  *
306  * This method can be called for matrices finalized or not.
307  */
308  void SetWidth(int width_ = -1);
309 
310  /// Returns the actual Width of the matrix.
311  /*! This method can be called for matrices finalized or not. */
312  int ActualWidth() const;
313 
314  /// Sort the column indices corresponding to each row.
315  void SortColumnIndices();
316 
317  /** @brief Move the diagonal entry to the first position in each row,
318  preserving the order of the rest of the columns. */
319  void MoveDiagonalFirst();
320 
321  /// Returns reference to a_{ij}.
322  virtual double &Elem(int i, int j);
323 
324  /// Returns constant reference to a_{ij}.
325  virtual const double &Elem(int i, int j) const;
326 
327  /// Returns reference to A[i][j].
328  double &operator()(int i, int j);
329 
330  /// Returns reference to A[i][j].
331  const double &operator()(int i, int j) const;
332 
333  /// Returns the Diagonal of A
334  void GetDiag(Vector & d) const;
335 
336  /// Produces a DenseMatrix from a SparseMatrix
337  DenseMatrix *ToDenseMatrix() const;
338 
339  /// Produces a DenseMatrix from a SparseMatrix
340  void ToDenseMatrix(DenseMatrix & B) const;
341 
342  virtual MemoryClass GetMemoryClass() const
343  {
344  return Finalized() ?
346  }
347 
348  /// Matrix vector multiplication.
349  virtual void Mult(const Vector &x, Vector &y) const;
350 
351  /// y += A * x (default) or y += a * A * x
352  virtual void AddMult(const Vector &x, Vector &y,
353  const double a = 1.0) const;
354 
355  /// Multiply a vector with the transposed matrix. y = At * x
356  /** If the matrix is modified, call ResetTranspose() and optionally
357  EnsureMultTranspose() to make sure this method uses the correct updated
358  transpose. */
359  virtual void MultTranspose(const Vector &x, Vector &y) const;
360 
361  /// y += At * x (default) or y += a * At * x
362  /** If the matrix is modified, call ResetTranspose() and optionally
363  EnsureMultTranspose() to make sure this method uses the correct updated
364  transpose. */
365  virtual void AddMultTranspose(const Vector &x, Vector &y,
366  const double a = 1.0) const;
367 
368  /** @brief Build and store internally the transpose of this matrix which will
369  be used in the methods AddMultTranspose(), MultTranspose(), and
370  AbsMultTranspose(). */
371  /** If this method has been called, the internal transpose matrix will be
372  used to perform the action of the transpose matrix in AddMultTranspose(),
373  MultTranspose(), and AbsMultTranspose().
374 
375  Warning: any changes in this matrix will invalidate the internal
376  transpose. To rebuild the transpose, call ResetTranspose() followed by
377  (optionally) a call to this method. If the internal transpose is already
378  built, this method has no effect.
379 
380  When any non-serial-CPU backend is enabled, i.e. the call
381  Device::Allows(~ Backend::CPU_MASK) returns true, the above methods
382  require the internal transpose to be built. If that is not the case (i.e.
383  the internal transpose is not built), these methods will automatically
384  call EnsureMultTranspose(). When using any backend from
385  Backend::CPU_MASK, calling this method is optional.
386 
387  This method can only be used when the sparse matrix is finalized.
388 
389  @sa EnsureMultTranspose(), ResetTranspose(). */
390  void BuildTranspose() const;
391 
392  /** Reset (destroy) the internal transpose matrix. See BuildTranspose() for
393  more details. */
394  void ResetTranspose() const;
395 
396  /** @brief Ensures that the matrix is capable of performing MultTranspose(),
397  AddMultTranspose(), and AbsMultTranspose(). */
398  /** For non-serial-CPU backends (e.g. GPU, OpenMP), multiplying by the
399  transpose requires that the internal transpose matrix be already built.
400  When such a backend is enabled, this function will build the internal
401  transpose matrix, see BuildTranspose().
402 
403  For the serial CPU backends, the internal transpose is not required, and
404  this function is a no-op. This allows for significant memory savings
405  when the internal transpose matrix is not required. */
406  void EnsureMultTranspose() const;
407 
408  void PartMult(const Array<int> &rows, const Vector &x, Vector &y) const;
409  void PartAddMult(const Array<int> &rows, const Vector &x, Vector &y,
410  const double a=1.0) const;
411 
412  /// y = A * x, treating all entries as booleans (zero=false, nonzero=true).
413  /** The actual values stored in the data array, #A, are not used - this means
414  that all entries in the sparsity pattern are considered to be true by
415  this method. */
416  void BooleanMult(const Array<int> &x, Array<int> &y) const;
417 
418  /// y = At * x, treating all entries as booleans (zero=false, nonzero=true).
419  /** The actual values stored in the data array, #A, are not used - this means
420  that all entries in the sparsity pattern are considered to be true by
421  this method. */
422  void BooleanMultTranspose(const Array<int> &x, Array<int> &y) const;
423 
424  /// y = |A| * x, using entry-wise absolute values of matrix A
425  void AbsMult(const Vector &x, Vector &y) const;
426 
427  /// y = |At| * x, using entry-wise absolute values of the transpose of matrix A
428  /** If the matrix is modified, call ResetTranspose() and optionally
429  EnsureMultTranspose() to make sure this method uses the correct updated
430  transpose. */
431  void AbsMultTranspose(const Vector &x, Vector &y) const;
432 
433  /// Compute y^t A x
434  double InnerProduct(const Vector &x, const Vector &y) const;
435 
436  /// For all i compute \f$ x_i = \sum_j A_{ij} \f$
437  void GetRowSums(Vector &x) const;
438  /// For i = irow compute \f$ x_i = \sum_j | A_{i, j} | \f$
439  double GetRowNorml1(int irow) const;
440 
441  /// This virtual method is not supported: it always returns NULL.
442  virtual MatrixInverse *Inverse() const;
443 
444  /// Eliminates a column from the transpose matrix.
445  void EliminateRow(int row, const double sol, Vector &rhs);
446 
447  /// Eliminates a row from the matrix.
448  /*!
449  * - If @a dpolicy = #DIAG_ZERO, all the entries in the row will be set to 0.
450  * - If @a dpolicy = #DIAG_ONE (matrix must be square), the diagonal entry
451  * will be set equal to 1 and all other entries in the row to 0.
452  * - The policy #DIAG_KEEP is not supported.
453  */
454  void EliminateRow(int row, DiagonalPolicy dpolicy = DIAG_ZERO);
455 
456  /// Eliminates the column @a col from the matrix.
457  /** - If @a dpolicy = #DIAG_ZERO, all entries in the column will be set to 0.
458  - If @a dpolicy = #DIAG_ONE (matrix must be square), the diagonal entry
459  will be set equal to 1 and all other entries in the column to 0.
460  - The policy #DIAG_KEEP is not supported. */
461  void EliminateCol(int col, DiagonalPolicy dpolicy = DIAG_ZERO);
462 
463  /// Eliminate all columns i for which @a cols[i] != 0.
464  /** Elimination of a column means that all entries in the column are set to
465  zero. In addition, if the pointers @a x and @a b are not NULL, the
466  eliminated matrix entries are multiplied by the corresponding solution
467  value in @a *x and subtracted from the r.h.s. vector, @a *b. */
468  void EliminateCols(const Array<int> &cols, const Vector *x = NULL,
469  Vector *b = NULL);
470 
471  /** @brief Similar to EliminateCols + save the eliminated entries into
472  @a Ae so that (*this) + Ae is equal to the original matrix. */
473  void EliminateCols(const Array<int> &col_marker, SparseMatrix &Ae);
474 
475  /// Eliminate row @a rc and column @a rc and modify the @a rhs using @a sol.
476  /** Eliminates the column @a rc to the @a rhs, deletes the row @a rc and
477  replaces the element (rc,rc) with 1.0; assumes that element (i,rc)
478  is assembled if and only if the element (rc,i) is assembled.
479  By default, elements (rc,rc) are set to 1.0, although this behavior
480  can be adjusted by changing the @a dpolicy parameter. */
481  void EliminateRowCol(int rc, const double sol, Vector &rhs,
482  DiagonalPolicy dpolicy = DIAG_ONE);
483 
484  /** @brief Similar to
485  EliminateRowCol(int, const double, Vector &, DiagonalPolicy), but
486  multiple values for eliminated unknowns are accepted, and accordingly
487  multiple right-hand-sides are used. */
488  void EliminateRowColMultipleRHS(int rc, const Vector &sol,
489  DenseMatrix &rhs,
490  DiagonalPolicy dpolicy = DIAG_ONE);
491 
492  /// Perform elimination and set the diagonal entry to the given value
493  void EliminateRowColDiag(int rc, double value);
494 
495  /// Eliminate row @a rc and column @a rc.
496  void EliminateRowCol(int rc, DiagonalPolicy dpolicy = DIAG_ONE);
497 
498  /** @brief Similar to EliminateRowCol(int, DiagonalPolicy) + save the
499  eliminated entries into @a Ae so that (*this) + Ae is equal to the
500  original matrix */
501  void EliminateRowCol(int rc, SparseMatrix &Ae,
502  DiagonalPolicy dpolicy = DIAG_ONE);
503 
504  /** @brief Eliminate essential (Dirichlet) boundary conditions.
505 
506  @param[in] ess_dofs indices of the degrees of freedom belonging to the
507  essential boundary conditions.
508  @param[in] diag_policy policy for diagonal entries. */
509  void EliminateBC(const Array<int> &ess_dofs,
510  DiagonalPolicy diag_policy);
511 
512  /// If a row contains only one diag entry of zero, set it to 1.
513  void SetDiagIdentity();
514  /// If a row contains only zeros, set its diagonal to 1.
515  virtual void EliminateZeroRows(const double threshold = 1e-12);
516 
517  /// Gauss-Seidel forward and backward iterations over a vector x.
518  void Gauss_Seidel_forw(const Vector &x, Vector &y) const;
519  void Gauss_Seidel_back(const Vector &x, Vector &y) const;
520 
521  /// Determine appropriate scaling for Jacobi iteration
522  double GetJacobiScaling() const;
523  /** One scaled Jacobi iteration for the system A x = b.
524  x1 = x0 + sc D^{-1} (b - A x0) where D is the diag of A.
525  Absolute values of D are used when use_abs_diag = true. */
526  void Jacobi(const Vector &b, const Vector &x0, Vector &x1,
527  double sc, bool use_abs_diag = false) const;
528 
529  /// x = sc b / A_ii. When use_abs_diag = true, |A_ii| is used.
530  void DiagScale(const Vector &b, Vector &x,
531  double sc = 1.0, bool use_abs_diag = false) const;
532 
533  /** x1 = x0 + sc D^{-1} (b - A x0) where \f$ D_{ii} = \sum_j |A_{ij}| \f$. */
534  void Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
535  double sc = 1.0) const;
536 
537  /** x1 = x0 + sc D^{-1} (b - A x0) where \f$ D_{ii} = \sum_j A_{ij} \f$. */
538  void Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
539  double sc = 1.0) const;
540 
541  /** @brief Finalize the matrix initialization, switching the storage format
542  from LIL to CSR. */
543  /** This method should be called once, after the matrix has been initialized.
544  Internally, this method converts the matrix from row-wise linked list
545  (LIL) format into CSR (compressed sparse row) format. */
546  virtual void Finalize(int skip_zeros = 1) { Finalize(skip_zeros, false); }
547 
548  /// A slightly more general version of the Finalize(int) method.
549  void Finalize(int skip_zeros, bool fix_empty_rows);
550 
551  /// Returns whether or not CSR format has been finalized.
552  bool Finalized() const { return !A.Empty(); }
553  /// Returns whether or not the columns are sorted.
554  bool ColumnsAreSorted() const { return isSorted; }
555 
556  /** @brief Remove entries smaller in absolute value than a given tolerance
557  @a tol. If @a fix_empty_rows is true, a zero value is inserted in the
558  diagonal entry (for square matrices only) */
559  void Threshold(double tol, bool fix_empty_rows = false);
560 
561  /** Split the matrix into M x N blocks of sparse matrices in CSR format.
562  The 'blocks' array is M x N (i.e. M and N are determined by its
563  dimensions) and its entries are overwritten by the new blocks. */
564  void GetBlocks(Array2D<SparseMatrix *> &blocks) const;
565 
566  void GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
567  DenseMatrix &subm) const;
568 
569  /** @brief Initialize the SparseMatrix for fast access to the entries of the
570  given @a row which becomes the "current row". */
571  /** Fast access to the entries of the "current row" can be performed using
572  the methods: SearchRow(const int), _Add_(const int, const double),
573  _Set_(const int, const double), and _Get_(const int). */
574  inline void SetColPtr(const int row) const;
575  /** @brief Reset the "current row" set by calling SetColPtr(). This method
576  must be called between any two calls to SetColPtr(). */
577  inline void ClearColPtr() const;
578  /// Perform a fast search for an entry in the "current row". See SetColPtr().
579  /** If the matrix is not finalized and the entry is not found in the
580  SparseMatrix, it will be added to the sparsity pattern initialized with
581  zero. If the matrix is finalized and the entry is not found, an error
582  will be generated. */
583  inline double &SearchRow(const int col);
584  /// Add a value to an entry in the "current row". See SetColPtr().
585  inline void _Add_(const int col, const double a)
586  { SearchRow(col) += a; }
587  /// Set an entry in the "current row". See SetColPtr().
588  inline void _Set_(const int col, const double a)
589  { SearchRow(col) = a; }
590  /// Read the value of an entry in the "current row". See SetColPtr().
591  inline double _Get_(const int col) const;
592 
593  inline double &SearchRow(const int row, const int col);
594  inline void _Add_(const int row, const int col, const double a)
595  { SearchRow(row, col) += a; }
596  inline void _Set_(const int row, const int col, const double a)
597  { SearchRow(row, col) = a; }
598 
599  void Set(const int i, const int j, const double val);
600  void Add(const int i, const int j, const double val);
601 
602  void SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
603  const DenseMatrix &subm, int skip_zeros = 1);
604 
605  void SetSubMatrixTranspose(const Array<int> &rows, const Array<int> &cols,
606  const DenseMatrix &subm, int skip_zeros = 1);
607 
608  /** Insert the DenseMatrix into this SparseMatrix at the specified rows and
609  columns. If \c skip_zeros==0 , all entries from the DenseMatrix are
610  added including zeros. If \c skip_zeros==2 , no zeros are added to the
611  SparseMatrix regardless of their position in the matrix. Otherwise, the
612  default \c skip_zeros behavior is to omit the zero from the SparseMatrix
613  unless it would break the symmetric structure of the SparseMatrix. */
614  void AddSubMatrix(const Array<int> &rows, const Array<int> &cols,
615  const DenseMatrix &subm, int skip_zeros = 1);
616 
617  bool RowIsEmpty(const int row) const;
618 
619  /// Extract all column indices and values from a given row.
620  /** If the matrix is finalized (i.e. in CSR format), @a cols and @a srow will
621  simply be references to the specific portion of the #J and #A arrays.
622  As required by the AbstractSparseMatrix interface this method returns:
623  - 0, if @a cols and @a srow are copies of the values in the matrix, i.e.
624  when the matrix is open.
625  - 1, if @a cols and @a srow are views of the values in the matrix, i.e.
626  when the matrix is finalized.
627  @warning This method breaks the const-ness when the matrix is finalized
628  because it gives write access to the #J and #A arrays. */
629  virtual int GetRow(const int row, Array<int> &cols, Vector &srow) const;
630 
631  void SetRow(const int row, const Array<int> &cols, const Vector &srow);
632  void AddRow(const int row, const Array<int> &cols, const Vector &srow);
633 
634  void ScaleRow(const int row, const double scale);
635  /// this = diag(sl) * this;
636  void ScaleRows(const Vector & sl);
637  /// this = this * diag(sr);
638  void ScaleColumns(const Vector & sr);
639 
640  /** @brief Add the sparse matrix 'B' to '*this'. This operation will cause an
641  error if '*this' is finalized and 'B' has larger sparsity pattern. */
643 
644  /** @brief Add the sparse matrix 'B' scaled by the scalar 'a' into '*this'.
645  Only entries in the sparsity pattern of '*this' are added. */
646  void Add(const double a, const SparseMatrix &B);
647 
648  SparseMatrix &operator=(double a);
649 
650  SparseMatrix &operator*=(double a);
651 
652  /// Prints matrix to stream out.
653  void Print(std::ostream &out = mfem::out, int width_ = 4) const;
654 
655  /// Prints matrix in matlab format.
656  virtual void PrintMatlab(std::ostream &out = mfem::out) const;
657 
658  /// Prints matrix in Matrix Market sparse format.
659  void PrintMM(std::ostream &out = mfem::out) const;
660 
661  /// Prints matrix to stream out in hypre_CSRMatrix format.
662  void PrintCSR(std::ostream &out) const;
663 
664  /// Prints a sparse matrix to stream out in CSR format.
665  void PrintCSR2(std::ostream &out) const;
666 
667  /// Print various sparse matrix statistics.
668  void PrintInfo(std::ostream &out) const;
669 
670  /// Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix
671  double IsSymmetric() const;
672 
673  /// (*this) = 1/2 ((*this) + (*this)^t)
674  void Symmetrize();
675 
676  /// Returns the number of the nonzero elements in the matrix
677  virtual int NumNonZeroElems() const;
678 
679  double MaxNorm() const;
680 
681  /// Count the number of entries with |a_ij| <= tol.
682  int CountSmallElems(double tol) const;
683 
684  /// Count the number of entries that are NOT finite, i.e. Inf or Nan.
685  int CheckFinite() const;
686 
687  /// Set the graph ownership flag (I and J arrays).
688  void SetGraphOwner(bool ownij)
689  { I.SetHostPtrOwner(ownij); J.SetHostPtrOwner(ownij); }
690 
691  /// Set the data ownership flag (A array).
692  void SetDataOwner(bool owna) { A.SetHostPtrOwner(owna); }
693 
694  /// Get the graph ownership flag (I and J arrays).
695  bool OwnsGraph() const { return I.OwnsHostPtr() && J.OwnsHostPtr(); }
696 
697  /// Get the data ownership flag (A array).
698  bool OwnsData() const { return A.OwnsHostPtr(); }
699 
700  /// Lose the ownership of the graph (I, J) and data (A) arrays.
701  void LoseData() { SetGraphOwner(false); SetDataOwner(false); }
702 
703  void Swap(SparseMatrix &other);
704 
705  /// Destroys sparse matrix.
706  virtual ~SparseMatrix();
707 
708  Type GetType() const { return MFEM_SPARSEMAT; }
709 };
710 
711 inline std::ostream& operator<<(std::ostream& os, SparseMatrix const& mat)
712 {
713  mat.Print(os);
714  return os;
715 }
716 
717 /// Applies f() to each element of the matrix (after it is finalized).
718 void SparseMatrixFunction(SparseMatrix &S, double (*f)(double));
719 
720 
721 /// Transpose of a sparse matrix. A must be finalized.
722 SparseMatrix *Transpose(const SparseMatrix &A);
723 /// Transpose of a sparse matrix. A does not need to be a CSR matrix.
724 SparseMatrix *TransposeAbstractSparseMatrix (const AbstractSparseMatrix &A,
725  int useActualWidth);
726 
727 /// Matrix product A.B.
728 /** If @a OAB is not NULL, we assume it has the structure of A.B and store the
729  result in @a OAB. If @a OAB is NULL, we create a new SparseMatrix to store
730  the result and return a pointer to it.
731 
732  All matrices must be finalized. */
733 SparseMatrix *Mult(const SparseMatrix &A, const SparseMatrix &B,
734  SparseMatrix *OAB = NULL);
735 
736 /// C = A^T B
737 SparseMatrix *TransposeMult(const SparseMatrix &A, const SparseMatrix &B);
738 
739 /// Matrix product of sparse matrices. A and B do not need to be CSR matrices
740 SparseMatrix *MultAbstractSparseMatrix (const AbstractSparseMatrix &A,
741  const AbstractSparseMatrix &B);
742 
743 /// Matrix product A.B
744 DenseMatrix *Mult(const SparseMatrix &A, DenseMatrix &B);
745 
746 /// RAP matrix product (with R=P^T)
747 DenseMatrix *RAP(const SparseMatrix &A, DenseMatrix &P);
748 
749 /// RAP matrix product (with R=P^T)
750 DenseMatrix *RAP(DenseMatrix &A, const SparseMatrix &P);
751 
752 /** RAP matrix product (with P=R^T). ORAP is like OAB above.
753  All matrices must be finalized. */
754 SparseMatrix *RAP(const SparseMatrix &A, const SparseMatrix &R,
755  SparseMatrix *ORAP = NULL);
756 
757 /// General RAP with given R^T, A and P
758 SparseMatrix *RAP(const SparseMatrix &Rt, const SparseMatrix &A,
759  const SparseMatrix &P);
760 
761 /// Matrix multiplication A^t D A. All matrices must be finalized.
762 SparseMatrix *Mult_AtDA(const SparseMatrix &A, const Vector &D,
763  SparseMatrix *OAtDA = NULL);
764 
765 
766 /// Matrix addition result = A + B.
767 SparseMatrix * Add(const SparseMatrix & A, const SparseMatrix & B);
768 /// Matrix addition result = a*A + b*B
769 SparseMatrix * Add(double a, const SparseMatrix & A, double b,
770  const SparseMatrix & B);
771 /// Matrix addition result = sum_i A_i
772 SparseMatrix * Add(Array<SparseMatrix *> & Ai);
773 
774 /// B += alpha * A
775 void Add(const SparseMatrix &A, double alpha, DenseMatrix &B);
776 
777 /// Produces a block matrix with blocks A_{ij}*B
778 DenseMatrix *OuterProduct(const DenseMatrix &A, const DenseMatrix &B);
779 
780 /// Produces a block matrix with blocks A_{ij}*B
781 SparseMatrix *OuterProduct(const DenseMatrix &A, const SparseMatrix &B);
782 
783 /// Produces a block matrix with blocks A_{ij}*B
784 SparseMatrix *OuterProduct(const SparseMatrix &A, const DenseMatrix &B);
785 
786 /// Produces a block matrix with blocks A_{ij}*B
787 SparseMatrix *OuterProduct(const SparseMatrix &A, const SparseMatrix &B);
788 
789 
790 // Inline methods
791 
792 inline void SparseMatrix::SetColPtr(const int row) const
793 {
794  if (Rows)
795  {
796  if (ColPtrNode == NULL)
797  {
798  ColPtrNode = new RowNode *[width];
799  for (int i = 0; i < width; i++)
800  {
801  ColPtrNode[i] = NULL;
802  }
803  }
804  for (RowNode *node_p = Rows[row]; node_p != NULL; node_p = node_p->Prev)
805  {
806  ColPtrNode[node_p->Column] = node_p;
807  }
808  }
809  else
810  {
811  if (ColPtrJ == NULL)
812  {
813  ColPtrJ = new int[width];
814  for (int i = 0; i < width; i++)
815  {
816  ColPtrJ[i] = -1;
817  }
818  }
819  for (int j = I[row], end = I[row+1]; j < end; j++)
820  {
821  ColPtrJ[J[j]] = j;
822  }
823  }
824  current_row = row;
825 }
826 
827 inline void SparseMatrix::ClearColPtr() const
828 {
829  if (Rows)
830  {
831  for (RowNode *node_p = Rows[current_row]; node_p != NULL;
832  node_p = node_p->Prev)
833  {
834  ColPtrNode[node_p->Column] = NULL;
835  }
836  }
837  else
838  {
839  for (int j = I[current_row], end = I[current_row+1]; j < end; j++)
840  {
841  ColPtrJ[J[j]] = -1;
842  }
843  }
844 }
845 
846 inline double &SparseMatrix::SearchRow(const int col)
847 {
848  if (Rows)
849  {
850  RowNode *node_p = ColPtrNode[col];
851  if (node_p == NULL)
852  {
853 #ifdef MFEM_USE_MEMALLOC
854  node_p = NodesMem->Alloc();
855 #else
856  node_p = new RowNode;
857 #endif
858  node_p->Prev = Rows[current_row];
859  node_p->Column = col;
860  node_p->Value = 0.0;
861  Rows[current_row] = ColPtrNode[col] = node_p;
862  }
863  return node_p->Value;
864  }
865  else
866  {
867  const int j = ColPtrJ[col];
868  MFEM_VERIFY(j != -1, "Entry for column " << col << " is not allocated.");
869  return A[j];
870  }
871 }
872 
873 inline double SparseMatrix::_Get_(const int col) const
874 {
875  if (Rows)
876  {
877  RowNode *node_p = ColPtrNode[col];
878  return (node_p == NULL) ? 0.0 : node_p->Value;
879  }
880  else
881  {
882  const int j = ColPtrJ[col];
883  return (j == -1) ? 0.0 : A[j];
884  }
885 }
886 
887 inline double &SparseMatrix::SearchRow(const int row, const int col)
888 {
889  if (Rows)
890  {
891  RowNode *node_p;
892 
893  for (node_p = Rows[row]; 1; node_p = node_p->Prev)
894  {
895  if (node_p == NULL)
896  {
897 #ifdef MFEM_USE_MEMALLOC
898  node_p = NodesMem->Alloc();
899 #else
900  node_p = new RowNode;
901 #endif
902  node_p->Prev = Rows[row];
903  node_p->Column = col;
904  node_p->Value = 0.0;
905  Rows[row] = node_p;
906  break;
907  }
908  else if (node_p->Column == col)
909  {
910  break;
911  }
912  }
913  return node_p->Value;
914  }
915  else
916  {
917  int *Ip = I+row, *Jp = J;
918  for (int k = Ip[0], end = Ip[1]; k < end; k++)
919  {
920  if (Jp[k] == col)
921  {
922  return A[k];
923  }
924  }
925  MFEM_ABORT("Could not find entry for row = " << row << ", col = " << col);
926  }
927  return A[0];
928 }
929 
930 /// Specialization of the template function Swap<> for class SparseMatrix
931 template<> inline void Swap<SparseMatrix>(SparseMatrix &a, SparseMatrix &b)
932 {
933  a.Swap(b);
934 }
935 
936 } // namespace mfem
937 
938 #endif
RowNode ** ColPtrNode
Definition: sparsemat.hpp:77
Memory< int > I
Array with size (height+1) containing the row offsets.
Definition: sparsemat.hpp:62
Elem * Alloc()
Definition: mem_alloc.hpp:166
void _Add_(const int row, const int col, const double a)
Definition: sparsemat.hpp:594
cusparseStatus_t status
Definition: sparsemat.hpp:106
const double * HostReadData() const
Definition: sparsemat.hpp:277
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:733
SparseMatrix * TransposeAbstractSparseMatrix(const AbstractSparseMatrix &A, int useActualWidth)
Transpose of a sparse matrix. A does not need to be a CSR matrix.
Definition: sparsemat.cpp:3562
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:1588
void _Add_(const int col, const double a)
Add a value to an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:585
void Clear()
Clear the contents of the SparseMatrix.
Definition: sparsemat.hpp:209
void EliminateCol(int col, DiagonalPolicy dpolicy=DIAG_ZERO)
Eliminates the column col from the matrix.
Definition: sparsemat.cpp:1741
SparseMatrix * At
Transpose of A. Owned. Used to perform MultTranspose() on devices.
Definition: sparsemat.hpp:80
int Size() const
For backward compatibility, define Size() to be synonym of Height().
Definition: sparsemat.hpp:206
static MemoryClass GetHostMemoryClass()
Get the current Host MemoryClass. This is the MemoryClass used by most MFEM host Memory objects...
Definition: device.hpp:268
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:669
void PrintMM(std::ostream &out=mfem::out) const
Prints matrix in Matrix Market sparse format.
Definition: sparsemat.cpp:3314
void _Set_(const int row, const int col, const double a)
Definition: sparsemat.hpp:596
bool OwnsData() const
Get the data ownership flag (A array).
Definition: sparsemat.hpp:698
static size_t bufferSize
Definition: sparsemat.hpp:101
double GetRowNorml1(int irow) const
For i = irow compute .
Definition: sparsemat.cpp:1237
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:546
cusparseMatDescr_t descr
Definition: sparsemat.hpp:108
void MakeRef(const SparseMatrix &master)
Clear the contents of the SparseMatrix and make it a reference to master.
Definition: sparsemat.cpp:313
double & SearchRow(const int col)
Perform a fast search for an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:846
const Memory< int > & GetMemoryI() const
Definition: sparsemat.hpp:238
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:391
void PrintCSR(std::ostream &out) const
Prints matrix to stream out in hypre_CSRMatrix format.
Definition: sparsemat.cpp:3336
double GetJacobiScaling() const
Determine appropriate scaling for Jacobi iteration.
Definition: sparsemat.cpp:2555
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:344
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void SetColPtr(const int row) const
Initialize the SparseMatrix for fast access to the entries of the given row which becomes the "curren...
Definition: sparsemat.hpp:792
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:227
Abstract data type for sparse matrices.
Definition: matrix.hpp:73
double IsSymmetric() const
Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix.
Definition: sparsemat.cpp:1527
SparseMatrix & operator+=(const SparseMatrix &B)
Add the sparse matrix &#39;B&#39; to &#39;*this&#39;. This operation will cause an error if &#39;*this&#39; is finalized and ...
Definition: sparsemat.cpp:3141
T * Write(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for write access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true...
Definition: device.hpp:336
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
Abstract data type for matrix inverse.
Definition: matrix.hpp:62
static void * dBuffer
Definition: sparsemat.hpp:102
void EliminateRowColMultipleRHS(int rc, const Vector &sol, DenseMatrix &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Similar to EliminateRowCol(int, const double, Vector &, DiagonalPolicy), but multiple values for elim...
Definition: sparsemat.cpp:1953
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
y += A * x (default) or y += a * A * x
Definition: sparsemat.cpp:740
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Extract all column indices and values from a given row.
Definition: sparsemat.cpp:2935
void ClearGPUSparse()
Clear the cuSPARSE/hipSPARSE descriptors. This must be called after releasing the device memory of A...
Definition: sparsemat.cpp:81
double MaxNorm() const
Definition: sparsemat.cpp:1611
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: sparsemat.cpp:601
void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: sparsemat.cpp:3244
void LoseData()
Lose the ownership of the graph (I, J) and data (A) arrays.
Definition: sparsemat.hpp:701
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
Definition: sparsemat.cpp:2886
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:3022
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:405
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:457
Memory< double > & GetMemoryData()
Definition: sparsemat.hpp:269
Memory< int > & GetMemoryI()
Definition: sparsemat.hpp:237
cusparseDnVecDescr_t vecY_descr
Definition: sparsemat.hpp:113
double * HostReadWriteData()
Definition: sparsemat.hpp:281
void _Set_(const int col, const double a)
Set an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:588
SparseMatrix * Mult_AtDA(const SparseMatrix &A, const Vector &D, SparseMatrix *OAtDA)
Matrix multiplication A^t D A. All matrices must be finalized.
Definition: sparsemat.cpp:3930
const Memory< double > & GetMemoryData() const
Definition: sparsemat.hpp:270
void Gauss_Seidel_back(const Vector &x, Vector &y) const
Definition: sparsemat.cpp:2472
void ScaleRow(const int row, const double scale)
Definition: sparsemat.cpp:3054
RowNodeAlloc * NodesMem
Definition: sparsemat.hpp:84
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:232
virtual MatrixInverse * Inverse() const
This virtual method is not supported: it always returns NULL.
Definition: sparsemat.cpp:1688
Memory< int > & GetMemoryJ()
Definition: sparsemat.hpp:253
DenseMatrix * ToDenseMatrix() const
Produces a DenseMatrix from a SparseMatrix.
Definition: sparsemat.cpp:703
void UseGPUSparse(bool useGPUSparse_=true)
Runtime option to use cuSPARSE or hipSPARSE. Only valid when using a CUDA or HIP backend.
Definition: sparsemat.hpp:191
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
Definition: sparsemat.cpp:1569
void Jacobi3(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:2715
void MoveDiagonalFirst()
Move the diagonal entry to the first position in each row, preserving the order of the rest of the co...
Definition: sparsemat.cpp:577
void SetGraphOwner(bool ownij)
Set the graph ownership flag (I and J arrays).
Definition: sparsemat.hpp:688
ID for class SparseMatrix.
Definition: operator.hpp:286
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2321
void Set(const int i, const int j, const double val)
Definition: sparsemat.cpp:2768
hipsparseDnVecDescr_t vecY_descr
Definition: sparsemat.hpp:125
void GetBlocks(Array2D< SparseMatrix *> &blocks) const
Definition: sparsemat.cpp:1437
MFEM_DEPRECATED void ClearCuSparse()
Deprecated equivalent of ClearGPUSparse().
Definition: sparsemat.hpp:216
void ScaleColumns(const Vector &sr)
this = this * diag(sr);
Definition: sparsemat.cpp:3113
void GetRowSums(Vector &x) const
For all i compute .
Definition: sparsemat.cpp:1214
virtual ~SparseMatrix()
Destroys sparse matrix.
Definition: sparsemat.cpp:4202
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
bool Empty() const
Check if the SparseMatrix is empty.
Definition: sparsemat.hpp:219
bool isSorted
Are the columns sorted already.
Definition: sparsemat.hpp:88
void ClearColPtr() const
Reset the "current row" set by calling SetColPtr(). This method must be called between any two calls ...
Definition: sparsemat.hpp:827
Memory< double > A
Array with size I[height], containing the actual entries of the sparse matrix, as indexed by the I ar...
Definition: sparsemat.hpp:68
bool Empty() const
Return true if the Memory object is empty, see Reset().
SparseMatrix()
Create an empty SparseMatrix.
Definition: sparsemat.hpp:131
Data type sparse matrix.
Definition: sparsemat.hpp:50
static MemoryClass GetDeviceMemoryClass()
Get the current Device MemoryClass. This is the MemoryClass used by most MFEM device kernels to acces...
Definition: device.hpp:281
hipsparseStatus_t status
Definition: sparsemat.hpp:119
hipsparseSpMatDescr_t matA_descr
Definition: sparsemat.hpp:123
void SetWidth(int width_=-1)
Change the width of a SparseMatrix.
Definition: sparsemat.cpp:419
double b
Definition: lissajous.cpp:42
void EnsureMultTranspose() const
Ensures that the matrix is capable of performing MultTranspose(), AddMultTranspose(), and AbsMultTranspose().
Definition: sparsemat.cpp:974
const int * HostReadJ() const
Definition: sparsemat.hpp:261
virtual void PrintMatlab(std::ostream &out=mfem::out) const
Prints matrix in matlab format.
Definition: sparsemat.cpp:3292
int * ReadWriteI(bool on_dev=true)
Definition: sparsemat.hpp:243
int CountSmallElems(double tol) const
Count the number of entries with |a_ij| <= tol.
Definition: sparsemat.cpp:1636
static int SparseMatrixCount
Definition: sparsemat.hpp:100
const double * ReadData(bool on_dev=true) const
Definition: sparsemat.hpp:271
double InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Definition: sparsemat.cpp:1173
bool OwnsGraph() const
Get the graph ownership flag (I and J arrays).
Definition: sparsemat.hpp:695
void SetDiagIdentity()
If a row contains only one diag entry of zero, set it to 1.
Definition: sparsemat.cpp:2358
MFEM_DEPRECATED void UseCuSparse(bool useCuSparse_=true)
Deprecated equivalent of UseGPUSparse().
Definition: sparsemat.hpp:194
void ScaleRows(const Vector &sl)
this = diag(sl) * this;
Definition: sparsemat.cpp:3082
static cusparseHandle_t handle
Definition: sparsemat.hpp:107
void Swap< SparseMatrix >(SparseMatrix &a, SparseMatrix &b)
Specialization of the template function Swap<> for class SparseMatrix.
Definition: sparsemat.hpp:931
DenseMatrix * OuterProduct(const DenseMatrix &A, const DenseMatrix &B)
Produces a block matrix with blocks A_{ij}*B.
Definition: sparsemat.cpp:4074
void DiagScale(const Vector &b, Vector &x, double sc=1.0, bool use_abs_diag=false) const
x = sc b / A_ii. When use_abs_diag = true, |A_ii| is used.
Definition: sparsemat.cpp:2620
int * WriteI(bool on_dev=true)
Definition: sparsemat.hpp:241
int Capacity() const
Return the size of the allocated memory.
Set the diagonal value to one.
Definition: operator.hpp:50
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3252
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true...
Definition: device.hpp:319
bool RowIsEmpty(const int row) const
Definition: sparsemat.cpp:2914
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:552
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:2708
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2722
void Add(const int i, const int j, const double val)
Definition: sparsemat.cpp:2787
Dynamic 2D array using row-major layout.
Definition: array.hpp:354
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2806
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:908
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:283
void SparseMatrixFunction(SparseMatrix &S, double(*f)(double))
Applies f() to each element of the matrix (after it is finalized).
Definition: sparsemat.cpp:3495
const int * ReadI(bool on_dev=true) const
Definition: sparsemat.hpp:239
double _Get_(const int col) const
Read the value of an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:873
void EliminateRowColDiag(int rc, double value)
Perform elimination and set the diagonal entry to the given value.
Definition: sparsemat.cpp:2158
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
void SetSubMatrixTranspose(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2845
void Gauss_Seidel_forw(const Vector &x, Vector &y) const
Gauss-Seidel forward and backward iterations over a vector x.
Definition: sparsemat.cpp:2388
void SetHostPtrOwner(bool own) const
Set/clear the ownership flag for the host pointer. Ownership indicates whether the pointer will be de...
int ActualWidth() const
Returns the actual Width of the matrix.
Definition: sparsemat.cpp:3469
void EliminateRowCol(int rc, const double sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate row rc and column rc and modify the rhs using sol.
Definition: sparsemat.cpp:1854
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
void Threshold(double tol, bool fix_empty_rows=false)
Remove entries smaller in absolute value than a given tolerance tol. If fix_empty_rows is true...
Definition: sparsemat.cpp:1261
const int * GetI() const
Return the array I, const version.
Definition: sparsemat.hpp:224
double * HostWriteData()
Definition: sparsemat.hpp:279
const double * GetData() const
Return the element data, i.e. the array A, const version.
Definition: sparsemat.hpp:234
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int * HostReadWriteJ()
Definition: sparsemat.hpp:265
RowNode ** Rows
Array of linked lists, one for every row. This array represents the linked list (LIL) storage format...
Definition: sparsemat.hpp:73
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc, bool use_abs_diag=false) const
Definition: sparsemat.cpp:2588
const int * HostReadI() const
Definition: sparsemat.hpp:245
hipsparseDnVecDescr_t vecX_descr
Definition: sparsemat.hpp:124
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:1026
int * WriteJ(bool on_dev=true)
Definition: sparsemat.hpp:257
void EliminateRow(int row, const double sol, Vector &rhs)
Eliminates a column from the transpose matrix.
Definition: sparsemat.cpp:1693
double a
Definition: lissajous.cpp:41
void BuildTranspose() const
Build and store internally the transpose of this matrix which will be used in the methods AddMultTran...
Definition: sparsemat.cpp:960
MemAlloc< RowNode, 1024 > RowNodeAlloc
Definition: sparsemat.hpp:83
T * ReadWrite(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read+write access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true, or the mfem::Device&#39;s HostMemoryClass, otherwise.
Definition: device.hpp:353
SparseMatrix * TransposeMult(const SparseMatrix &A, const SparseMatrix &B)
C = A^T B.
Definition: sparsemat.cpp:3767
virtual void EliminateZeroRows(const double threshold=1e-12)
If a row contains only zeros, set its diagonal to 1.
Definition: sparsemat.cpp:2369
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: sparsemat.hpp:342
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:47
double & operator()(int i, int j)
Returns reference to A[i][j].
Definition: sparsemat.cpp:611
void AbsMult(const Vector &x, Vector &y) const
y = |A| * x, using entry-wise absolute values of matrix A
Definition: sparsemat.cpp:1079
void PrintInfo(std::ostream &out) const
Print various sparse matrix statistics.
Definition: sparsemat.cpp:3385
void EliminateCols(const Array< int > &cols, const Vector *x=NULL, Vector *b=NULL)
Eliminate all columns i for which cols[i] != 0.
Definition: sparsemat.cpp:1781
void OverrideSize(int height_, int width_)
Sets the height and width of the matrix.
Definition: sparsemat.cpp:297
const int * ReadJ(bool on_dev=true) const
Definition: sparsemat.hpp:255
bool ColumnsAreSorted() const
Returns whether or not the columns are sorted.
Definition: sparsemat.hpp:554
SparseMatrix & operator=(const SparseMatrix &rhs)
Assignment operator: deep copy.
Definition: sparsemat.cpp:303
void AbsMultTranspose(const Vector &x, Vector &y) const
y = |At| * x, using entry-wise absolute values of the transpose of matrix A
Definition: sparsemat.cpp:1128
SparseMatrix & operator*=(double a)
Definition: sparsemat.cpp:3220
void ResetTranspose() const
Definition: sparsemat.cpp:968
void SetDataOwner(bool owna)
Set the data ownership flag (A array).
Definition: sparsemat.hpp:692
void PrintCSR2(std::ostream &out) const
Prints a sparse matrix to stream out in CSR format.
Definition: sparsemat.cpp:3360
int MaxRowSize() const
Returns the maximum number of elements among all rows.
Definition: sparsemat.cpp:367
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
Definition: sparsemat.hpp:711
bool OwnsHostPtr() const
Return true if the host pointer is owned. Ownership indicates whether the pointer will be deleted by ...
const int * GetJ() const
Return the array J, const version.
Definition: sparsemat.hpp:229
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2975
cusparseSpMatDescr_t matA_descr
Definition: sparsemat.hpp:111
const double alpha
Definition: ex15.cpp:369
double * ReadWriteData(bool on_dev=true)
Definition: sparsemat.hpp:275
int * ReadWriteJ(bool on_dev=true)
Definition: sparsemat.hpp:259
static hipsparseHandle_t handle
Definition: sparsemat.hpp:120
Vector data type.
Definition: vector.hpp:60
void PartAddMult(const Array< int > &rows, const Vector &x, Vector &y, const double a=1.0) const
Definition: sparsemat.cpp:1008
cusparseDnVecDescr_t vecX_descr
Definition: sparsemat.hpp:112
int CheckFinite() const
Count the number of entries that are NOT finite, i.e. Inf or Nan.
Definition: sparsemat.cpp:1664
const Memory< int > & GetMemoryJ() const
Definition: sparsemat.hpp:254
int * HostReadWriteI()
Definition: sparsemat.hpp:249
cusparseMatDescr_t matA_descr
Definition: sparsemat.hpp:115
SparseMatrix * MultAbstractSparseMatrix(const AbstractSparseMatrix &A, const AbstractSparseMatrix &B)
Matrix product of sparse matrices. A and B do not need to be CSR matrices.
Definition: sparsemat.cpp:3775
Memory< int > J
Array with size I[height], containing the column indices for all matrix entries, as indexed by the I ...
Definition: sparsemat.hpp:65
void Swap(SparseMatrix &other)
Definition: sparsemat.cpp:4182
MemoryClass
Memory classes identify sets of memory types.
Definition: mem_manager.hpp:73
Type GetType() const
Definition: sparsemat.hpp:708
Set the diagonal value to zero.
Definition: operator.hpp:49
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void EliminateBC(const Array< int > &ess_dofs, DiagonalPolicy diag_policy)
Eliminate essential (Dirichlet) boundary conditions.
Definition: sparsemat.cpp:2315
virtual void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
y += At * x (default) or y += a * At * x
Definition: sparsemat.cpp:915
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Definition: sparsemat.cpp:982
double * WriteData(bool on_dev=true)
Definition: sparsemat.hpp:273
void BooleanMultTranspose(const Array< int > &x, Array< int > &y) const
y = At * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:1056