1 // Copyright (c) 2010-2021, 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 namespace mfem
26 {
27
28 class
29 #if defined(__alignas_is_defined)
30 alignas(double)
31 #endif
32 RowNode
33 {
34 public:
35 double Value;
36 RowNode *Prev;
37 int Column;
38 };
39
40 /// Data type sparse matrix
41 class SparseMatrix : public AbstractSparseMatrix
42 {
43 protected:
44 /// @name Arrays used by the CSR storage format.
45 /** */
46 ///@{
47 /// @brief %Array with size (#height+1) containing the row offsets.
48 /** The data for row r, 0 <= r < height, is at offsets j, I[r] <= j < I[r+1].
49 The offsets, j, are indices in the #J and #A arrays. The first entry in
50 this array is always zero, I[0] = 0, and the last entry, I[height], gives
51 the total number of entries stored (at a minimum, all nonzeros must be
52 represented) in the sparse matrix. */
53 Memory<int> I;
54 /** @brief %Array with size #I[#height], containing the column indices for
55 all matrix entries, as indexed by the #I array. */
56 Memory<int> J;
57 /** @brief %Array with size #I[#height], containing the actual entries of the
58 sparse matrix, as indexed by the #I array. */
59 Memory<double> A;
60 ///@}
61
62 /** @brief %Array of linked lists, one for every row. This array represents
63 the linked list (LIL) storage format. */
64 RowNode **Rows;
65
66 mutable int current_row;
67 mutable int* ColPtrJ;
68 mutable RowNode ** ColPtrNode;
69
70 /// Transpose of A. Owned. Used to perform MultTranspose() on devices.
71 mutable SparseMatrix *At;
72
73 #ifdef MFEM_USE_MEMALLOC
74 typedef MemAlloc <RowNode, 1024> RowNodeAlloc;
75 RowNodeAlloc * NodesMem;
76 #endif
77
78 /// Are the columns sorted already.
79 bool isSorted;
80
81 void Destroy(); // Delete all owned data
82 void SetEmpty(); // Init all entries with empty values
83
84 bool useCuSparse{true}; // Use cuSPARSE if available
85
86 // Initialize cuSPARSE
87 void InitCuSparse();
88
89 #ifdef MFEM_USE_CUDA
90 cusparseStatus_t status;
91 static cusparseHandle_t handle;
92 cusparseMatDescr_t descr=0;
93 static int SparseMatrixCount;
94 static size_t bufferSize;
95 static void *dBuffer;
96 mutable bool initBuffers{false};
97 #if CUDA_VERSION >= 10010
98 mutable cusparseSpMatDescr_t matA_descr;
99 mutable cusparseDnVecDescr_t vecX_descr;
100 mutable cusparseDnVecDescr_t vecY_descr;
101 #else
102 mutable cusparseMatDescr_t matA_descr;
103 #endif
104 #endif
105
106 public:
107 /// Create an empty SparseMatrix.
SparseMatrix()108 SparseMatrix()
109 {
110 SetEmpty();
111
112 InitCuSparse();
113 }
114
115 /** @brief Create a sparse matrix with flexible sparsity structure using a
116 row-wise linked list (LIL) format. */
117 /** New entries are added as needed by methods like AddSubMatrix(),
118 SetSubMatrix(), etc. Calling Finalize() will convert the SparseMatrix to
119 the more compact compressed sparse row (CSR) format. */
120 explicit SparseMatrix(int nrows, int ncols = -1);
121
122 /** @brief Create a sparse matrix in CSR format. Ownership of @a i, @a j, and
123 @a data is transferred to the SparseMatrix. */
124 SparseMatrix(int *i, int *j, double *data, int m, int n);
125
126 /** @brief Create a sparse matrix in CSR format. Ownership of @a i, @a j, and
127 @a data is optionally transferred to the SparseMatrix. */
128 /** If the parameter @a data is NULL, then the internal #A array is allocated
129 by this constructor (initializing it with zeros and taking ownership,
130 regardless of the parameter @a owna). */
131 SparseMatrix(int *i, int *j, double *data, int m, int n, bool ownij,
132 bool owna, bool issorted);
133
134 /** @brief Create a sparse matrix in CSR format where each row has space
135 allocated for exactly @a rowsize entries. */
136 /** SetRow() can then be called or the #I, #J, #A arrays can be used
137 directly. */
138 SparseMatrix(int nrows, int ncols, int rowsize);
139
140 /// Copy constructor (deep copy).
141 /** If @a mat is finalized and @a copy_graph is false, the #I and #J arrays
142 will use a shallow copy (copy the pointers only) without transferring
143 ownership.
144 If @a mt is MemoryType::PRESERVE the memory type of the resulting
145 SparseMatrix's #I, #J, and #A arrays will be the same as @a mat,
146 otherwise the type will be @a mt for those arrays that are deep
147 copied. */
148 SparseMatrix(const SparseMatrix &mat, bool copy_graph = true,
149 MemoryType mt = MemoryType::PRESERVE);
150
151 /// Create a SparseMatrix with diagonal @a v, i.e. A = Diag(v)
152 SparseMatrix(const Vector & v);
153
154 // Runtime option to use cuSPARSE. Only valid when using a CUDA backend.
UseCuSparse(bool useCuSparse_=true)155 void UseCuSparse(bool useCuSparse_ = true) { useCuSparse = useCuSparse_;}
156
157 /// Assignment operator: deep copy
158 SparseMatrix& operator=(const SparseMatrix &rhs);
159
160 /** @brief Clear the contents of the SparseMatrix and make it a reference to
161 @a master */
162 /** After this call, the matrix will point to the same data as @a master but
163 it will not own its data. The @a master must be finalized. */
164 void MakeRef(const SparseMatrix &master);
165
166 /// For backward compatibility, define Size() to be synonym of Height().
Size() const167 int Size() const { return Height(); }
168
169 /// Clear the contents of the SparseMatrix.
Clear()170 void Clear() { Destroy(); SetEmpty(); }
171
172 /** @brief Clear the CuSparse descriptors.
173 This must be called after releasing the device memory of A. */
174 void ClearCuSparse();
175
176 /// Check if the SparseMatrix is empty.
Empty() const177 bool Empty() const { return (A == NULL) && (Rows == NULL); }
178
179 /// Return the array #I.
GetI()180 inline int *GetI() { return I; }
181 /// Return the array #I, const version.
GetI() const182 inline const int *GetI() const { return I; }
183
184 /// Return the array #J.
GetJ()185 inline int *GetJ() { return J; }
186 /// Return the array #J, const version.
GetJ() const187 inline const int *GetJ() const { return J; }
188
189 /// Return the element data, i.e. the array #A.
GetData()190 inline double *GetData() { return A; }
191 /// Return the element data, i.e. the array #A, const version.
GetData() const192 inline const double *GetData() const { return A; }
193
194 // Memory access methods for the #I array.
GetMemoryI()195 Memory<int> &GetMemoryI() { return I; }
GetMemoryI() const196 const Memory<int> &GetMemoryI() const { return I; }
ReadI(bool on_dev=true) const197 const int *ReadI(bool on_dev = true) const
198 { return mfem::Read(I, Height()+1, on_dev); }
WriteI(bool on_dev=true)199 int *WriteI(bool on_dev = true)
200 { return mfem::Write(I, Height()+1, on_dev); }
ReadWriteI(bool on_dev=true)201 int *ReadWriteI(bool on_dev = true)
202 { return mfem::ReadWrite(I, Height()+1, on_dev); }
HostReadI() const203 const int *HostReadI() const
204 { return mfem::Read(I, Height()+1, false); }
HostWriteI()205 int *HostWriteI()
206 { return mfem::Write(I, Height()+1, false); }
HostReadWriteI()207 int *HostReadWriteI()
208 { return mfem::ReadWrite(I, Height()+1, false); }
209
210 // Memory access methods for the #J array.
GetMemoryJ()211 Memory<int> &GetMemoryJ() { return J; }
GetMemoryJ() const212 const Memory<int> &GetMemoryJ() const { return J; }
ReadJ(bool on_dev=true) const213 const int *ReadJ(bool on_dev = true) const
214 { return mfem::Read(J, J.Capacity(), on_dev); }
WriteJ(bool on_dev=true)215 int *WriteJ(bool on_dev = true)
216 { return mfem::Write(J, J.Capacity(), on_dev); }
ReadWriteJ(bool on_dev=true)217 int *ReadWriteJ(bool on_dev = true)
218 { return mfem::ReadWrite(J, J.Capacity(), on_dev); }
HostReadJ() const219 const int *HostReadJ() const
220 { return mfem::Read(J, J.Capacity(), false); }
HostWriteJ()221 int *HostWriteJ()
222 { return mfem::Write(J, J.Capacity(), false); }
HostReadWriteJ()223 int *HostReadWriteJ()
224 { return mfem::ReadWrite(J, J.Capacity(), false); }
225
226 // Memory access methods for the #A array.
GetMemoryData()227 Memory<double> &GetMemoryData() { return A; }
GetMemoryData() const228 const Memory<double> &GetMemoryData() const { return A; }
ReadData(bool on_dev=true) const229 const double *ReadData(bool on_dev = true) const
230 { return mfem::Read(A, A.Capacity(), on_dev); }
WriteData(bool on_dev=true)231 double *WriteData(bool on_dev = true)
232 { return mfem::Write(A, A.Capacity(), on_dev); }
ReadWriteData(bool on_dev=true)233 double *ReadWriteData(bool on_dev = true)
234 { return mfem::ReadWrite(A, A.Capacity(), on_dev); }
HostReadData() const235 const double *HostReadData() const
236 { return mfem::Read(A, A.Capacity(), false); }
HostWriteData()237 double *HostWriteData()
238 { return mfem::Write(A, A.Capacity(), false); }
HostReadWriteData()239 double *HostReadWriteData()
240 { return mfem::ReadWrite(A, A.Capacity(), false); }
241
242 /// Returns the number of elements in row @a i.
243 int RowSize(const int i) const;
244
245 /// Returns the maximum number of elements among all rows.
246 int MaxRowSize() const;
247
248 /// Return a pointer to the column indices in a row.
249 int *GetRowColumns(const int row);
250 /// Return a pointer to the column indices in a row, const version.
251 const int *GetRowColumns(const int row) const;
252
253 /// Return a pointer to the entries in a row.
254 double *GetRowEntries(const int row);
255 /// Return a pointer to the entries in a row, const version.
256 const double *GetRowEntries(const int row) const;
257
258 /// Change the width of a SparseMatrix.
259 /*!
260 * If width_ = -1 (DEFAULT), this routine will set the new width
261 * to the actual Width of the matrix awidth = max(J) + 1.
262 * Values 0 <= width_ < awidth are not allowed (error check in Debug Mode only)
263 *
264 * This method can be called for matrices finalized or not.
265 */
266 void SetWidth(int width_ = -1);
267
268 /// Returns the actual Width of the matrix.
269 /*! This method can be called for matrices finalized or not. */
270 int ActualWidth() const;
271
272 /// Sort the column indices corresponding to each row.
273 void SortColumnIndices();
274
275 /** @brief Move the diagonal entry to the first position in each row,
276 preserving the order of the rest of the columns. */
277 void MoveDiagonalFirst();
278
279 /// Returns reference to a_{ij}.
280 virtual double &Elem(int i, int j);
281
282 /// Returns constant reference to a_{ij}.
283 virtual const double &Elem(int i, int j) const;
284
285 /// Returns reference to A[i][j].
286 double &operator()(int i, int j);
287
288 /// Returns reference to A[i][j].
289 const double &operator()(int i, int j) const;
290
291 /// Returns the Diagonal of A
292 void GetDiag(Vector & d) const;
293
294 /// Produces a DenseMatrix from a SparseMatrix
295 DenseMatrix *ToDenseMatrix() const;
296
297 /// Produces a DenseMatrix from a SparseMatrix
298 void ToDenseMatrix(DenseMatrix & B) const;
299
GetMemoryClass() const300 virtual MemoryClass GetMemoryClass() const
301 {
302 return Finalized() ?
303 Device::GetDeviceMemoryClass() : Device::GetHostMemoryClass();
304 }
305
306 /// Matrix vector multiplication.
307 virtual void Mult(const Vector &x, Vector &y) const;
308
309 /// y += A * x (default) or y += a * A * x
310 void AddMult(const Vector &x, Vector &y, const double a = 1.0) const;
311
312 /// Multiply a vector with the transposed matrix. y = At * x
313 void MultTranspose(const Vector &x, Vector &y) const;
314
315 /// y += At * x (default) or y += a * At * x
316 void AddMultTranspose(const Vector &x, Vector &y,
317 const double a = 1.0) const;
318
319 /** @brief Build and store internally the transpose of this matrix which will
320 be used in the methods AddMultTranspose() and MultTranspose(). */
321 /** If this method has been called, the internal transpose matrix will be
322 used to perform the action of the transpose matrix in AddMultTranspose(),
323 and MultTranspose().
324
325 Warning: any changes in this matrix will invalidate the internal
326 transpose. To rebuild the transpose, call ResetTranspose() followed by a
327 call to this method. If the internal transpose is already built, this
328 method has no effect.
329
330 When any non-default backend is enabled, i.e. Device::IsEnabled() is
331 true, the methods AddMultTranspose(), and MultTranspose(), require the
332 internal transpose to be built. If that is not the case (i.e. the
333 internal transpose is not built), these methods will raise an error with
334 an appropriate message pointing to this method. When using the default
335 backend, calling this method is optional.
336
337 This method can only be used when the sparse matrix is finalized. */
338 void BuildTranspose() const;
339
340 /** Reset (destroy) the internal transpose matrix. See BuildTranspose() for
341 more details. */
342 void ResetTranspose() const;
343
344 void PartMult(const Array<int> &rows, const Vector &x, Vector &y) const;
345 void PartAddMult(const Array<int> &rows, const Vector &x, Vector &y,
346 const double a=1.0) const;
347
348 /// y = A * x, treating all entries as booleans (zero=false, nonzero=true).
349 /** The actual values stored in the data array, #A, are not used - this means
350 that all entries in the sparsity pattern are considered to be true by
351 this method. */
352 void BooleanMult(const Array<int> &x, Array<int> &y) const;
353
354 /// y = At * x, treating all entries as booleans (zero=false, nonzero=true).
355 /** The actual values stored in the data array, #A, are not used - this means
356 that all entries in the sparsity pattern are considered to be true by
357 this method. */
358 void BooleanMultTranspose(const Array<int> &x, Array<int> &y) const;
359
360 /// y = |A| * x, using entry-wise absolute values of matrix A
361 void AbsMult(const Vector &x, Vector &y) const;
362
363 /// y = |At| * x, using entry-wise absolute values of the transpose of matrix A
364 void AbsMultTranspose(const Vector &x, Vector &y) const;
365
366 /// Compute y^t A x
367 double InnerProduct(const Vector &x, const Vector &y) const;
368
369 /// For all i compute \f$ x_i = \sum_j A_{ij} \f$
370 void GetRowSums(Vector &x) const;
371 /// For i = irow compute \f$ x_i = \sum_j | A_{i, j} | \f$
372 double GetRowNorml1(int irow) const;
373
374 /// This virtual method is not supported: it always returns NULL.
375 virtual MatrixInverse *Inverse() const;
376
377 /// Eliminates a column from the transpose matrix.
378 void EliminateRow(int row, const double sol, Vector &rhs);
379
380 /// Eliminates a row from the matrix.
381 /*!
382 * - If @a dpolicy = #DIAG_ZERO, all the entries in the row will be set to 0.
383 * - If @a dpolicy = #DIAG_ONE (matrix must be square), the diagonal entry
384 * will be set equal to 1 and all other entries in the row to 0.
385 * - The policy #DIAG_KEEP is not supported.
386 */
387 void EliminateRow(int row, DiagonalPolicy dpolicy = DIAG_ZERO);
388
389 /// Eliminates the column @a col from the matrix.
390 /** - If @a dpolicy = #DIAG_ZERO, all entries in the column will be set to 0.
391 - If @a dpolicy = #DIAG_ONE (matrix must be square), the diagonal entry
392 will be set equal to 1 and all other entries in the column to 0.
393 - The policy #DIAG_KEEP is not supported. */
394 void EliminateCol(int col, DiagonalPolicy dpolicy = DIAG_ZERO);
395
396 /// Eliminate all columns i for which @a cols[i] != 0.
397 /** Elimination of a column means that all entries in the column are set to
398 zero. In addition, if the pointers @a x and @a b are not NULL, the
399 eliminated matrix entries are multiplied by the corresponding solution
400 value in @a *x and subtracted from the r.h.s. vector, @a *b. */
401 void EliminateCols(const Array<int> &cols, const Vector *x = NULL,
402 Vector *b = NULL);
403
404 /** @brief Similar to EliminateCols + save the eliminated entries into
405 @a Ae so that (*this) + Ae is equal to the original matrix. */
406 void EliminateCols(const Array<int> &col_marker, SparseMatrix &Ae);
407
408 /// Eliminate row @a rc and column @a rc and modify the @a rhs using @a sol.
409 /** Eliminates the column @a rc to the @a rhs, deletes the row @a rc and
410 replaces the element (rc,rc) with 1.0; assumes that element (i,rc)
411 is assembled if and only if the element (rc,i) is assembled.
412 By default, elements (rc,rc) are set to 1.0, although this behavior
413 can be adjusted by changing the @a dpolicy parameter. */
414 void EliminateRowCol(int rc, const double sol, Vector &rhs,
415 DiagonalPolicy dpolicy = DIAG_ONE);
416
417 /** @brief Similar to
418 EliminateRowCol(int, const double, Vector &, DiagonalPolicy), but
419 multiple values for eliminated unknowns are accepted, and accordingly
420 multiple right-hand-sides are used. */
421 void EliminateRowColMultipleRHS(int rc, const Vector &sol,
422 DenseMatrix &rhs,
423 DiagonalPolicy dpolicy = DIAG_ONE);
424
425 /// Perform elimination and set the diagonal entry to the given value
426 void EliminateRowColDiag(int rc, double value);
427
428 /// Eliminate row @a rc and column @a rc.
429 void EliminateRowCol(int rc, DiagonalPolicy dpolicy = DIAG_ONE);
430
431 /** @brief Similar to EliminateRowCol(int, DiagonalPolicy) + save the
432 eliminated entries into @a Ae so that (*this) + Ae is equal to the
433 original matrix */
434 void EliminateRowCol(int rc, SparseMatrix &Ae,
435 DiagonalPolicy dpolicy = DIAG_ONE);
436
437 /// If a row contains only one diag entry of zero, set it to 1.
438 void SetDiagIdentity();
439 /// If a row contains only zeros, set its diagonal to 1.
440 virtual void EliminateZeroRows(const double threshold = 1e-12);
441
442 /// Gauss-Seidel forward and backward iterations over a vector x.
443 void Gauss_Seidel_forw(const Vector &x, Vector &y) const;
444 void Gauss_Seidel_back(const Vector &x, Vector &y) const;
445
446 /// Determine appropriate scaling for Jacobi iteration
447 double GetJacobiScaling() const;
448 /** One scaled Jacobi iteration for the system A x = b.
449 x1 = x0 + sc D^{-1} (b - A x0) where D is the diag of A. */
450 void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const;
451
452 void DiagScale(const Vector &b, Vector &x, double sc = 1.0) const;
453
454 /** x1 = x0 + sc D^{-1} (b - A x0) where \f$ D_{ii} = \sum_j |A_{ij}| \f$. */
455 void Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
456 double sc = 1.0) const;
457
458 /** x1 = x0 + sc D^{-1} (b - A x0) where \f$ D_{ii} = \sum_j A_{ij} \f$. */
459 void Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
460 double sc = 1.0) const;
461
462 /** @brief Finalize the matrix initialization, switching the storage format
463 from LIL to CSR. */
464 /** This method should be called once, after the matrix has been initialized.
465 Internally, this method converts the matrix from row-wise linked list
466 (LIL) format into CSR (compressed sparse row) format. */
Finalize(int skip_zeros=1)467 virtual void Finalize(int skip_zeros = 1) { Finalize(skip_zeros, false); }
468
469 /// A slightly more general version of the Finalize(int) method.
470 void Finalize(int skip_zeros, bool fix_empty_rows);
471
472 /// Returns whether or not CSR format has been finalized.
Finalized() const473 bool Finalized() const { return !A.Empty(); }
474 /// Returns whether or not the columns are sorted.
ColumnsAreSorted() const475 bool ColumnsAreSorted() const { return isSorted; }
476
477 /** @brief Remove entries smaller in absolute value than a given tolerance
478 @a tol. If @a fix_empty_rows is true, a zero value is inserted in the
479 diagonal entry (for square matrices only) */
480 void Threshold(double tol, bool fix_empty_rows = false);
481
482 /** Split the matrix into M x N blocks of sparse matrices in CSR format.
483 The 'blocks' array is M x N (i.e. M and N are determined by its
484 dimensions) and its entries are overwritten by the new blocks. */
485 void GetBlocks(Array2D<SparseMatrix *> &blocks) const;
486
487 void GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
488 DenseMatrix &subm) const;
489
490 /** @brief Initialize the SparseMatrix for fast access to the entries of the
491 given @a row which becomes the "current row". */
492 /** Fast access to the entries of the "current row" can be performed using
493 the methods: SearchRow(const int), _Add_(const int, const double),
494 _Set_(const int, const double), and _Get_(const int). */
495 inline void SetColPtr(const int row) const;
496 /** @brief Reset the "current row" set by calling SetColPtr(). This method
497 must be called between any two calls to SetColPtr(). */
498 inline void ClearColPtr() const;
499 /// Perform a fast search for an entry in the "current row". See SetColPtr().
500 /** If the matrix is not finalized and the entry is not found in the
501 SparseMatrix, it will be added to the sparsity pattern initialized with
502 zero. If the matrix is finalized and the entry is not found, an error
503 will be generated. */
504 inline double &SearchRow(const int col);
505 /// Add a value to an entry in the "current row". See SetColPtr().
_Add_(const int col,const double a)506 inline void _Add_(const int col, const double a)
507 { SearchRow(col) += a; }
508 /// Set an entry in the "current row". See SetColPtr().
_Set_(const int col,const double a)509 inline void _Set_(const int col, const double a)
510 { SearchRow(col) = a; }
511 /// Read the value of an entry in the "current row". See SetColPtr().
512 inline double _Get_(const int col) const;
513
514 inline double &SearchRow(const int row, const int col);
_Add_(const int row,const int col,const double a)515 inline void _Add_(const int row, const int col, const double a)
516 { SearchRow(row, col) += a; }
_Set_(const int row,const int col,const double a)517 inline void _Set_(const int row, const int col, const double a)
518 { SearchRow(row, col) = a; }
519
520 void Set(const int i, const int j, const double a);
521 void Add(const int i, const int j, const double a);
522
523 void SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
524 const DenseMatrix &subm, int skip_zeros = 1);
525
526 void SetSubMatrixTranspose(const Array<int> &rows, const Array<int> &cols,
527 const DenseMatrix &subm, int skip_zeros = 1);
528
529 /** Insert the DenseMatrix into this SparseMatrix at the specified rows and
530 columns. If \c skip_zeros==0 , all entries from the DenseMatrix are
531 added including zeros. If \c skip_zeros==2 , no zeros are added to the
532 SparseMatrix regardless of their position in the matrix. Otherwise, the
533 default \c skip_zeros behavior is to omit the zero from the SparseMatrix
534 unless it would break the symmetric structure of the SparseMatrix. */
535 void AddSubMatrix(const Array<int> &rows, const Array<int> &cols,
536 const DenseMatrix &subm, int skip_zeros = 1);
537
538 bool RowIsEmpty(const int row) const;
539
540 /// Extract all column indices and values from a given row.
541 /** If the matrix is finalized (i.e. in CSR format), @a cols and @a srow will
542 simply be references to the specific portion of the #J and #A arrays.
543 As required by the AbstractSparseMatrix interface this method returns:
544 - 0, if @a cols and @a srow are copies of the values in the matrix, i.e.
545 when the matrix is open.
546 - 1, if @a cols and @a srow are views of the values in the matrix, i.e.
547 when the matrix is finalized.
548 @warning This method breaks the const-ness when the matrix is finalized
549 because it gives write access to the #J and #A arrays. */
550 virtual int GetRow(const int row, Array<int> &cols, Vector &srow) const;
551
552 void SetRow(const int row, const Array<int> &cols, const Vector &srow);
553 void AddRow(const int row, const Array<int> &cols, const Vector &srow);
554
555 void ScaleRow(const int row, const double scale);
556 /// this = diag(sl) * this;
557 void ScaleRows(const Vector & sl);
558 /// this = this * diag(sr);
559 void ScaleColumns(const Vector & sr);
560
561 /** @brief Add the sparse matrix 'B' to '*this'. This operation will cause an
562 error if '*this' is finalized and 'B' has larger sparsity pattern. */
563 SparseMatrix &operator+=(const SparseMatrix &B);
564
565 /** @brief Add the sparse matrix 'B' scaled by the scalar 'a' into '*this'.
566 Only entries in the sparsity pattern of '*this' are added. */
567 void Add(const double a, const SparseMatrix &B);
568
569 SparseMatrix &operator=(double a);
570
571 SparseMatrix &operator*=(double a);
572
573 /// Prints matrix to stream out.
574 void Print(std::ostream &out = mfem::out, int width_ = 4) const;
575
576 /// Prints matrix in matlab format.
577 void PrintMatlab(std::ostream &out = mfem::out) const;
578
579 /// Prints matrix in Matrix Market sparse format.
580 void PrintMM(std::ostream &out = mfem::out) const;
581
582 /// Prints matrix to stream out in hypre_CSRMatrix format.
583 void PrintCSR(std::ostream &out) const;
584
585 /// Prints a sparse matrix to stream out in CSR format.
586 void PrintCSR2(std::ostream &out) const;
587
588 /// Print various sparse matrix statistics.
589 void PrintInfo(std::ostream &out) const;
590
591 /// Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix
592 double IsSymmetric() const;
593
594 /// (*this) = 1/2 ((*this) + (*this)^t)
595 void Symmetrize();
596
597 /// Returns the number of the nonzero elements in the matrix
598 virtual int NumNonZeroElems() const;
599
600 double MaxNorm() const;
601
602 /// Count the number of entries with |a_ij| <= tol.
603 int CountSmallElems(double tol) const;
604
605 /// Count the number of entries that are NOT finite, i.e. Inf or Nan.
606 int CheckFinite() const;
607
608 /// Set the graph ownership flag (I and J arrays).
SetGraphOwner(bool ownij)609 void SetGraphOwner(bool ownij)
610 { I.SetHostPtrOwner(ownij); J.SetHostPtrOwner(ownij); }
611
612 /// Set the data ownership flag (A array).
SetDataOwner(bool owna)613 void SetDataOwner(bool owna) { A.SetHostPtrOwner(owna); }
614
615 /// Get the graph ownership flag (I and J arrays).
OwnsGraph() const616 bool OwnsGraph() const { return I.OwnsHostPtr() && J.OwnsHostPtr(); }
617
618 /// Get the data ownership flag (A array).
OwnsData() const619 bool OwnsData() const { return A.OwnsHostPtr(); }
620
621 /// Lose the ownership of the graph (I, J) and data (A) arrays.
LoseData()622 void LoseData() { SetGraphOwner(false); SetDataOwner(false); }
623
624 void Swap(SparseMatrix &other);
625
626 /// Destroys sparse matrix.
~SparseMatrix()627 virtual ~SparseMatrix()
628 {
629 Destroy();
630 #ifdef MFEM_USE_CUDA
631 if (useCuSparse)
632 {
633 if (SparseMatrixCount==1)
634 {
635 if (handle)
636 {
637 cusparseDestroy(handle);
638 handle = nullptr;
639 }
640 if (dBuffer)
641 {
642 CuMemFree(dBuffer);
643 dBuffer = nullptr;
644 bufferSize = 0;
645 }
646 }
647 SparseMatrixCount--;
648 }
649 #endif
650 }
651
GetType() const652 Type GetType() const { return MFEM_SPARSEMAT; }
653 };
654
operator <<(std::ostream & os,SparseMatrix const & mat)655 inline std::ostream& operator<<(std::ostream& os, SparseMatrix const& mat)
656 {
657 mat.Print(os);
658 return os;
659 }
660
661 /// Applies f() to each element of the matrix (after it is finalized).
662 void SparseMatrixFunction(SparseMatrix &S, double (*f)(double));
663
664
665 /// Transpose of a sparse matrix. A must be finalized.
666 SparseMatrix *Transpose(const SparseMatrix &A);
667 /// Transpose of a sparse matrix. A does not need to be a CSR matrix.
668 SparseMatrix *TransposeAbstractSparseMatrix (const AbstractSparseMatrix &A,
669 int useActualWidth);
670
671 /// Matrix product A.B.
672 /** If @a OAB is not NULL, we assume it has the structure of A.B and store the
673 result in @a OAB. If @a OAB is NULL, we create a new SparseMatrix to store
674 the result and return a pointer to it.
675
676 All matrices must be finalized. */
677 SparseMatrix *Mult(const SparseMatrix &A, const SparseMatrix &B,
678 SparseMatrix *OAB = NULL);
679
680 /// C = A^T B
681 SparseMatrix *TransposeMult(const SparseMatrix &A, const SparseMatrix &B);
682
683 /// Matrix product of sparse matrices. A and B do not need to be CSR matrices
684 SparseMatrix *MultAbstractSparseMatrix (const AbstractSparseMatrix &A,
685 const AbstractSparseMatrix &B);
686
687 /// Matrix product A.B
688 DenseMatrix *Mult(const SparseMatrix &A, DenseMatrix &B);
689
690 /// RAP matrix product (with R=P^T)
691 DenseMatrix *RAP(const SparseMatrix &A, DenseMatrix &P);
692
693 /// RAP matrix product (with R=P^T)
694 DenseMatrix *RAP(DenseMatrix &A, const SparseMatrix &P);
695
696 /** RAP matrix product (with P=R^T). ORAP is like OAB above.
697 All matrices must be finalized. */
698 SparseMatrix *RAP(const SparseMatrix &A, const SparseMatrix &R,
699 SparseMatrix *ORAP = NULL);
700
701 /// General RAP with given R^T, A and P
702 SparseMatrix *RAP(const SparseMatrix &Rt, const SparseMatrix &A,
703 const SparseMatrix &P);
704
705 /// Matrix multiplication A^t D A. All matrices must be finalized.
706 SparseMatrix *Mult_AtDA(const SparseMatrix &A, const Vector &D,
707 SparseMatrix *OAtDA = NULL);
708
709
710 /// Matrix addition result = A + B.
711 SparseMatrix * Add(const SparseMatrix & A, const SparseMatrix & B);
712 /// Matrix addition result = a*A + b*B
713 SparseMatrix * Add(double a, const SparseMatrix & A, double b,
714 const SparseMatrix & B);
715 /// Matrix addition result = sum_i A_i
716 SparseMatrix * Add(Array<SparseMatrix *> & Ai);
717
718 /// B += alpha * A
719 void Add(const SparseMatrix &A, double alpha, DenseMatrix &B);
720
721 /// Produces a block matrix with blocks A_{ij}*B
722 DenseMatrix *OuterProduct(const DenseMatrix &A, const DenseMatrix &B);
723
724 /// Produces a block matrix with blocks A_{ij}*B
725 SparseMatrix *OuterProduct(const DenseMatrix &A, const SparseMatrix &B);
726
727 /// Produces a block matrix with blocks A_{ij}*B
728 SparseMatrix *OuterProduct(const SparseMatrix &A, const DenseMatrix &B);
729
730 /// Produces a block matrix with blocks A_{ij}*B
731 SparseMatrix *OuterProduct(const SparseMatrix &A, const SparseMatrix &B);
732
733
734 // Inline methods
735
SetColPtr(const int row) const736 inline void SparseMatrix::SetColPtr(const int row) const
737 {
738 if (Rows)
739 {
740 if (ColPtrNode == NULL)
741 {
742 ColPtrNode = new RowNode *[width];
743 for (int i = 0; i < width; i++)
744 {
745 ColPtrNode[i] = NULL;
746 }
747 }
748 for (RowNode *node_p = Rows[row]; node_p != NULL; node_p = node_p->Prev)
749 {
750 ColPtrNode[node_p->Column] = node_p;
751 }
752 }
753 else
754 {
755 if (ColPtrJ == NULL)
756 {
757 ColPtrJ = new int[width];
758 for (int i = 0; i < width; i++)
759 {
760 ColPtrJ[i] = -1;
761 }
762 }
763 for (int j = I[row], end = I[row+1]; j < end; j++)
764 {
765 ColPtrJ[J[j]] = j;
766 }
767 }
768 current_row = row;
769 }
770
ClearColPtr() const771 inline void SparseMatrix::ClearColPtr() const
772 {
773 if (Rows)
774 {
775 for (RowNode *node_p = Rows[current_row]; node_p != NULL;
776 node_p = node_p->Prev)
777 {
778 ColPtrNode[node_p->Column] = NULL;
779 }
780 }
781 else
782 {
783 for (int j = I[current_row], end = I[current_row+1]; j < end; j++)
784 {
785 ColPtrJ[J[j]] = -1;
786 }
787 }
788 }
789
SearchRow(const int col)790 inline double &SparseMatrix::SearchRow(const int col)
791 {
792 if (Rows)
793 {
794 RowNode *node_p = ColPtrNode[col];
795 if (node_p == NULL)
796 {
797 #ifdef MFEM_USE_MEMALLOC
798 node_p = NodesMem->Alloc();
799 #else
800 node_p = new RowNode;
801 #endif
802 node_p->Prev = Rows[current_row];
803 node_p->Column = col;
804 node_p->Value = 0.0;
805 Rows[current_row] = ColPtrNode[col] = node_p;
806 }
807 return node_p->Value;
808 }
809 else
810 {
811 const int j = ColPtrJ[col];
812 MFEM_VERIFY(j != -1, "Entry for column " << col << " is not allocated.");
813 return A[j];
814 }
815 }
816
_Get_(const int col) const817 inline double SparseMatrix::_Get_(const int col) const
818 {
819 if (Rows)
820 {
821 RowNode *node_p = ColPtrNode[col];
822 return (node_p == NULL) ? 0.0 : node_p->Value;
823 }
824 else
825 {
826 const int j = ColPtrJ[col];
827 return (j == -1) ? 0.0 : A[j];
828 }
829 }
830
SearchRow(const int row,const int col)831 inline double &SparseMatrix::SearchRow(const int row, const int col)
832 {
833 if (Rows)
834 {
835 RowNode *node_p;
836
837 for (node_p = Rows[row]; 1; node_p = node_p->Prev)
838 {
839 if (node_p == NULL)
840 {
841 #ifdef MFEM_USE_MEMALLOC
842 node_p = NodesMem->Alloc();
843 #else
844 node_p = new RowNode;
845 #endif
846 node_p->Prev = Rows[row];
847 node_p->Column = col;
848 node_p->Value = 0.0;
849 Rows[row] = node_p;
850 break;
851 }
852 else if (node_p->Column == col)
853 {
854 break;
855 }
856 }
857 return node_p->Value;
858 }
859 else
860 {
861 int *Ip = I+row, *Jp = J;
862 for (int k = Ip[0], end = Ip[1]; k < end; k++)
863 {
864 if (Jp[k] == col)
865 {
866 return A[k];
867 }
868 }
869 MFEM_ABORT("Could not find entry for row = " << row << ", col = " << col);
870 }
871 return A[0];
872 }
873
874 /// Specialization of the template function Swap<> for class SparseMatrix
Swap(SparseMatrix & a,SparseMatrix & b)875 template<> inline void Swap<SparseMatrix>(SparseMatrix &a, SparseMatrix &b)
876 {
877 a.Swap(b);
878 }
879
880 } // namespace mfem
881
882 #endif
883