MFEM  v3.3.2
Finite element discretization library
sparsemat.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #ifndef MFEM_SPARSEMAT_HPP
13 #define MFEM_SPARSEMAT_HPP
14 
15 // Data types for sparse matrix
16 
17 #include "../general/mem_alloc.hpp"
18 #include "../general/table.hpp"
19 #include "../general/globals.hpp"
20 #include "densemat.hpp"
21 
22 namespace mfem
23 {
24 
25 class
26 #if defined(__alignas_is_defined)
27  alignas(double)
28 #endif
29  RowNode
30 {
31 public:
32  double Value;
33  RowNode *Prev;
34  int Column;
35 };
36 
37 /// Data type sparse matrix
39 {
40 protected:
41  /// @name Arrays used by the CSR storage format.
42  /** */
43  ///@{
44  /// @brief %Array with size (#height+1) containing the row offsets.
45  /** The data for row r, 0 <= r < height, is at offsets j, I[r] <= j < I[r+1].
46  The offsets, j, are indices in the #J and #A arrays. The first entry in
47  this array is always zero, I[0] = 0, and the last entry, I[height], gives
48  the total number of entries stored (at a minimum, all nonzeros must be
49  represented) in the sparse matrix. */
50  int *I;
51  /** @brief %Array with size #I[#height], containing the column indices for
52  all matrix entries, as indexed by the #I array. */
53  int *J;
54  /** @brief %Array with size #I[#height], containing the actual entries of the
55  sparse matrix, as indexed by the #I array. */
56  double *A;
57  ///@}
58 
59  /** @brief %Array of linked lists, one for every row. This array represents
60  the linked list (LIL) storage format. */
61  RowNode **Rows;
62 
63  mutable int current_row;
64  mutable int* ColPtrJ;
65  mutable RowNode ** ColPtrNode;
66 
67 #ifdef MFEM_USE_MEMALLOC
69  RowNodeAlloc * NodesMem;
70 #endif
71 
72  /// Say whether we own the pointers for I and J (should we free them?).
73  bool ownGraph;
74  /// Say whether we own the pointers for A (should we free them?).
75  bool ownData;
76 
77  /// Are the columns sorted already.
78  bool isSorted;
79 
80  void Destroy(); // Delete all owned data
81  void SetEmpty(); // Init all entries with empty values
82 
83 public:
84  /// Create an empty SparseMatrix.
86 
87  /** @brief Create a sparse matrix with flexible sparsity structure using a
88  row-wise linked list (LIL) format. */
89  /** New entries are added as needed by methods like AddSubMatrix(),
90  SetSubMatrix(), etc. Calling Finalize() will convert the SparseMatrix to
91  the more compact compressed sparse row (CSR) format. */
92  explicit SparseMatrix(int nrows, int ncols = -1);
93 
94  /** @brief Create a sparse matrix in CSR format. Ownership of @a i, @a j, and
95  @a data is transferred to the SparseMatrix. */
96  SparseMatrix(int *i, int *j, double *data, int m, int n);
97 
98  /** @brief Create a sparse matrix in CSR format. Ownership of @a i, @a j, and
99  @a data is optionally transferred to the SparseMatrix. */
100  SparseMatrix(int *i, int *j, double *data, int m, int n, bool ownij,
101  bool owna, bool issorted);
102 
103  /** @brief Create a sparse matrix in CSR format where each row has space
104  allocated for exactly @a rowsize entries. */
105  /** SetRow() can then be called or the #I, #J, #A arrays can be used
106  directly. */
107  SparseMatrix(int nrows, int ncols, int rowsize);
108 
109  /// Copy constructor (deep copy).
110  /** If @a mat is finalized and @a copy_graph is false, the #I and #J arrays
111  will use a shallow copy (copy the pointers only) without transferring
112  ownership. */
113  SparseMatrix(const SparseMatrix &mat, bool copy_graph = true);
114 
115  /** @brief Clear the contents of the SparseMatrix and make it a reference to
116  @a master */
117  /** After this call, the matrix will point to the same data as @a master but
118  it will not own its data. The @a master must be finalized. */
119  void MakeRef(const SparseMatrix &master);
120 
121  /// For backward compatibility, define Size() to be synonym of Height().
122  int Size() const { return Height(); }
123 
124  /// Clear the contents of the SparseMatrix.
125  void Clear() { Destroy(); SetEmpty(); }
126 
127  /// Check if the SparseMatrix is empty.
128  bool Empty() const { return (A == NULL) && (Rows == NULL); }
129 
130  /// Return the array #I
131  inline int *GetI() const { return I; }
132  /// Return the array #J
133  inline int *GetJ() const { return J; }
134  /// Return element data, i.e. array #A
135  inline double *GetData() const { return A; }
136  /// Returns the number of elements in row @a i
137  int RowSize(const int i) const;
138  /// Returns the maximum number of elements among all rows
139  int MaxRowSize() const;
140  /// Return a pointer to the column indices in a row
141  int *GetRowColumns(const int row);
142  const int *GetRowColumns(const int row) const;
143  /// Return a pointer to the entries in a row
144  double *GetRowEntries(const int row);
145  const double *GetRowEntries(const int row) const;
146 
147  /// Change the width of a SparseMatrix.
148  /*!
149  * If width_ = -1 (DEFAULT), this routine will set the new width
150  * to the actual Width of the matrix awidth = max(J) + 1.
151  * Values 0 <= width_ < awidth are not allowed (error check in Debug Mode only)
152  *
153  * This method can be called for matrices finalized or not.
154  */
155  void SetWidth(int width_ = -1);
156 
157  /// Returns the actual Width of the matrix.
158  /*! This method can be called for matrices finalized or not. */
159  int ActualWidth();
160 
161  /// Sort the column indices corresponding to each row.
162  void SortColumnIndices();
163 
164  /** @brief Move the diagonal entry to the first position in each row,
165  preserving the order of the rest of the columns. */
166  void MoveDiagonalFirst();
167 
168  /// Returns reference to a_{ij}.
169  virtual double &Elem(int i, int j);
170 
171  /// Returns constant reference to a_{ij}.
172  virtual const double &Elem(int i, int j) const;
173 
174  /// Returns reference to A[i][j].
175  double &operator()(int i, int j);
176 
177  /// Returns reference to A[i][j].
178  const double &operator()(int i, int j) const;
179 
180  /// Returns the Diagonal of A
181  void GetDiag(Vector & d) const;
182 
183  /// Matrix vector multiplication.
184  virtual void Mult(const Vector &x, Vector &y) const;
185 
186  /// y += A * x (default) or y += a * A * x
187  void AddMult(const Vector &x, Vector &y, const double a = 1.0) const;
188 
189  /// Multiply a vector with the transposed matrix. y = At * x
190  void MultTranspose(const Vector &x, Vector &y) const;
191 
192  /// y += At * x (default) or y += a * At * x
193  void AddMultTranspose(const Vector &x, Vector &y,
194  const double a = 1.0) const;
195 
196  void PartMult(const Array<int> &rows, const Vector &x, Vector &y) const;
197  void PartAddMult(const Array<int> &rows, const Vector &x, Vector &y,
198  const double a=1.0) const;
199 
200  /// y = A * x, but treat all elements as booleans (zero=false, nonzero=true).
201  void BooleanMult(const Array<int> &x, Array<int> &y) const;
202  /// y = At * x, but treat all elements as booleans (zero=false, nonzero=true).
203  void BooleanMultTranspose(const Array<int> &x, Array<int> &y) const;
204 
205  /// Compute y^t A x
206  double InnerProduct(const Vector &x, const Vector &y) const;
207 
208  /// For all i compute \f$ x_i = \sum_j A_{ij} \f$
209  void GetRowSums(Vector &x) const;
210  /// For i = irow compute \f$ x_i = \sum_j | A_{i, j} | \f$
211  double GetRowNorml1(int irow) const;
212 
213  /// Returns a pointer to approximation of the matrix inverse.
214  virtual MatrixInverse *Inverse() const;
215 
216  /// Eliminates a column from the transpose matrix.
217  void EliminateRow(int row, const double sol, Vector &rhs);
218  /// Eliminates a row from the matrix.
219  /*!
220  * If setOneDiagonal = 0, all the entries in the row will be set to 0.
221  * If setOneDiagonal = 1 (matrix must be square),
222  * the diagonal entry will be set equal to 1
223  * and all the others entries to 0.
224  */
225  void EliminateRow(int row, int setOneDiagonal = 0);
226  void EliminateCol(int col);
227  /// Eliminate all columns 'i' for which cols[i] != 0
228  void EliminateCols(Array<int> &cols, Vector *x = NULL, Vector *b = NULL);
229 
230  /** Eliminates the column 'rc' to the 'rhs', deletes the row 'rc' and
231  replaces the element (rc,rc) with 1.0; assumes that element (i,rc)
232  is assembled if and only if the element (rc,i) is assembled.
233  If d != 0 then the element (rc,rc) remains the same. */
234  void EliminateRowCol(int rc, const double sol, Vector &rhs, int d = 0);
235 
236  /** Like previous one, but multiple values for eliminated unknowns are
237  accepted, and accordingly multiple right-hand-sides are used. */
238  void EliminateRowColMultipleRHS(int rc, const Vector &sol,
239  DenseMatrix &rhs, int d = 0);
240 
241  void EliminateRowCol(int rc, int d = 0);
242  /// Perform elimination and set the diagonal entry to the given value
243  void EliminateRowColDiag(int rc, double value);
244  // Same as above + save the eliminated entries in Ae so that
245  // (*this) + Ae is the original matrix
246  void EliminateRowCol(int rc, SparseMatrix &Ae, int d = 0);
247 
248  /// If a row contains only one diag entry of zero, set it to 1.
249  void SetDiagIdentity();
250  /// If a row contains only zeros, set its diagonal to 1.
251  void EliminateZeroRows();
252 
253  /// Gauss-Seidel forward and backward iterations over a vector x.
254  void Gauss_Seidel_forw(const Vector &x, Vector &y) const;
255  void Gauss_Seidel_back(const Vector &x, Vector &y) const;
256 
257  /// Determine appropriate scaling for Jacobi iteration
258  double GetJacobiScaling() const;
259  /** One scaled Jacobi iteration for the system A x = b.
260  x1 = x0 + sc D^{-1} (b - A x0) where D is the diag of A. */
261  void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const;
262 
263  void DiagScale(const Vector &b, Vector &x, double sc = 1.0) const;
264 
265  /** x1 = x0 + sc D^{-1} (b - A x0) where \f$ D_{ii} = \sum_j |A_{ij}| \f$. */
266  void Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
267  double sc = 1.0) const;
268 
269  /** x1 = x0 + sc D^{-1} (b - A x0) where \f$ D_{ii} = \sum_j A_{ij} \f$. */
270  void Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
271  double sc = 1.0) const;
272 
273  /** @brief Finalize the matrix initialization, switching the storage format
274  from LIL to CSR. */
275  /** This method should be called once, after the matrix has been initialized.
276  Internally, this method converts the matrix from row-wise linked list
277  (LIL) format into CSR (compressed sparse row) format. */
278  virtual void Finalize(int skip_zeros = 1) { Finalize(skip_zeros, false); }
279 
280  /// A slightly more general version of the Finalize(int) method.
281  void Finalize(int skip_zeros, bool fix_empty_rows);
282 
283  bool Finalized() const { return (A != NULL); }
284  bool areColumnsSorted() const { return isSorted; }
285 
286  /** Split the matrix into M x N blocks of sparse matrices in CSR format.
287  The 'blocks' array is M x N (i.e. M and N are determined by its
288  dimensions) and its entries are overwritten by the new blocks. */
289  void GetBlocks(Array2D<SparseMatrix *> &blocks) const;
290 
291  void GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
292  DenseMatrix &subm);
293 
294  inline void SetColPtr(const int row) const;
295  inline void ClearColPtr() const;
296  inline double &SearchRow(const int col);
297  inline void _Add_(const int col, const double a)
298  { SearchRow(col) += a; }
299  inline void _Set_(const int col, const double a)
300  { SearchRow(col) = a; }
301  inline double _Get_(const int col) const;
302 
303  inline double &SearchRow(const int row, const int col);
304  inline void _Add_(const int row, const int col, const double a)
305  { SearchRow(row, col) += a; }
306  inline void _Set_(const int row, const int col, const double a)
307  { SearchRow(row, col) = a; }
308 
309  void Set(const int i, const int j, const double a);
310  void Add(const int i, const int j, const double a);
311 
312  void SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
313  const DenseMatrix &subm, int skip_zeros = 1);
314 
315  void SetSubMatrixTranspose(const Array<int> &rows, const Array<int> &cols,
316  const DenseMatrix &subm, int skip_zeros = 1);
317 
318  void AddSubMatrix(const Array<int> &rows, const Array<int> &cols,
319  const DenseMatrix &subm, int skip_zeros = 1);
320 
321  bool RowIsEmpty(const int row) const;
322 
323  /// Extract all column indices and values from a given row.
324  /** If the matrix is finalized (i.e. in CSR format), @a cols and @a srow will
325  simply be references to the specific portion of the #J and #A arrays.
326  As required by the AbstractSparseMatrix interface this method returns:
327  - 0, if @a cols and @a srow are copies of the values in the matrix, i.e.
328  when the matrix is open.
329  - 1, if @a cols and @a srow are views of the values in the matrix, i.e.
330  when the matrix is finalized. */
331  virtual int GetRow(const int row, Array<int> &cols, Vector &srow) const;
332 
333  void SetRow(const int row, const Array<int> &cols, const Vector &srow);
334  void AddRow(const int row, const Array<int> &cols, const Vector &srow);
335 
336  void ScaleRow(const int row, const double scale);
337  // this = diag(sl) * this;
338  void ScaleRows(const Vector & sl);
339  // this = this * diag(sr);
340  void ScaleColumns(const Vector & sr);
341 
342  /** Add the sparse matrix 'B' to '*this'. This operation will cause an error
343  if '*this' is finalized and 'B' has larger sparsity pattern. */
345 
346  /** Add the sparse matrix 'B' scaled by the scalar 'a' into '*this'.
347  Only entries in the sparsity pattern of '*this' are added. */
348  void Add(const double a, const SparseMatrix &B);
349 
350  SparseMatrix &operator=(double a);
351 
352  SparseMatrix &operator*=(double a);
353 
354  /// Prints matrix to stream out.
355  void Print(std::ostream &out = mfem::out, int width_ = 4) const;
356 
357  /// Prints matrix in matlab format.
358  void PrintMatlab(std::ostream &out = mfem::out) const;
359 
360  /// Prints matrix in Matrix Market sparse format.
361  void PrintMM(std::ostream &out = mfem::out) const;
362 
363  /// Prints matrix to stream out in hypre_CSRMatrix format.
364  void PrintCSR(std::ostream &out) const;
365 
366  /// Prints a sparse matrix to stream out in CSR format.
367  void PrintCSR2(std::ostream &out) const;
368 
369  /// Print various sparse matrix staticstics.
370  void PrintInfo(std::ostream &out) const;
371 
372  /// Walks the sparse matrix
373  int Walk(int &i, int &j, double &a);
374 
375  /// Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix
376  double IsSymmetric() const;
377 
378  /// (*this) = 1/2 ((*this) + (*this)^t)
379  void Symmetrize();
380 
381  /// Returns the number of the nonzero elements in the matrix
382  virtual int NumNonZeroElems() const;
383 
384  double MaxNorm() const;
385 
386  /// Count the number of entries with |a_ij| <= tol.
387  int CountSmallElems(double tol) const;
388 
389  /// Count the number of entries that are NOT finite, i.e. Inf or Nan.
390  int CheckFinite() const;
391 
392  /// Set the graph ownership flag (I and J arrays).
393  void SetGraphOwner(bool ownij) { ownGraph = ownij; }
394  /// Set the data ownership flag (A array).
395  void SetDataOwner(bool owna) { ownData = owna; }
396  /// Get the graph ownership flag (I and J arrays).
397  bool OwnsGraph() const { return ownGraph; }
398  /// Get the data ownership flag (A array).
399  bool OwnsData() const { return ownData; }
400  /// Lose the ownership of the graph (I, J) and data (A) arrays.
401  void LoseData() { ownGraph = ownData = false; }
402 
403  void Swap(SparseMatrix &other);
404 
405  /// Destroys sparse matrix.
406  virtual ~SparseMatrix() { Destroy(); }
407 
408  Type GetType() const { return MFEM_SPARSEMAT; }
409 };
410 
411 /// Applies f() to each element of the matrix (after it is finalized).
412 void SparseMatrixFunction(SparseMatrix &S, double (*f)(double));
413 
414 
415 /// Transpose of a sparse matrix. A must be finalized.
416 SparseMatrix *Transpose(const SparseMatrix &A);
417 /// Transpose of a sparse matrix. A does not need to be a CSR matrix.
418 SparseMatrix *TransposeAbstractSparseMatrix (const AbstractSparseMatrix &A,
419  int useActualWidth);
420 
421 /** Matrix product A.B.
422  If OAB is not NULL, we assume it has the structure
423  of A.B and store the result in OAB.
424  If OAB is NULL, we create a new SparseMatrix to store
425  the result and return a pointer to it.
426  All matrices must be finalized. */
427 SparseMatrix *Mult(const SparseMatrix &A, const SparseMatrix &B,
428  SparseMatrix *OAB = NULL);
429 
430 /// Matrix product of sparse matrices. A and B do not need to be CSR matrices
431 SparseMatrix *MultAbstractSparseMatrix (const AbstractSparseMatrix &A,
432  const AbstractSparseMatrix &B);
433 
434 /// Matrix product A.B
435 DenseMatrix *Mult(const SparseMatrix &A, DenseMatrix &B);
436 
437 /// RAP matrix product (with R=P^T)
438 DenseMatrix *RAP(const SparseMatrix &A, DenseMatrix &P);
439 
440 /** RAP matrix product (with P=R^T). ORAP is like OAB above.
441  All matrices must be finalized. */
442 SparseMatrix *RAP(const SparseMatrix &A, const SparseMatrix &R,
443  SparseMatrix *ORAP = NULL);
444 
445 /// General RAP with given R^T, A and P
446 SparseMatrix *RAP(const SparseMatrix &Rt, const SparseMatrix &A,
447  const SparseMatrix &P);
448 
449 /// Matrix multiplication A^t D A. All matrices must be finalized.
450 SparseMatrix *Mult_AtDA(const SparseMatrix &A, const Vector &D,
451  SparseMatrix *OAtDA = NULL);
452 
453 
454 /// Matrix addition result = A + B.
455 SparseMatrix * Add(const SparseMatrix & A, const SparseMatrix & B);
456 /// Matrix addition result = a*A + b*B
457 SparseMatrix * Add(double a, const SparseMatrix & A, double b,
458  const SparseMatrix & B);
459 /// Matrix addition result = sum_i A_i
460 SparseMatrix * Add(Array<SparseMatrix *> & Ai);
461 
462 
463 // Inline methods
464 
465 inline void SparseMatrix::SetColPtr(const int row) const
466 {
467  if (Rows)
468  {
469  if (ColPtrNode == NULL)
470  {
471  ColPtrNode = new RowNode *[width];
472  for (int i = 0; i < width; i++)
473  {
474  ColPtrNode[i] = NULL;
475  }
476  }
477  for (RowNode *node_p = Rows[row]; node_p != NULL; node_p = node_p->Prev)
478  {
479  ColPtrNode[node_p->Column] = node_p;
480  }
481  }
482  else
483  {
484  if (ColPtrJ == NULL)
485  {
486  ColPtrJ = new int[width];
487  for (int i = 0; i < width; i++)
488  {
489  ColPtrJ[i] = -1;
490  }
491  }
492  for (int j = I[row], end = I[row+1]; j < end; j++)
493  {
494  ColPtrJ[J[j]] = j;
495  }
496  }
497  current_row = row;
498 }
499 
500 inline void SparseMatrix::ClearColPtr() const
501 {
502  if (Rows)
503  for (RowNode *node_p = Rows[current_row]; node_p != NULL;
504  node_p = node_p->Prev)
505  {
506  ColPtrNode[node_p->Column] = NULL;
507  }
508  else
509  for (int j = I[current_row], end = I[current_row+1]; j < end; j++)
510  {
511  ColPtrJ[J[j]] = -1;
512  }
513 }
514 
515 inline double &SparseMatrix::SearchRow(const int col)
516 {
517  if (Rows)
518  {
519  RowNode *node_p = ColPtrNode[col];
520  if (node_p == NULL)
521  {
522 #ifdef MFEM_USE_MEMALLOC
523  node_p = NodesMem->Alloc();
524 #else
525  node_p = new RowNode;
526 #endif
527  node_p->Prev = Rows[current_row];
528  node_p->Column = col;
529  node_p->Value = 0.0;
530  Rows[current_row] = ColPtrNode[col] = node_p;
531  }
532  return node_p->Value;
533  }
534  else
535  {
536  const int j = ColPtrJ[col];
537  MFEM_VERIFY(j != -1, "Entry for column " << col << " is not allocated.");
538  return A[j];
539  }
540 }
541 
542 inline double SparseMatrix::_Get_(const int col) const
543 {
544  if (Rows)
545  {
546  RowNode *node_p = ColPtrNode[col];
547  return (node_p == NULL) ? 0.0 : node_p->Value;
548  }
549  else
550  {
551  const int j = ColPtrJ[col];
552  return (j == -1) ? 0.0 : A[j];
553  }
554 }
555 
556 inline double &SparseMatrix::SearchRow(const int row, const int col)
557 {
558  if (Rows)
559  {
560  RowNode *node_p;
561 
562  for (node_p = Rows[row]; 1; node_p = node_p->Prev)
563  {
564  if (node_p == NULL)
565  {
566 #ifdef MFEM_USE_MEMALLOC
567  node_p = NodesMem->Alloc();
568 #else
569  node_p = new RowNode;
570 #endif
571  node_p->Prev = Rows[row];
572  node_p->Column = col;
573  node_p->Value = 0.0;
574  Rows[row] = node_p;
575  break;
576  }
577  else if (node_p->Column == col)
578  {
579  break;
580  }
581  }
582  return node_p->Value;
583  }
584  else
585  {
586  int *Ip = I+row, *Jp = J;
587  for (int k = Ip[0], end = Ip[1]; k < end; k++)
588  {
589  if (Jp[k] == col)
590  {
591  return A[k];
592  }
593  }
594  MFEM_ABORT("Could not find entry for row = " << row << ", col = " << col);
595  }
596  return A[0];
597 }
598 
599 /// Specialization of the template function Swap<> for class SparseMatrix
600 template<> inline void Swap<SparseMatrix>(SparseMatrix &a, SparseMatrix &b)
601 {
602  a.Swap(b);
603 }
604 
605 }
606 
607 #endif
RowNode ** ColPtrNode
Definition: sparsemat.hpp:65
Elem * Alloc()
Definition: mem_alloc.hpp:146
void _Add_(const int row, const int col, const double a)
Definition: sparsemat.hpp:304
double GetJacobiScaling() const
Determine appropriate scaling for Jacobi iteration.
Definition: sparsemat.cpp:1706
int CheckFinite() const
Count the number of entries that are NOT finite, i.e. Inf or Nan.
Definition: sparsemat.cpp:1045
SparseMatrix * TransposeAbstractSparseMatrix(const AbstractSparseMatrix &A, int useActualWidth)
Transpose of a sparse matrix. A does not need to be a CSR matrix.
Definition: sparsemat.cpp:2659
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:213
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:972
void Add(const int i, const int j, const double a)
Definition: sparsemat.cpp:1907
void _Add_(const int col, const double a)
Definition: sparsemat.hpp:297
void Clear()
Clear the contents of the SparseMatrix.
Definition: sparsemat.hpp:125
Type GetType() const
Definition: sparsemat.hpp:408
void _Set_(const int row, const int col, const double a)
Definition: sparsemat.hpp:306
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1476
void DiagScale(const Vector &b, Vector &x, double sc=1.0) const
Definition: sparsemat.cpp:1770
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:278
void MakeRef(const SparseMatrix &master)
Clear the contents of the SparseMatrix and make it a reference to master.
Definition: sparsemat.cpp:187
double & SearchRow(const int col)
Definition: sparsemat.hpp:515
bool Empty() const
Check if the SparseMatrix is empty.
Definition: sparsemat.hpp:128
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:260
int * I
Array with size (height+1) containing the row offsets.
Definition: sparsemat.hpp:50
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, but treat all elements as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:633
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:468
Abstract data type for sparse matrices.
Definition: matrix.hpp:69
SparseMatrix & operator+=(const SparseMatrix &B)
Definition: sparsemat.cpp:2249
void EliminateRowCol(int rc, const double sol, Vector &rhs, int d=0)
Definition: sparsemat.cpp:1162
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
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:481
double MaxNorm() const
Definition: sparsemat.cpp:994
virtual ~SparseMatrix()
Destroys sparse matrix.
Definition: sparsemat.hpp:406
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:558
Abstract data type for matrix inverse.
Definition: matrix.hpp:58
void EliminateRowColMultipleRHS(int rc, const Vector &sol, DenseMatrix &rhs, int d=0)
Definition: sparsemat.cpp:1242
void PrintInfo(std::ostream &out) const
Print various sparse matrix staticstics.
Definition: sparsemat.cpp:2469
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: sparsemat.cpp:379
void LoseData()
Lose the ownership of the graph (I, J) and data (A) arrays.
Definition: sparsemat.hpp:401
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2130
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:274
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:326
void PrintMM(std::ostream &out=mfem::out) const
Prints matrix in Matrix Market sparse format.
Definition: sparsemat.cpp:2400
void _Set_(const int col, const double a)
Definition: sparsemat.hpp:299
void BooleanMultTranspose(const Array< int > &x, Array< int > &y) const
y = At * x, but treat all elements as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:656
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:3004
void PartAddMult(const Array< int > &rows, const Vector &x, Vector &y, const double a=1.0) const
Definition: sparsemat.cpp:615
bool RowIsEmpty(const int row) const
Definition: sparsemat.cpp:2022
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Definition: sparsemat.cpp:825
void ScaleRow(const int row, const double scale)
Definition: sparsemat.cpp:2162
RowNodeAlloc * NodesMem
Definition: sparsemat.hpp:69
int ActualWidth()
Returns the actual Width of the matrix.
Definition: sparsemat.cpp:2568
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
Definition: sparsemat.cpp:957
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:355
void SetGraphOwner(bool ownij)
Set the graph ownership flag (I and J arrays).
Definition: sparsemat.hpp:393
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2942
bool OwnsData() const
Get the data ownership flag (A array).
Definition: sparsemat.hpp:399
void ClearColPtr() const
Definition: sparsemat.hpp:500
void EliminateCol(int col)
Definition: sparsemat.cpp:1117
void ScaleColumns(const Vector &sr)
Definition: sparsemat.cpp:2221
void PrintCSR2(std::ostream &out) const
Prints a sparse matrix to stream out in CSR format.
Definition: sparsemat.cpp:2444
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:2043
bool isSorted
Are the columns sorted already.
Definition: sparsemat.hpp:78
SparseMatrix()
Create an empty SparseMatrix.
Definition: sparsemat.hpp:85
Data type sparse matrix.
Definition: sparsemat.hpp:38
bool ownData
Say whether we own the pointers for A (should we free them?).
Definition: sparsemat.hpp:75
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
void EliminateCols(Array< int > &cols, Vector *x=NULL, Vector *b=NULL)
Eliminate all columns 'i' for which cols[i] != 0.
Definition: sparsemat.cpp:1131
bool OwnsGraph() const
Get the graph ownership flag (I and J arrays).
Definition: sparsemat.hpp:397
void SetWidth(int width_=-1)
Change the width of a SparseMatrix.
Definition: sparsemat.cpp:288
bool areColumnsSorted() const
Definition: sparsemat.hpp:284
double IsSymmetric() const
Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix.
Definition: sparsemat.cpp:915
void SetDiagIdentity()
If a row contains only one diag entry of zero, set it to 1.
Definition: sparsemat.cpp:1520
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:408
void ScaleRows(const Vector &sl)
Definition: sparsemat.cpp:2190
void Swap< SparseMatrix >(SparseMatrix &a, SparseMatrix &b)
Specialization of the template function Swap<> for class SparseMatrix.
Definition: sparsemat.hpp:600
void PrintMatlab(std::ostream &out=mfem::out) const
Prints matrix in matlab format.
Definition: sparsemat.cpp:2382
void EliminateZeroRows()
If a row contains only zeros, set its diagonal to 1.
Definition: sparsemat.cpp:1529
double * GetData() const
Return element data, i.e. array A.
Definition: sparsemat.hpp:135
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1850
int * GetI() const
Return the array I.
Definition: sparsemat.hpp:131
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const
Definition: sparsemat.cpp:1739
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1926
bool Finalized() const
Definition: sparsemat.hpp:283
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:123
void SparseMatrixFunction(SparseMatrix &S, double(*f)(double))
Applies f() to each element of the matrix (after it is finalized).
Definition: sparsemat.cpp:2594
virtual MatrixInverse * Inverse() const
Returns a pointer to approximation of the matrix inverse.
Definition: sparsemat.cpp:1069
void EliminateRowColDiag(int rc, double value)
Perform elimination and set the diagonal entry to the given value.
Definition: sparsemat.cpp:1392
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:475
void SetColPtr(const int row) const
Definition: sparsemat.hpp:465
SparseMatrix & operator=(double a)
Definition: sparsemat.cpp:2302
void SetSubMatrixTranspose(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1959
void Set(const int i, const int j, const double a)
Definition: sparsemat.cpp:1888
int * J
Array with size I[height], containing the column indices for all matrix entries, as indexed by the I ...
Definition: sparsemat.hpp:53
int Size() const
For backward compatibility, define Size() to be synonym of Height().
Definition: sparsemat.hpp:122
void Gauss_Seidel_back(const Vector &x, Vector &y) const
Definition: sparsemat.cpp:1631
void Gauss_Seidel_forw(const Vector &x, Vector &y) const
Gauss-Seidel forward and backward iterations over a vector x.
Definition: sparsemat.cpp:1556
RowNode ** Rows
Array of linked lists, one for every row. This array represents the linked list (LIL) storage format...
Definition: sparsemat.hpp:61
void EliminateRow(int row, const double sol, Vector &rhs)
Eliminates a column from the transpose matrix.
Definition: sparsemat.cpp:1074
double * A
Array with size I[height], containing the actual entries of the sparse matrix, as indexed by the I ar...
Definition: sparsemat.hpp:56
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:447
MemAlloc< RowNode, 1024 > RowNodeAlloc
Definition: sparsemat.hpp:68
double GetRowNorml1(int irow) const
For i = irow compute .
Definition: sparsemat.cpp:724
double & operator()(int i, int j)
Returns reference to A[i][j].
Definition: sparsemat.cpp:389
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:552
bool ownGraph
Say whether we own the pointers for I and J (should we free them?).
Definition: sparsemat.hpp:73
SparseMatrix & operator*=(double a)
Definition: sparsemat.cpp:2320
void SetDataOwner(bool owna)
Set the data ownership flag (A array).
Definition: sparsemat.hpp:395
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2082
Vector data type.
Definition: vector.hpp:41
int Walk(int &i, int &j, double &a)
Walks the sparse matrix.
void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: sparsemat.cpp:2338
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm)
Definition: sparsemat.cpp:1994
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:1802
int MaxRowSize() const
Returns the maximum number of elements among all rows.
Definition: sparsemat.cpp:236
void Jacobi3(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:1826
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:2862
double _Get_(const int col) const
Definition: sparsemat.hpp:542
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:64
void Swap(SparseMatrix &other)
Definition: sparsemat.cpp:3132
int * GetJ() const
Return the array J.
Definition: sparsemat.hpp:133
void PrintCSR(std::ostream &out) const
Prints matrix to stream out in hypre_CSRMatrix format.
Definition: sparsemat.cpp:2420
double InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Definition: sparsemat.cpp:679
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25
int CountSmallElems(double tol) const
Count the number of entries with |a_ij| <= tol.
Definition: sparsemat.cpp:1017
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Definition: sparsemat.cpp:597
void GetRowSums(Vector &x) const
For all i compute .
Definition: sparsemat.cpp:705