MFEM  v3.4
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
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  /// Assignment operator: deep copy
116  SparseMatrix& operator=(const SparseMatrix &rhs);
117 
118  /** @brief Clear the contents of the SparseMatrix and make it a reference to
119  @a master */
120  /** After this call, the matrix will point to the same data as @a master but
121  it will not own its data. The @a master must be finalized. */
122  void MakeRef(const SparseMatrix &master);
123 
124  /// For backward compatibility, define Size() to be synonym of Height().
125  int Size() const { return Height(); }
126 
127  /// Clear the contents of the SparseMatrix.
128  void Clear() { Destroy(); SetEmpty(); }
129 
130  /// Check if the SparseMatrix is empty.
131  bool Empty() const { return (A == NULL) && (Rows == NULL); }
132 
133  /// Return the array #I
134  inline int *GetI() const { return I; }
135  /// Return the array #J
136  inline int *GetJ() const { return J; }
137  /// Return element data, i.e. array #A
138  inline double *GetData() const { return A; }
139  /// Returns the number of elements in row @a i
140  int RowSize(const int i) const;
141  /// Returns the maximum number of elements among all rows
142  int MaxRowSize() const;
143  /// Return a pointer to the column indices in a row
144  int *GetRowColumns(const int row);
145  const int *GetRowColumns(const int row) const;
146  /// Return a pointer to the entries in a row
147  double *GetRowEntries(const int row);
148  const double *GetRowEntries(const int row) const;
149 
150  /// Change the width of a SparseMatrix.
151  /*!
152  * If width_ = -1 (DEFAULT), this routine will set the new width
153  * to the actual Width of the matrix awidth = max(J) + 1.
154  * Values 0 <= width_ < awidth are not allowed (error check in Debug Mode only)
155  *
156  * This method can be called for matrices finalized or not.
157  */
158  void SetWidth(int width_ = -1);
159 
160  /// Returns the actual Width of the matrix.
161  /*! This method can be called for matrices finalized or not. */
162  int ActualWidth();
163 
164  /// Sort the column indices corresponding to each row.
165  void SortColumnIndices();
166 
167  /** @brief Move the diagonal entry to the first position in each row,
168  preserving the order of the rest of the columns. */
169  void MoveDiagonalFirst();
170 
171  /// Returns reference to a_{ij}.
172  virtual double &Elem(int i, int j);
173 
174  /// Returns constant reference to a_{ij}.
175  virtual const double &Elem(int i, int j) const;
176 
177  /// Returns reference to A[i][j].
178  double &operator()(int i, int j);
179 
180  /// Returns reference to A[i][j].
181  const double &operator()(int i, int j) const;
182 
183  /// Returns the Diagonal of A
184  void GetDiag(Vector & d) const;
185 
186  /// Matrix vector multiplication.
187  virtual void Mult(const Vector &x, Vector &y) const;
188 
189  /// y += A * x (default) or y += a * A * x
190  void AddMult(const Vector &x, Vector &y, const double a = 1.0) const;
191 
192  /// Multiply a vector with the transposed matrix. y = At * x
193  void MultTranspose(const Vector &x, Vector &y) const;
194 
195  /// y += At * x (default) or y += a * At * x
196  void AddMultTranspose(const Vector &x, Vector &y,
197  const double a = 1.0) const;
198 
199  void PartMult(const Array<int> &rows, const Vector &x, Vector &y) const;
200  void PartAddMult(const Array<int> &rows, const Vector &x, Vector &y,
201  const double a=1.0) const;
202 
203  /// y = A * x, but treat all elements as booleans (zero=false, nonzero=true).
204  void BooleanMult(const Array<int> &x, Array<int> &y) const;
205  /// y = At * x, but treat all elements as booleans (zero=false, nonzero=true).
206  void BooleanMultTranspose(const Array<int> &x, Array<int> &y) const;
207 
208  /// Compute y^t A x
209  double InnerProduct(const Vector &x, const Vector &y) const;
210 
211  /// For all i compute \f$ x_i = \sum_j A_{ij} \f$
212  void GetRowSums(Vector &x) const;
213  /// For i = irow compute \f$ x_i = \sum_j | A_{i, j} | \f$
214  double GetRowNorml1(int irow) const;
215 
216  /// This virtual method is not supported: it always returns NULL.
217  virtual MatrixInverse *Inverse() const;
218 
219  /// Eliminates a column from the transpose matrix.
220  void EliminateRow(int row, const double sol, Vector &rhs);
221 
222  /// Eliminates a row from the matrix.
223  /*!
224  * - If @a dpolicy = #DIAG_ZERO, all the entries in the row will be set to 0.
225  * - If @a dpolicy = #DIAG_ONE (matrix must be square), the diagonal entry
226  * will be set equal to 1 and all other entries in the row to 0.
227  * - The policy #DIAG_KEEP is not supported.
228  */
229  void EliminateRow(int row, DiagonalPolicy dpolicy = DIAG_ZERO);
230 
231  /// Eliminates the column @a col from the matrix.
232  /** - If @a dpolicy = #DIAG_ZERO, all entries in the column will be set to 0.
233  - If @a dpolicy = #DIAG_ONE (matrix must be square), the diagonal entry
234  will be set equal to 1 and all other entries in the column to 0.
235  - The policy #DIAG_KEEP is not supported. */
236  void EliminateCol(int col, DiagonalPolicy dpolicy = DIAG_ZERO);
237 
238  /// Eliminate all columns i for which @a cols[i] != 0.
239  /** Elimination of a column means that all entries in the column are set to
240  zero. In addition, if the pointers @a x and @a b are not NULL, the
241  eliminated matrix entries are multiplied by the corresponding solution
242  value in @a *x and subtracted from the r.h.s. vector, @a *b. */
243  void EliminateCols(const Array<int> &cols, const Vector *x = NULL,
244  Vector *b = NULL);
245 
246  /// Eliminate row @a rc and column @a rc and modify the @a rhs using @a sol.
247  /** Eliminates the column @a rc to the @a rhs, deletes the row @a rc and
248  replaces the element (rc,rc) with 1.0; assumes that element (i,rc)
249  is assembled if and only if the element (rc,i) is assembled.
250  By default, elements (rc,rc) are set to 1.0, although this behavior
251  can be adjusted by changing the @a dpolicy parameter. */
252  void EliminateRowCol(int rc, const double sol, Vector &rhs,
253  DiagonalPolicy dpolicy = DIAG_ONE);
254 
255  /** @brief Similar to
256  EliminateRowCol(int, const double, Vector &, DiagonalPolicy), but
257  multiple values for eliminated unknowns are accepted, and accordingly
258  multiple right-hand-sides are used. */
259  void EliminateRowColMultipleRHS(int rc, const Vector &sol,
260  DenseMatrix &rhs,
261  DiagonalPolicy dpolicy = DIAG_ONE);
262 
263  /// Perform elimination and set the diagonal entry to the given value
264  void EliminateRowColDiag(int rc, double value);
265 
266  /// Eliminate row @a rc and column @a rc.
267  void EliminateRowCol(int rc, DiagonalPolicy dpolicy = DIAG_ONE);
268 
269  /** @brief Similar to EliminateRowCol(int, DiagonalPolicy) + save the
270  eliminated entries into @a Ae so that (*this) + Ae is equal to the
271  original matrix */
272  void EliminateRowCol(int rc, SparseMatrix &Ae,
273  DiagonalPolicy dpolicy = DIAG_ONE);
274 
275  /// If a row contains only one diag entry of zero, set it to 1.
276  void SetDiagIdentity();
277  /// If a row contains only zeros, set its diagonal to 1.
278  virtual void EliminateZeroRows(const double threshold = 1e-12);
279 
280  /// Gauss-Seidel forward and backward iterations over a vector x.
281  void Gauss_Seidel_forw(const Vector &x, Vector &y) const;
282  void Gauss_Seidel_back(const Vector &x, Vector &y) const;
283 
284  /// Determine appropriate scaling for Jacobi iteration
285  double GetJacobiScaling() const;
286  /** One scaled Jacobi iteration for the system A x = b.
287  x1 = x0 + sc D^{-1} (b - A x0) where D is the diag of A. */
288  void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const;
289 
290  void DiagScale(const Vector &b, Vector &x, double sc = 1.0) const;
291 
292  /** x1 = x0 + sc D^{-1} (b - A x0) where \f$ D_{ii} = \sum_j |A_{ij}| \f$. */
293  void Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
294  double sc = 1.0) const;
295 
296  /** x1 = x0 + sc D^{-1} (b - A x0) where \f$ D_{ii} = \sum_j A_{ij} \f$. */
297  void Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
298  double sc = 1.0) const;
299 
300  /** @brief Finalize the matrix initialization, switching the storage format
301  from LIL to CSR. */
302  /** This method should be called once, after the matrix has been initialized.
303  Internally, this method converts the matrix from row-wise linked list
304  (LIL) format into CSR (compressed sparse row) format. */
305  virtual void Finalize(int skip_zeros = 1) { Finalize(skip_zeros, false); }
306 
307  /// A slightly more general version of the Finalize(int) method.
308  void Finalize(int skip_zeros, bool fix_empty_rows);
309 
310  bool Finalized() const { return (A != NULL); }
311  bool areColumnsSorted() const { return isSorted; }
312 
313  /** Split the matrix into M x N blocks of sparse matrices in CSR format.
314  The 'blocks' array is M x N (i.e. M and N are determined by its
315  dimensions) and its entries are overwritten by the new blocks. */
316  void GetBlocks(Array2D<SparseMatrix *> &blocks) const;
317 
318  void GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
319  DenseMatrix &subm) const;
320 
321  inline void SetColPtr(const int row) const;
322  inline void ClearColPtr() const;
323  inline double &SearchRow(const int col);
324  inline void _Add_(const int col, const double a)
325  { SearchRow(col) += a; }
326  inline void _Set_(const int col, const double a)
327  { SearchRow(col) = a; }
328  inline double _Get_(const int col) const;
329 
330  inline double &SearchRow(const int row, const int col);
331  inline void _Add_(const int row, const int col, const double a)
332  { SearchRow(row, col) += a; }
333  inline void _Set_(const int row, const int col, const double a)
334  { SearchRow(row, col) = a; }
335 
336  void Set(const int i, const int j, const double a);
337  void Add(const int i, const int j, const double a);
338 
339  void SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
340  const DenseMatrix &subm, int skip_zeros = 1);
341 
342  void SetSubMatrixTranspose(const Array<int> &rows, const Array<int> &cols,
343  const DenseMatrix &subm, int skip_zeros = 1);
344 
345  void AddSubMatrix(const Array<int> &rows, const Array<int> &cols,
346  const DenseMatrix &subm, int skip_zeros = 1);
347 
348  bool RowIsEmpty(const int row) const;
349 
350  /// Extract all column indices and values from a given row.
351  /** If the matrix is finalized (i.e. in CSR format), @a cols and @a srow will
352  simply be references to the specific portion of the #J and #A arrays.
353  As required by the AbstractSparseMatrix interface this method returns:
354  - 0, if @a cols and @a srow are copies of the values in the matrix, i.e.
355  when the matrix is open.
356  - 1, if @a cols and @a srow are views of the values in the matrix, i.e.
357  when the matrix is finalized. */
358  virtual int GetRow(const int row, Array<int> &cols, Vector &srow) const;
359 
360  void SetRow(const int row, const Array<int> &cols, const Vector &srow);
361  void AddRow(const int row, const Array<int> &cols, const Vector &srow);
362 
363  void ScaleRow(const int row, const double scale);
364  // this = diag(sl) * this;
365  void ScaleRows(const Vector & sl);
366  // this = this * diag(sr);
367  void ScaleColumns(const Vector & sr);
368 
369  /** Add the sparse matrix 'B' to '*this'. This operation will cause an error
370  if '*this' is finalized and 'B' has larger sparsity pattern. */
372 
373  /** Add the sparse matrix 'B' scaled by the scalar 'a' into '*this'.
374  Only entries in the sparsity pattern of '*this' are added. */
375  void Add(const double a, const SparseMatrix &B);
376 
377  SparseMatrix &operator=(double a);
378 
379  SparseMatrix &operator*=(double a);
380 
381  /// Prints matrix to stream out.
382  void Print(std::ostream &out = mfem::out, int width_ = 4) const;
383 
384  /// Prints matrix in matlab format.
385  void PrintMatlab(std::ostream &out = mfem::out) const;
386 
387  /// Prints matrix in Matrix Market sparse format.
388  void PrintMM(std::ostream &out = mfem::out) const;
389 
390  /// Prints matrix to stream out in hypre_CSRMatrix format.
391  void PrintCSR(std::ostream &out) const;
392 
393  /// Prints a sparse matrix to stream out in CSR format.
394  void PrintCSR2(std::ostream &out) const;
395 
396  /// Print various sparse matrix staticstics.
397  void PrintInfo(std::ostream &out) const;
398 
399  /// Walks the sparse matrix
400  int Walk(int &i, int &j, double &a);
401 
402  /// Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix
403  double IsSymmetric() const;
404 
405  /// (*this) = 1/2 ((*this) + (*this)^t)
406  void Symmetrize();
407 
408  /// Returns the number of the nonzero elements in the matrix
409  virtual int NumNonZeroElems() const;
410 
411  double MaxNorm() const;
412 
413  /// Count the number of entries with |a_ij| <= tol.
414  int CountSmallElems(double tol) const;
415 
416  /// Count the number of entries that are NOT finite, i.e. Inf or Nan.
417  int CheckFinite() const;
418 
419  /// Set the graph ownership flag (I and J arrays).
420  void SetGraphOwner(bool ownij) { ownGraph = ownij; }
421  /// Set the data ownership flag (A array).
422  void SetDataOwner(bool owna) { ownData = owna; }
423  /// Get the graph ownership flag (I and J arrays).
424  bool OwnsGraph() const { return ownGraph; }
425  /// Get the data ownership flag (A array).
426  bool OwnsData() const { return ownData; }
427  /// Lose the ownership of the graph (I, J) and data (A) arrays.
428  void LoseData() { ownGraph = ownData = false; }
429 
430  void Swap(SparseMatrix &other);
431 
432  /// Destroys sparse matrix.
433  virtual ~SparseMatrix() { Destroy(); }
434 
435  Type GetType() const { return MFEM_SPARSEMAT; }
436 };
437 
438 /// Applies f() to each element of the matrix (after it is finalized).
439 void SparseMatrixFunction(SparseMatrix &S, double (*f)(double));
440 
441 
442 /// Transpose of a sparse matrix. A must be finalized.
443 SparseMatrix *Transpose(const SparseMatrix &A);
444 /// Transpose of a sparse matrix. A does not need to be a CSR matrix.
445 SparseMatrix *TransposeAbstractSparseMatrix (const AbstractSparseMatrix &A,
446  int useActualWidth);
447 
448 /** Matrix product A.B.
449  If OAB is not NULL, we assume it has the structure
450  of A.B and store the result in OAB.
451  If OAB is NULL, we create a new SparseMatrix to store
452  the result and return a pointer to it.
453  All matrices must be finalized. */
454 SparseMatrix *Mult(const SparseMatrix &A, const SparseMatrix &B,
455  SparseMatrix *OAB = NULL);
456 
457 /// Matrix product of sparse matrices. A and B do not need to be CSR matrices
458 SparseMatrix *MultAbstractSparseMatrix (const AbstractSparseMatrix &A,
459  const AbstractSparseMatrix &B);
460 
461 /// Matrix product A.B
462 DenseMatrix *Mult(const SparseMatrix &A, DenseMatrix &B);
463 
464 /// RAP matrix product (with R=P^T)
465 DenseMatrix *RAP(const SparseMatrix &A, DenseMatrix &P);
466 
467 /** RAP matrix product (with P=R^T). ORAP is like OAB above.
468  All matrices must be finalized. */
469 SparseMatrix *RAP(const SparseMatrix &A, const SparseMatrix &R,
470  SparseMatrix *ORAP = NULL);
471 
472 /// General RAP with given R^T, A and P
473 SparseMatrix *RAP(const SparseMatrix &Rt, const SparseMatrix &A,
474  const SparseMatrix &P);
475 
476 /// Matrix multiplication A^t D A. All matrices must be finalized.
477 SparseMatrix *Mult_AtDA(const SparseMatrix &A, const Vector &D,
478  SparseMatrix *OAtDA = NULL);
479 
480 
481 /// Matrix addition result = A + B.
482 SparseMatrix * Add(const SparseMatrix & A, const SparseMatrix & B);
483 /// Matrix addition result = a*A + b*B
484 SparseMatrix * Add(double a, const SparseMatrix & A, double b,
485  const SparseMatrix & B);
486 /// Matrix addition result = sum_i A_i
487 SparseMatrix * Add(Array<SparseMatrix *> & Ai);
488 
489 
490 // Inline methods
491 
492 inline void SparseMatrix::SetColPtr(const int row) const
493 {
494  if (Rows)
495  {
496  if (ColPtrNode == NULL)
497  {
498  ColPtrNode = new RowNode *[width];
499  for (int i = 0; i < width; i++)
500  {
501  ColPtrNode[i] = NULL;
502  }
503  }
504  for (RowNode *node_p = Rows[row]; node_p != NULL; node_p = node_p->Prev)
505  {
506  ColPtrNode[node_p->Column] = node_p;
507  }
508  }
509  else
510  {
511  if (ColPtrJ == NULL)
512  {
513  ColPtrJ = new int[width];
514  for (int i = 0; i < width; i++)
515  {
516  ColPtrJ[i] = -1;
517  }
518  }
519  for (int j = I[row], end = I[row+1]; j < end; j++)
520  {
521  ColPtrJ[J[j]] = j;
522  }
523  }
524  current_row = row;
525 }
526 
527 inline void SparseMatrix::ClearColPtr() const
528 {
529  if (Rows)
530  for (RowNode *node_p = Rows[current_row]; node_p != NULL;
531  node_p = node_p->Prev)
532  {
533  ColPtrNode[node_p->Column] = NULL;
534  }
535  else
536  for (int j = I[current_row], end = I[current_row+1]; j < end; j++)
537  {
538  ColPtrJ[J[j]] = -1;
539  }
540 }
541 
542 inline double &SparseMatrix::SearchRow(const int col)
543 {
544  if (Rows)
545  {
546  RowNode *node_p = ColPtrNode[col];
547  if (node_p == NULL)
548  {
549 #ifdef MFEM_USE_MEMALLOC
550  node_p = NodesMem->Alloc();
551 #else
552  node_p = new RowNode;
553 #endif
554  node_p->Prev = Rows[current_row];
555  node_p->Column = col;
556  node_p->Value = 0.0;
557  Rows[current_row] = ColPtrNode[col] = node_p;
558  }
559  return node_p->Value;
560  }
561  else
562  {
563  const int j = ColPtrJ[col];
564  MFEM_VERIFY(j != -1, "Entry for column " << col << " is not allocated.");
565  return A[j];
566  }
567 }
568 
569 inline double SparseMatrix::_Get_(const int col) const
570 {
571  if (Rows)
572  {
573  RowNode *node_p = ColPtrNode[col];
574  return (node_p == NULL) ? 0.0 : node_p->Value;
575  }
576  else
577  {
578  const int j = ColPtrJ[col];
579  return (j == -1) ? 0.0 : A[j];
580  }
581 }
582 
583 inline double &SparseMatrix::SearchRow(const int row, const int col)
584 {
585  if (Rows)
586  {
587  RowNode *node_p;
588 
589  for (node_p = Rows[row]; 1; node_p = node_p->Prev)
590  {
591  if (node_p == NULL)
592  {
593 #ifdef MFEM_USE_MEMALLOC
594  node_p = NodesMem->Alloc();
595 #else
596  node_p = new RowNode;
597 #endif
598  node_p->Prev = Rows[row];
599  node_p->Column = col;
600  node_p->Value = 0.0;
601  Rows[row] = node_p;
602  break;
603  }
604  else if (node_p->Column == col)
605  {
606  break;
607  }
608  }
609  return node_p->Value;
610  }
611  else
612  {
613  int *Ip = I+row, *Jp = J;
614  for (int k = Ip[0], end = Ip[1]; k < end; k++)
615  {
616  if (Jp[k] == col)
617  {
618  return A[k];
619  }
620  }
621  MFEM_ABORT("Could not find entry for row = " << row << ", col = " << col);
622  }
623  return A[0];
624 }
625 
626 /// Specialization of the template function Swap<> for class SparseMatrix
627 template<> inline void Swap<SparseMatrix>(SparseMatrix &a, SparseMatrix &b)
628 {
629  a.Swap(b);
630 }
631 
632 }
633 
634 #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:331
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:486
SparseMatrix * TransposeAbstractSparseMatrix(const AbstractSparseMatrix &A, int useActualWidth)
Transpose of a sparse matrix. A does not need to be a CSR matrix.
Definition: sparsemat.cpp:2787
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:983
void Add(const int i, const int j, const double a)
Definition: sparsemat.cpp:2035
int * GetI() const
Return the array I.
Definition: sparsemat.hpp:134
void _Add_(const int col, const double a)
Definition: sparsemat.hpp:324
void Clear()
Clear the contents of the SparseMatrix.
Definition: sparsemat.hpp:128
void EliminateCol(int col, DiagonalPolicy dpolicy=DIAG_ZERO)
Eliminates the column col from the matrix.
Definition: sparsemat.cpp:1133
int Size() const
For backward compatibility, define Size() to be synonym of Height().
Definition: sparsemat.hpp:125
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:458
void PrintMM(std::ostream &out=mfem::out) const
Prints matrix in Matrix Market sparse format.
Definition: sparsemat.cpp:2528
void _Set_(const int row, const int col, const double a)
Definition: sparsemat.hpp:333
bool OwnsData() const
Get the data ownership flag (A array).
Definition: sparsemat.hpp:426
HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1558
double GetRowNorml1(int irow) const
For i = irow compute .
Definition: sparsemat.cpp:735
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:305
void MakeRef(const SparseMatrix &master)
Clear the contents of the SparseMatrix and make it a reference to master.
Definition: sparsemat.cpp:197
double & SearchRow(const int col)
Definition: sparsemat.hpp:542
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:271
void PrintCSR(std::ostream &out) const
Prints matrix to stream out in hypre_CSRMatrix format.
Definition: sparsemat.cpp:2548
int * I
Array with size (height+1) containing the row offsets.
Definition: sparsemat.hpp:50
double GetJacobiScaling() const
Determine appropriate scaling for Jacobi iteration.
Definition: sparsemat.cpp:1834
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:224
Set the diagonal value to zero.
Definition: matrix.hpp:34
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
void SetColPtr(const int row) const
Definition: sparsemat.hpp:492
Abstract data type for sparse matrices.
Definition: matrix.hpp:77
double IsSymmetric() const
Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix.
Definition: sparsemat.cpp:926
SparseMatrix & operator+=(const SparseMatrix &B)
Definition: sparsemat.cpp:2377
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual ~SparseMatrix()
Destroys sparse matrix.
Definition: sparsemat.hpp:433
Abstract data type for matrix inverse.
Definition: matrix.hpp:66
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:1300
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:492
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:2171
double MaxNorm() const
Definition: sparsemat.cpp:1005
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: sparsemat.cpp:390
void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: sparsemat.cpp:2466
void LoseData()
Lose the ownership of the graph (I, J) and data (A) arrays.
Definition: sparsemat.hpp:428
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
Definition: sparsemat.cpp:2122
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2258
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:285
void SortColumnIndices()
Sort the column indices corresponding to each row.
Definition: sparsemat.cpp:337
void _Set_(const int col, const double a)
Definition: sparsemat.hpp:326
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:3132
void Gauss_Seidel_back(const Vector &x, Vector &y) const
Definition: sparsemat.cpp:1759
void ScaleRow(const int row, const double scale)
Definition: sparsemat.cpp:2290
RowNodeAlloc * NodesMem
Definition: sparsemat.hpp:69
virtual MatrixInverse * Inverse() const
This virtual method is not supported: it always returns NULL.
Definition: sparsemat.cpp:1080
int ActualWidth()
Returns the actual Width of the matrix.
Definition: sparsemat.cpp:2696
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
Definition: sparsemat.cpp:968
void Jacobi3(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:1954
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:366
void SetGraphOwner(bool ownij)
Set the graph ownership flag (I and J arrays).
Definition: sparsemat.hpp:420
ID for class SparseMatrix.
Definition: operator.hpp:127
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2942
void GetBlocks(Array2D< SparseMatrix *> &blocks) const
Definition: sparsemat.cpp:836
void ScaleColumns(const Vector &sr)
Definition: sparsemat.cpp:2349
void GetRowSums(Vector &x) const
For all i compute .
Definition: sparsemat.cpp:716
bool Empty() const
Check if the SparseMatrix is empty.
Definition: sparsemat.hpp:131
bool isSorted
Are the columns sorted already.
Definition: sparsemat.hpp:78
void ClearColPtr() const
Definition: sparsemat.hpp:527
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
void SetWidth(int width_=-1)
Change the width of a SparseMatrix.
Definition: sparsemat.cpp:299
void PrintMatlab(std::ostream &out=mfem::out) const
Prints matrix in matlab format.
Definition: sparsemat.cpp:2510
int CountSmallElems(double tol) const
Count the number of entries with |a_ij| <= tol.
Definition: sparsemat.cpp:1028
double InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Definition: sparsemat.cpp:690
void DiagScale(const Vector &b, Vector &x, double sc=1.0) const
Definition: sparsemat.cpp:1898
bool OwnsGraph() const
Get the graph ownership flag (I and J arrays).
Definition: sparsemat.hpp:424
void SetDiagIdentity()
If a row contains only one diag entry of zero, set it to 1.
Definition: sparsemat.cpp:1648
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:418
void ScaleRows(const Vector &sl)
Definition: sparsemat.cpp:2318
void Swap< SparseMatrix >(SparseMatrix &a, SparseMatrix &b)
Specialization of the template function Swap<> for class SparseMatrix.
Definition: sparsemat.hpp:627
bool RowIsEmpty(const int row) const
Definition: sparsemat.cpp:2150
bool Finalized() const
Definition: sparsemat.hpp:310
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1, double sc=1.0) const
Definition: sparsemat.cpp:1930
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1978
Set the diagonal value to one.
Definition: matrix.hpp:35
Dynamic 2D array using row-major layout.
Definition: array.hpp:289
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2054
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:563
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:124
void SparseMatrixFunction(SparseMatrix &S, double(*f)(double))
Applies f() to each element of the matrix (after it is finalized).
Definition: sparsemat.cpp:2722
double _Get_(const int col) const
Definition: sparsemat.hpp:569
void EliminateRowColDiag(int rc, double value)
Perform elimination and set the diagonal entry to the given value.
Definition: sparsemat.cpp:1497
void SetSubMatrixTranspose(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2087
void Gauss_Seidel_forw(const Vector &x, Vector &y) const
Gauss-Seidel forward and backward iterations over a vector x.
Definition: sparsemat.cpp:1684
void Set(const int i, const int j, const double a)
Definition: sparsemat.cpp:2016
void EliminateRowCol(int rc, const double sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate row rc and column rc and modify the rhs using sol.
Definition: sparsemat.cpp:1204
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 Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
RowNode ** Rows
Array of linked lists, one for every row. This array represents the linked list (LIL) storage format...
Definition: sparsemat.hpp:61
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:644
void EliminateRow(int row, const double sol, Vector &rhs)
Eliminates a column from the transpose matrix.
Definition: sparsemat.cpp:1085
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
MemAlloc< RowNode, 1024 > RowNodeAlloc
Definition: sparsemat.hpp:68
int * GetJ() const
Return the array J.
Definition: sparsemat.hpp:136
virtual void EliminateZeroRows(const double threshold=1e-12)
If a row contains only zeros, set its diagonal to 1.
Definition: sparsemat.cpp:1657
double & operator()(int i, int j)
Returns reference to A[i][j].
Definition: sparsemat.cpp:400
void PrintInfo(std::ostream &out) const
Print various sparse matrix staticstics.
Definition: sparsemat.cpp:2597
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:1173
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const
Definition: sparsemat.cpp:1867
SparseMatrix & operator=(const SparseMatrix &rhs)
Assignment operator: deep copy.
Definition: sparsemat.cpp:187
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:2448
void SetDataOwner(bool owna)
Set the data ownership flag (A array).
Definition: sparsemat.hpp:422
void PrintCSR2(std::ostream &out) const
Prints a sparse matrix to stream out in CSR format.
Definition: sparsemat.cpp:2572
int MaxRowSize() const
Returns the maximum number of elements among all rows.
Definition: sparsemat.cpp:247
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
Definition: sparsemat.cpp:2210
Vector data type.
Definition: vector.hpp:48
int Walk(int &i, int &j, double &a)
Walks the sparse matrix.
void PartAddMult(const Array< int > &rows, const Vector &x, Vector &y, const double a=1.0) const
Definition: sparsemat.cpp:626
int CheckFinite() const
Count the number of entries that are NOT finite, i.e. Inf or Nan.
Definition: sparsemat.cpp:1056
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:2990
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:3260
double * GetData() const
Return element data, i.e. array A.
Definition: sparsemat.hpp:138
Type GetType() const
Definition: sparsemat.hpp:435
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25
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:569
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Definition: sparsemat.cpp:608
bool areColumnsSorted() const
Definition: sparsemat.hpp:311
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:667