MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
sparsemat.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
34namespace mfem
35{
36
37class
38#if defined(__alignas_is_defined)
39 alignas(real_t)
40#endif
42{
43public:
46 int Column;
47};
48
49/// Data type sparse matrix
51{
52protected:
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. */
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
83 typedef MemAlloc <RowNode, 1024> RowNodeAlloc;
85#endif
86
87 /// Are the columns sorted already.
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
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
129public:
130 /// Create an empty SparseMatrix.
132 {
133 SetEmpty();
134
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, real_t *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, real_t *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
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 real_t *GetData() { return A; }
233 /// Return the element data, i.e. the array #A, const version.
234 inline const real_t *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); }
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); }
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<real_t> &GetMemoryData() const { return A; }
271 const real_t *ReadData(bool on_dev = true) const
272 { return mfem::Read(A, A.Capacity(), on_dev); }
273 real_t *WriteData(bool on_dev = true)
274 { return mfem::Write(A, A.Capacity(), on_dev); }
275 real_t *ReadWriteData(bool on_dev = true)
276 { return mfem::ReadWrite(A, A.Capacity(), on_dev); }
277 const real_t *HostReadData() const
278 { return mfem::Read(A, A.Capacity(), false); }
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 real_t *GetRowEntries(const int row);
297 /// Return a pointer to the entries in a row, const version.
298 const real_t *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 real_t &Elem(int i, int j);
323
324 /// Returns constant reference to a_{ij}.
325 virtual const real_t &Elem(int i, int j) const;
326
327 /// Returns reference to A[i][j].
328 real_t &operator()(int i, int j);
329
330 /// Returns reference to A[i][j].
331 const real_t &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
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 real_t 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 real_t 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 real_t 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 real_t InnerProduct(const Vector &x, const Vector &y) const;
435
436 /// For all i compute $ x_i = \sum_j A_{ij} $
437 void GetRowSums(Vector &x) const;
438 /// For i = irow compute $ x_i = \sum_j | A_{i, j} | $
439 real_t 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 real_t 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 real_t 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, real_t 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 real_t 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 real_t 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 real_t 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 real_t sc = 1.0, bool use_abs_diag = false) const;
532
533 /** x1 = x0 + sc D^{-1} (b - A x0) where $ D_{ii} = \sum_j |A_{ij}| $. */
534 void Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
535 real_t sc = 1.0) const;
536
537 /** x1 = x0 + sc D^{-1} (b - A x0) where $ D_{ii} = \sum_j A_{ij} $. */
538 void Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
539 real_t 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(real_t 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 real_t &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 real_t a)
586 { SearchRow(col) += a; }
587 /// Set an entry in the "current row". See SetColPtr().
588 inline void _Set_(const int col, const real_t a)
589 { SearchRow(col) = a; }
590 /// Read the value of an entry in the "current row". See SetColPtr().
591 inline real_t _Get_(const int col) const;
592
593 inline real_t &SearchRow(const int row, const int col);
594 inline void _Add_(const int row, const int col, const real_t a)
595 { SearchRow(row, col) += a; }
596 inline void _Set_(const int row, const int col, const real_t a)
597 { SearchRow(row, col) = a; }
598
599 void Set(const int i, const int j, const real_t val);
600 void Add(const int i, const int j, const real_t 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 real_t 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 real_t a, const SparseMatrix &B);
647
649
651
652 /// Prints matrix to stream out.
653 /** @note The host in synchronized when the finalized matrix is on the device. */
654 void Print(std::ostream &out = mfem::out, int width_ = 4) const;
655
656 /// Prints matrix in matlab format.
657 /** @note The host in synchronized when the finalized matrix is on the device. */
658 virtual void PrintMatlab(std::ostream &out = mfem::out) const;
659
660 /// Prints matrix in Matrix Market sparse format.
661 /** @note The host in synchronized when the finalized matrix is on the device. */
662 void PrintMM(std::ostream &out = mfem::out) const;
663
664 /// Prints matrix to stream out in hypre_CSRMatrix format.
665 /** @note The host in synchronized when the finalized matrix is on the device. */
666 void PrintCSR(std::ostream &out) const;
667
668 /// Prints a sparse matrix to stream out in CSR format.
669 /** @note The host in synchronized when the finalized matrix is on the device. */
670 void PrintCSR2(std::ostream &out) const;
671
672 /// Print various sparse matrix statistics.
673 void PrintInfo(std::ostream &out) const;
674
675 /// Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix
676 real_t IsSymmetric() const;
677
678 /// (*this) = 1/2 ((*this) + (*this)^t)
679 void Symmetrize();
680
681 /// Returns the number of the nonzero elements in the matrix
682 virtual int NumNonZeroElems() const;
683
684 real_t MaxNorm() const;
685
686 /// Count the number of entries with |a_ij| <= tol.
687 int CountSmallElems(real_t tol) const;
688
689 /// Count the number of entries that are NOT finite, i.e. Inf or Nan.
690 int CheckFinite() const;
691
692 /// Set the graph ownership flag (I and J arrays).
693 void SetGraphOwner(bool ownij)
694 { I.SetHostPtrOwner(ownij); J.SetHostPtrOwner(ownij); }
695
696 /// Set the data ownership flag (A array).
697 void SetDataOwner(bool owna) { A.SetHostPtrOwner(owna); }
698
699 /// Get the graph ownership flag (I and J arrays).
700 bool OwnsGraph() const { return I.OwnsHostPtr() && J.OwnsHostPtr(); }
701
702 /// Get the data ownership flag (A array).
703 bool OwnsData() const { return A.OwnsHostPtr(); }
704
705 /// Lose the ownership of the graph (I, J) and data (A) arrays.
706 void LoseData() { SetGraphOwner(false); SetDataOwner(false); }
707
708 void Swap(SparseMatrix &other);
709
710 /// Destroys sparse matrix.
711 virtual ~SparseMatrix();
712
713 Type GetType() const { return MFEM_SPARSEMAT; }
714};
715
716inline std::ostream& operator<<(std::ostream& os, SparseMatrix const& mat)
717{
718 mat.Print(os);
719 return os;
720}
721
722/// Applies f() to each element of the matrix (after it is finalized).
723void SparseMatrixFunction(SparseMatrix &S, real_t (*f)(real_t));
724
725
726/// Transpose of a sparse matrix. A must be finalized.
727SparseMatrix *Transpose(const SparseMatrix &A);
728/// Transpose of a sparse matrix. A does not need to be a CSR matrix.
729SparseMatrix *TransposeAbstractSparseMatrix (const AbstractSparseMatrix &A,
730 int useActualWidth);
731
732/// Matrix product A.B.
733/** If @a OAB is not NULL, we assume it has the structure of A.B and store the
734 result in @a OAB. If @a OAB is NULL, we create a new SparseMatrix to store
735 the result and return a pointer to it.
736
737 All matrices must be finalized. */
738SparseMatrix *Mult(const SparseMatrix &A, const SparseMatrix &B,
739 SparseMatrix *OAB = NULL);
740
741/// C = A^T B
742SparseMatrix *TransposeMult(const SparseMatrix &A, const SparseMatrix &B);
743
744/// Matrix product of sparse matrices. A and B do not need to be CSR matrices
745SparseMatrix *MultAbstractSparseMatrix (const AbstractSparseMatrix &A,
746 const AbstractSparseMatrix &B);
747
748/// Matrix product A.B
749DenseMatrix *Mult(const SparseMatrix &A, DenseMatrix &B);
750
751/// RAP matrix product (with R=P^T)
752DenseMatrix *RAP(const SparseMatrix &A, DenseMatrix &P);
753
754/// RAP matrix product (with R=P^T)
755DenseMatrix *RAP(DenseMatrix &A, const SparseMatrix &P);
756
757/** RAP matrix product (with P=R^T). ORAP is like OAB above.
758 All matrices must be finalized. */
759SparseMatrix *RAP(const SparseMatrix &A, const SparseMatrix &R,
760 SparseMatrix *ORAP = NULL);
761
762/// General RAP with given R^T, A and P
763SparseMatrix *RAP(const SparseMatrix &Rt, const SparseMatrix &A,
764 const SparseMatrix &P);
765
766/// Matrix multiplication A^t D A. All matrices must be finalized.
767SparseMatrix *Mult_AtDA(const SparseMatrix &A, const Vector &D,
768 SparseMatrix *OAtDA = NULL);
769
770
771/// Matrix addition result = A + B.
772SparseMatrix * Add(const SparseMatrix & A, const SparseMatrix & B);
773/// Matrix addition result = a*A + b*B
774SparseMatrix * Add(real_t a, const SparseMatrix & A, real_t b,
775 const SparseMatrix & B);
776/// Matrix addition result = sum_i A_i
777SparseMatrix * Add(Array<SparseMatrix *> & Ai);
778
779/// B += alpha * A
780void Add(const SparseMatrix &A, real_t alpha, DenseMatrix &B);
781
782/// Produces a block matrix with blocks A_{ij}*B
783DenseMatrix *OuterProduct(const DenseMatrix &A, const DenseMatrix &B);
784
785/// Produces a block matrix with blocks A_{ij}*B
786SparseMatrix *OuterProduct(const DenseMatrix &A, const SparseMatrix &B);
787
788/// Produces a block matrix with blocks A_{ij}*B
789SparseMatrix *OuterProduct(const SparseMatrix &A, const DenseMatrix &B);
790
791/// Produces a block matrix with blocks A_{ij}*B
792SparseMatrix *OuterProduct(const SparseMatrix &A, const SparseMatrix &B);
793
794
795// Inline methods
796
797inline void SparseMatrix::SetColPtr(const int row) const
798{
799 if (Rows)
800 {
801 if (ColPtrNode == NULL)
802 {
803 ColPtrNode = new RowNode *[width];
804 for (int i = 0; i < width; i++)
805 {
806 ColPtrNode[i] = NULL;
807 }
808 }
809 for (RowNode *node_p = Rows[row]; node_p != NULL; node_p = node_p->Prev)
810 {
811 ColPtrNode[node_p->Column] = node_p;
812 }
813 }
814 else
815 {
816 if (ColPtrJ == NULL)
817 {
818 ColPtrJ = new int[width];
819 for (int i = 0; i < width; i++)
820 {
821 ColPtrJ[i] = -1;
822 }
823 }
824 for (int j = I[row], end = I[row+1]; j < end; j++)
825 {
826 ColPtrJ[J[j]] = j;
827 }
828 }
829 current_row = row;
830}
831
832inline void SparseMatrix::ClearColPtr() const
833{
834 if (Rows)
835 {
836 for (RowNode *node_p = Rows[current_row]; node_p != NULL;
837 node_p = node_p->Prev)
838 {
839 ColPtrNode[node_p->Column] = NULL;
840 }
841 }
842 else
843 {
844 for (int j = I[current_row], end = I[current_row+1]; j < end; j++)
845 {
846 ColPtrJ[J[j]] = -1;
847 }
848 }
849}
850
851inline real_t &SparseMatrix::SearchRow(const int col)
852{
853 if (Rows)
854 {
855 RowNode *node_p = ColPtrNode[col];
856 if (node_p == NULL)
857 {
858#ifdef MFEM_USE_MEMALLOC
859 node_p = NodesMem->Alloc();
860#else
861 node_p = new RowNode;
862#endif
863 node_p->Prev = Rows[current_row];
864 node_p->Column = col;
865 node_p->Value = 0.0;
866 Rows[current_row] = ColPtrNode[col] = node_p;
867 }
868 return node_p->Value;
869 }
870 else
871 {
872 const int j = ColPtrJ[col];
873 MFEM_VERIFY(j != -1, "Entry for column " << col << " is not allocated.");
874 return A[j];
875 }
876}
877
878inline real_t SparseMatrix::_Get_(const int col) const
879{
880 if (Rows)
881 {
882 RowNode *node_p = ColPtrNode[col];
883 return (node_p == NULL) ? 0.0 : node_p->Value;
884 }
885 else
886 {
887 const int j = ColPtrJ[col];
888 return (j == -1) ? 0.0 : A[j];
889 }
890}
891
892inline real_t &SparseMatrix::SearchRow(const int row, const int col)
893{
894 if (Rows)
895 {
896 RowNode *node_p;
897
898 for (node_p = Rows[row]; 1; node_p = node_p->Prev)
899 {
900 if (node_p == NULL)
901 {
902#ifdef MFEM_USE_MEMALLOC
903 node_p = NodesMem->Alloc();
904#else
905 node_p = new RowNode;
906#endif
907 node_p->Prev = Rows[row];
908 node_p->Column = col;
909 node_p->Value = 0.0;
910 Rows[row] = node_p;
911 break;
912 }
913 else if (node_p->Column == col)
914 {
915 break;
916 }
917 }
918 return node_p->Value;
919 }
920 else
921 {
922 int *Ip = I+row, *Jp = J;
923 for (int k = Ip[0], end = Ip[1]; k < end; k++)
924 {
925 if (Jp[k] == col)
926 {
927 return A[k];
928 }
929 }
930 MFEM_ABORT("Could not find entry for row = " << row << ", col = " << col);
931 }
932 return A[0];
933}
934
935/// Specialization of the template function Swap<> for class SparseMatrix
937{
938 a.Swap(b);
939}
940
941} // namespace mfem
942
943#endif
Abstract data type for sparse matrices.
Definition: matrix.hpp:74
Dynamic 2D array using row-major layout.
Definition: array.hpp:372
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
static MemoryClass GetHostMemoryClass()
Get the current Host MemoryClass. This is the MemoryClass used by most MFEM host Memory objects.
Definition: device.hpp:269
static MemoryClass GetDeviceMemoryClass()
Get the current Device MemoryClass. This is the MemoryClass used by most MFEM device kernels to acces...
Definition: device.hpp:282
Abstract data type for matrix inverse.
Definition: matrix.hpp:63
Elem * Alloc()
Definition: mem_alloc.hpp:166
Class used by MFEM to store pointers to host and/or device memory.
void SetHostPtrOwner(bool own) const
Set/clear the ownership flag for the host pointer. Ownership indicates whether the pointer will be de...
int Capacity() const
Return the size of the allocated memory.
bool OwnsHostPtr() const
Return true if the host pointer is owned. Ownership indicates whether the pointer will be deleted by ...
bool Empty() const
Return true if the Memory object is empty, see Reset().
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:48
@ DIAG_ONE
Set the diagonal value to one.
Definition: operator.hpp:50
@ DIAG_ZERO
Set the diagonal value to zero.
Definition: operator.hpp:49
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:284
@ MFEM_SPARSEMAT
ID for class SparseMatrix.
Definition: operator.hpp:286
RowNode * Prev
Definition: sparsemat.hpp:45
real_t Value
Definition: sparsemat.hpp:44
Data type sparse matrix.
Definition: sparsemat.hpp:51
void PrintCSR2(std::ostream &out) const
Prints a sparse matrix to stream out in CSR format.
Definition: sparsemat.cpp:3457
real_t MaxNorm() const
Definition: sparsemat.cpp:1665
const int * HostReadJ() const
Definition: sparsemat.hpp:261
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:691
int * ReadWriteI(bool on_dev=true)
Definition: sparsemat.hpp:243
real_t GetRowNorml1(int irow) const
For i = irow compute .
Definition: sparsemat.cpp:1291
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:954
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Definition: sparsemat.cpp:1036
bool isSorted
Are the columns sorted already.
Definition: sparsemat.hpp:88
RowNodeAlloc * NodesMem
Definition: sparsemat.hpp:84
bool ColumnsAreSorted() const
Returns whether or not the columns are sorted.
Definition: sparsemat.hpp:554
void Gauss_Seidel_back(const Vector &x, Vector &y) const
Definition: sparsemat.cpp:2527
bool RowIsEmpty(const int row) const
Definition: sparsemat.cpp:2969
void ClearColPtr() const
Reset the "current row" set by calling SetColPtr(). This method must be called between any two calls ...
Definition: sparsemat.hpp:832
void Jacobi3(const Vector &b, const Vector &x0, Vector &x1, real_t sc=1.0) const
Definition: sparsemat.cpp:2770
real_t * ReadWriteData(bool on_dev=true)
Definition: sparsemat.hpp:275
Memory< real_t > A
Array with size I[height], containing the actual entries of the sparse matrix, as indexed by the I ar...
Definition: sparsemat.hpp:68
void EliminateRowCol(int rc, const real_t sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate row rc and column rc and modify the rhs using sol.
Definition: sparsemat.cpp:1909
void SetDataOwner(bool owna)
Set the data ownership flag (A array).
Definition: sparsemat.hpp:697
SparseMatrix & operator*=(real_t a)
Definition: sparsemat.cpp:3275
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:2990
hipsparseDnVecDescr_t vecY_descr
Definition: sparsemat.hpp:125
void SetWidth(int width_=-1)
Change the width of a SparseMatrix.
Definition: sparsemat.cpp:419
int * HostReadWriteI()
Definition: sparsemat.hpp:249
MemAlloc< RowNode, 1024 > RowNodeAlloc
Definition: sparsemat.hpp:83
bool Empty() const
Check if the SparseMatrix is empty.
Definition: sparsemat.hpp:219
static hipsparseHandle_t handle
Definition: sparsemat.hpp:120
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
Definition: sparsemat.cpp:2941
SparseMatrix * At
Transpose of A. Owned. Used to perform MultTranspose() on devices.
Definition: sparsemat.hpp:80
virtual MatrixInverse * Inverse() const
This virtual method is not supported: it always returns NULL.
Definition: sparsemat.cpp:1742
void PartAddMult(const Array< int > &rows, const Vector &x, Vector &y, const real_t a=1.0) const
Definition: sparsemat.cpp:1062
const int * GetJ() const
Return the array J, const version.
Definition: sparsemat.hpp:229
void GetRowSums(Vector &x) const
For all i compute .
Definition: sparsemat.cpp:1268
cusparseMatDescr_t descr
Definition: sparsemat.hpp:108
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:552
Type GetType() const
Definition: sparsemat.hpp:713
void EliminateBC(const Array< int > &ess_dofs, DiagonalPolicy diag_policy)
Eliminate essential (Dirichlet) boundary conditions.
Definition: sparsemat.cpp:2370
void PrintInfo(std::ostream &out) const
Print various sparse matrix statistics.
Definition: sparsemat.cpp:3486
virtual void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const
y += At * x (default) or y += a * At * x
Definition: sparsemat.cpp:961
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:599
void Add(const int i, const int j, const real_t val)
Definition: sparsemat.cpp:2842
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:2008
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, real_t sc, bool use_abs_diag=false) const
Definition: sparsemat.cpp:2643
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:797
void SetDiagIdentity()
If a row contains only one diag entry of zero, set it to 1.
Definition: sparsemat.cpp:2413
MFEM_DEPRECATED void UseCuSparse(bool useCuSparse_=true)
Deprecated equivalent of UseGPUSparse().
Definition: sparsemat.hpp:194
const Memory< int > & GetMemoryJ() const
Definition: sparsemat.hpp:254
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:1182
void Gauss_Seidel_forw(const Vector &x, Vector &y) const
Gauss-Seidel forward and backward iterations over a vector x.
Definition: sparsemat.cpp:2443
void MakeRef(const SparseMatrix &master)
Clear the contents of the SparseMatrix and make it a reference to master.
Definition: sparsemat.cpp:313
int CheckFinite() const
Count the number of entries that are NOT finite, i.e. Inf or Nan.
Definition: sparsemat.cpp:1718
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
int * WriteJ(bool on_dev=true)
Definition: sparsemat.hpp:257
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
cusparseMatDescr_t matA_descr
Definition: sparsemat.hpp:115
void EliminateCol(int col, DiagonalPolicy dpolicy=DIAG_ZERO)
Eliminates the column col from the matrix.
Definition: sparsemat.cpp:1795
void LoseData()
Lose the ownership of the graph (I, J) and data (A) arrays.
Definition: sparsemat.hpp:706
void EnsureMultTranspose() const
Ensures that the matrix is capable of performing MultTranspose(), AddMultTranspose(),...
Definition: sparsemat.cpp:1028
void _Add_(const int col, const real_t a)
Add a value to an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:585
SparseMatrix()
Create an empty SparseMatrix.
Definition: sparsemat.hpp:131
cusparseSpMatDescr_t matA_descr
Definition: sparsemat.hpp:111
int * HostReadWriteJ()
Definition: sparsemat.hpp:265
void AbsMult(const Vector &x, Vector &y) const
y = |A| * x, using entry-wise absolute values of matrix A
Definition: sparsemat.cpp:1133
void ClearGPUSparse()
Clear the cuSPARSE/hipSPARSE descriptors. This must be called after releasing the device memory of A.
Definition: sparsemat.cpp:81
void Threshold(real_t 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:1315
Memory< int > & GetMemoryI()
Definition: sparsemat.hpp:237
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2861
const int * ReadI(bool on_dev=true) const
Definition: sparsemat.hpp:239
real_t InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Definition: sparsemat.cpp:1227
Memory< int > I
Array with size (height+1) containing the row offsets.
Definition: sparsemat.hpp:62
real_t & operator()(int i, int j)
Returns reference to A[i][j].
Definition: sparsemat.cpp:633
const real_t * HostReadData() const
Definition: sparsemat.hpp:277
void PrintMM(std::ostream &out=mfem::out) const
Prints matrix in Matrix Market sparse format.
Definition: sparsemat.cpp:3388
void ScaleColumns(const Vector &sr)
this = this * diag(sr);
Definition: sparsemat.cpp:3168
void _Add_(const int row, const int col, const real_t a)
Definition: sparsemat.hpp:594
virtual void PrintMatlab(std::ostream &out=mfem::out) const
Prints matrix in matlab format.
Definition: sparsemat.cpp:3347
void PrintCSR(std::ostream &out) const
Prints matrix to stream out in hypre_CSRMatrix format.
Definition: sparsemat.cpp:3429
void Swap(SparseMatrix &other)
Definition: sparsemat.cpp:4283
MFEM_DEPRECATED void ClearCuSparse()
Deprecated equivalent of ClearGPUSparse().
Definition: sparsemat.hpp:216
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: sparsemat.hpp:342
const Memory< real_t > & GetMemoryData() const
Definition: sparsemat.hpp:270
int * ReadWriteJ(bool on_dev=true)
Definition: sparsemat.hpp:259
cusparseDnVecDescr_t vecY_descr
Definition: sparsemat.hpp:113
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:1110
static size_t bufferSize
Definition: sparsemat.hpp:101
RowNode ** Rows
Array of linked lists, one for every row. This array represents the linked list (LIL) storage format.
Definition: sparsemat.hpp:73
void SetGraphOwner(bool ownij)
Set the graph ownership flag (I and J arrays).
Definition: sparsemat.hpp:693
real_t & SearchRow(const int col)
Perform a fast search for an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:851
int * WriteI(bool on_dev=true)
Definition: sparsemat.hpp:241
const Memory< int > & GetMemoryI() const
Definition: sparsemat.hpp:238
static int SparseMatrixCount
Definition: sparsemat.hpp:100
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2777
void _Set_(const int row, const int col, const real_t a)
Definition: sparsemat.hpp:596
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
Definition: sparsemat.cpp:1623
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:1080
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:391
real_t GetJacobiScaling() const
Determine appropriate scaling for Jacobi iteration.
Definition: sparsemat.cpp:2610
RowNode ** ColPtrNode
Definition: sparsemat.hpp:77
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:344
real_t * HostReadWriteData()
Definition: sparsemat.hpp:281
int CountSmallElems(real_t tol) const
Count the number of entries with |a_ij| <= tol.
Definition: sparsemat.cpp:1690
const int * ReadJ(bool on_dev=true) const
Definition: sparsemat.hpp:255
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:755
virtual real_t & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: sparsemat.cpp:623
int MaxRowSize() const
Returns the maximum number of elements among all rows.
Definition: sparsemat.cpp:367
real_t _Get_(const int col) const
Read the value of an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:878
virtual ~SparseMatrix()
Destroys sparse matrix.
Definition: sparsemat.cpp:4303
SparseMatrix & operator=(const SparseMatrix &rhs)
Assignment operator: deep copy.
Definition: sparsemat.cpp:303
const real_t * GetData() const
Return the element data, i.e. the array A, const version.
Definition: sparsemat.hpp:234
virtual void EliminateZeroRows(const real_t threshold=1e-12)
If a row contains only zeros, set its diagonal to 1.
Definition: sparsemat.cpp:2424
void OverrideSize(int height_, int width_)
Sets the height and width of the matrix.
Definition: sparsemat.cpp:297
real_t * HostWriteData()
Definition: sparsemat.hpp:279
cusparseStatus_t status
Definition: sparsemat.hpp:106
static cusparseHandle_t handle
Definition: sparsemat.hpp:107
void EliminateRow(int row, const real_t sol, Vector &rhs)
Eliminates a column from the transpose matrix.
Definition: sparsemat.cpp:1747
bool OwnsGraph() const
Get the graph ownership flag (I and J arrays).
Definition: sparsemat.hpp:700
hipsparseSpMatDescr_t matA_descr
Definition: sparsemat.hpp:123
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Definition: sparsemat.cpp:1491
SparseMatrix & operator+=(const SparseMatrix &B)
Add the sparse matrix 'B' to '*this'. This operation will cause an error if '*this' is finalized and ...
Definition: sparsemat.cpp:3196
hipsparseStatus_t status
Definition: sparsemat.hpp:119
void Clear()
Clear the contents of the SparseMatrix.
Definition: sparsemat.hpp:209
void ResetTranspose() const
Definition: sparsemat.cpp:1022
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR.
Definition: sparsemat.hpp:546
void EliminateRowColDiag(int rc, real_t value)
Perform elimination and set the diagonal entry to the given value.
Definition: sparsemat.cpp:2213
real_t * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:405
void ScaleRows(const Vector &sl)
this = diag(sl) * this;
Definition: sparsemat.cpp:3137
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:3030
bool OwnsData() const
Get the data ownership flag (A array).
Definition: sparsemat.hpp:703
const int * GetI() const
Return the array I, const version.
Definition: sparsemat.hpp:224
virtual void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const
y += A * x (default) or y += a * A * x
Definition: sparsemat.cpp:762
void DiagScale(const Vector &b, Vector &x, real_t 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:2675
real_t * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:232
void SetSubMatrixTranspose(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2900
void BuildTranspose() const
Build and store internally the transpose of this matrix which will be used in the methods AddMultTran...
Definition: sparsemat.cpp:1014
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:457
real_t IsSymmetric() const
Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix.
Definition: sparsemat.cpp:1581
DenseMatrix * ToDenseMatrix() const
Produces a DenseMatrix from a SparseMatrix.
Definition: sparsemat.cpp:725
void _Set_(const int col, const real_t a)
Set an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:588
static void * dBuffer
Definition: sparsemat.hpp:102
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:1836
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:227
Memory< int > & GetMemoryJ()
Definition: sparsemat.hpp:253
Memory< real_t > & GetMemoryData()
Definition: sparsemat.hpp:269
cusparseDnVecDescr_t vecX_descr
Definition: sparsemat.hpp:112
const real_t * ReadData(bool on_dev=true) const
Definition: sparsemat.hpp:271
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:1642
void ScaleRow(const int row, const real_t scale)
Definition: sparsemat.cpp:3109
real_t * WriteData(bool on_dev=true)
Definition: sparsemat.hpp:273
hipsparseDnVecDescr_t vecX_descr
Definition: sparsemat.hpp:124
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:3077
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1, real_t sc=1.0) const
Definition: sparsemat.cpp:2763
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
int ActualWidth() const
Returns the actual Width of the matrix.
Definition: sparsemat.cpp:3570
const int * HostReadI() const
Definition: sparsemat.hpp:245
void Set(const int i, const int j, const real_t val)
Definition: sparsemat.cpp:2823
int Size() const
For backward compatibility, define Size() to be synonym of Height().
Definition: sparsemat.hpp:206
void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: sparsemat.cpp:3299
Vector data type.
Definition: vector.hpp:80
const real_t alpha
Definition: ex15.cpp:369
real_t b
Definition: lissajous.cpp:42
real_t a
Definition: lissajous.cpp:41
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
Definition: sparsemat.hpp:716
DenseMatrix * OuterProduct(const DenseMatrix &A, const DenseMatrix &B)
Produces a block matrix with blocks A_{ij}*B.
Definition: sparsemat.cpp:4175
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's DeviceMemoryClass, if on_dev = true,...
Definition: device.hpp:320
void Swap< SparseMatrix >(SparseMatrix &a, SparseMatrix &b)
Specialization of the template function Swap<> for class SparseMatrix.
Definition: sparsemat.hpp:936
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
T * Write(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for write access to mem with the mfem::Device's DeviceMemoryClass, if on_dev = true,...
Definition: device.hpp:337
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
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:414
MemoryClass
Memory classes identify sets of memory types.
Definition: mem_manager.hpp:81
T * ReadWrite(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read+write access to mem with the mfem::Device's DeviceMemoryClass,...
Definition: device.hpp:354
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3464
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:3876
SparseMatrix * TransposeAbstractSparseMatrix(const AbstractSparseMatrix &A, int useActualWidth)
Transpose of a sparse matrix. A does not need to be a CSR matrix.
Definition: sparsemat.cpp:3663
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:4031
float real_t
Definition: config.hpp:43
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:39
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition: lor_mms.hpp:30
SparseMatrix * TransposeMult(const SparseMatrix &A, const SparseMatrix &B)
C = A^T B.
Definition: sparsemat.cpp:3868
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2433
void SparseMatrixFunction(SparseMatrix &S, real_t(*f)(real_t))
Applies f() to each element of the matrix (after it is finalized).
Definition: sparsemat.cpp:3596