MFEM  v4.5.2 Finite element discretization library
densemat.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
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_DENSEMAT
13 #define MFEM_DENSEMAT
14
15 #include "../config/config.hpp"
16 #include "../general/globals.hpp"
17 #include "matrix.hpp"
18
19 namespace mfem
20 {
21
22 /// Data type dense matrix using column-major storage
23 class DenseMatrix : public Matrix
24 {
25  friend class DenseTensor;
26  friend class DenseMatrixInverse;
27
28 private:
29  Memory<double> data;
30
31  void Eigensystem(Vector &ev, DenseMatrix *evect = NULL);
32
33  void Eigensystem(DenseMatrix &b, Vector &ev, DenseMatrix *evect = NULL);
34
35  // Auxiliary method used in FNorm2() and FNorm()
36  void FNorm(double &scale_factor, double &scaled_fnorm2) const;
37
38 public:
39  /** Default constructor for DenseMatrix.
40  Sets data = NULL and height = width = 0. */
41  DenseMatrix();
42
43  /// Copy constructor
44  DenseMatrix(const DenseMatrix &);
45
46  /// Creates square matrix of size s.
47  explicit DenseMatrix(int s);
48
49  /// Creates rectangular matrix of size m x n.
50  DenseMatrix(int m, int n);
51
52  /// Creates rectangular matrix equal to the transpose of mat.
53  DenseMatrix(const DenseMatrix &mat, char ch);
54
55  /// Construct a DenseMatrix using an existing data array.
56  /** The DenseMatrix does not assume ownership of the data array, i.e. it will
57  not delete the array. */
58  DenseMatrix(double *d, int h, int w)
59  : Matrix(h, w) { UseExternalData(d, h, w); }
60
61  /// Create a dense matrix using a braced initializer list
62  /// The inner lists correspond to rows of the matrix
63  template <int M, int N>
64  explicit DenseMatrix(const double (&values)[M][N]) : DenseMatrix(M, N)
65  {
66  // DenseMatrix is column-major so copies have to be element-wise
67  for (int i = 0; i < M; i++)
68  {
69  for (int j = 0; j < N; j++)
70  {
71  (*this)(i,j) = values[i][j];
72  }
73  }
74  }
75
76  /// Change the data array and the size of the DenseMatrix.
77  /** The DenseMatrix does not assume ownership of the data array, i.e. it will
78  not delete the data array @a d. This method should not be used with
79  DenseMatrix that owns its current data array. */
80  void UseExternalData(double *d, int h, int w)
81  {
82  data.Wrap(d, h*w, false);
83  height = h; width = w;
84  }
85
86  /// Change the data array and the size of the DenseMatrix.
87  /** The DenseMatrix does not assume ownership of the data array, i.e. it will
88  not delete the new array @a d. This method will delete the current data
89  array, if owned. */
90  void Reset(double *d, int h, int w)
91  { if (OwnsData()) { data.Delete(); } UseExternalData(d, h, w); }
92
93  /** Clear the data array and the dimensions of the DenseMatrix. This method
94  should not be used with DenseMatrix that owns its current data array. */
95  void ClearExternalData() { data.Reset(); height = width = 0; }
96
97  /// Delete the matrix data array (if owned) and reset the matrix state.
98  void Clear()
99  { if (OwnsData()) { data.Delete(); } ClearExternalData(); }
100
101  /// For backward compatibility define Size to be synonym of Width()
102  int Size() const { return Width(); }
103
104  /// Change the size of the DenseMatrix to s x s.
105  void SetSize(int s) { SetSize(s, s); }
106
107  /// Change the size of the DenseMatrix to h x w.
108  void SetSize(int h, int w);
109
110  /// Returns the matrix data array.
111  inline double *Data() const
112  { return const_cast<double*>((const double*)data);}
113
114  /// Returns the matrix data array.
115  inline double *GetData() const { return Data(); }
116
117  Memory<double> &GetMemory() { return data; }
118  const Memory<double> &GetMemory() const { return data; }
119
120  /// Return the DenseMatrix data (host pointer) ownership flag.
121  inline bool OwnsData() const { return data.OwnsHostPtr(); }
122
123  /// Returns reference to a_{ij}.
124  inline double &operator()(int i, int j);
125
126  /// Returns constant reference to a_{ij}.
127  inline const double &operator()(int i, int j) const;
128
129  /// Matrix inner product: tr(A^t B)
130  double operator*(const DenseMatrix &m) const;
131
132  /// Trace of a square matrix
133  double Trace() const;
134
135  /// Returns reference to a_{ij}.
136  virtual double &Elem(int i, int j);
137
138  /// Returns constant reference to a_{ij}.
139  virtual const double &Elem(int i, int j) const;
140
141  /// Matrix vector multiplication.
142  void Mult(const double *x, double *y) const;
143
144  /// Matrix vector multiplication.
145  void Mult(const double *x, Vector &y) const;
146
147  /// Matrix vector multiplication.
148  void Mult(const Vector &x, double *y) const;
149
150  /// Matrix vector multiplication.
151  virtual void Mult(const Vector &x, Vector &y) const;
152
153  /// Multiply a vector with the transpose matrix.
154  void MultTranspose(const double *x, double *y) const;
155
156  /// Multiply a vector with the transpose matrix.
157  void MultTranspose(const double *x, Vector &y) const;
158
159  /// Multiply a vector with the transpose matrix.
160  void MultTranspose(const Vector &x, double *y) const;
161
162  /// Multiply a vector with the transpose matrix.
163  virtual void MultTranspose(const Vector &x, Vector &y) const;
164
165  using Operator::Mult;
167
168  /// y += a * A.x
169  virtual void AddMult(const Vector &x, Vector &y, const double a = 1.0) const;
170
171  /// y += a * A^t x
172  virtual void AddMultTranspose(const Vector &x, Vector &y,
173  const double a = 1.0) const;
174
175  /// y += a * A.x
176  void AddMult_a(double a, const Vector &x, Vector &y) const;
177
178  /// y += a * A^t x
179  void AddMultTranspose_a(double a, const Vector &x, Vector &y) const;
180
181  /// Compute y^t A x
182  double InnerProduct(const double *x, const double *y) const;
183
184  /// LeftScaling this = diag(s) * this
185  void LeftScaling(const Vector & s);
186  /// InvLeftScaling this = diag(1./s) * this
187  void InvLeftScaling(const Vector & s);
188  /// RightScaling: this = this * diag(s);
189  void RightScaling(const Vector & s);
190  /// InvRightScaling: this = this * diag(1./s);
191  void InvRightScaling(const Vector & s);
192  /// SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
193  void SymmetricScaling(const Vector & s);
194  /// InvSymmetricScaling this = diag(sqrt(1./s)) * this * diag(sqrt(1./s))
195  void InvSymmetricScaling(const Vector & s);
196
197  /// Compute y^t A x
198  double InnerProduct(const Vector &x, const Vector &y) const
199  { return InnerProduct(x.GetData(), y.GetData()); }
200
201  /// Returns a pointer to the inverse matrix.
202  virtual MatrixInverse *Inverse() const;
203
204  /// Replaces the current matrix with its inverse
205  void Invert();
206
207  /// Replaces the current matrix with its square root inverse
208  void SquareRootInverse();
209
210  /// Calculates the determinant of the matrix
211  /// (optimized for 2x2, 3x3, and 4x4 matrices)
212  double Det() const;
213
214  double Weight() const;
215
216  /** @brief Set the matrix to alpha * A, assuming that A has the same
217  dimensions as the matrix and uses column-major layout. */
218  void Set(double alpha, const double *A);
219  /// Set the matrix to alpha * A.
220  void Set(double alpha, const DenseMatrix &A)
221  {
222  SetSize(A.Height(), A.Width());
223  Set(alpha, A.GetData());
224  }
225
226  /// Adds the matrix A multiplied by the number c to the matrix.
227  void Add(const double c, const DenseMatrix &A);
228
229  /// Adds the matrix A multiplied by the number c to the matrix,
230  /// assuming A has the same dimensions and uses column-major layout.
231  void Add(const double c, const double *A);
232
233  /// Sets the matrix elements equal to constant c
234  DenseMatrix &operator=(double c);
235
236  /// Copy the matrix entries from the given array
237  DenseMatrix &operator=(const double *d);
238
239  /// Sets the matrix size and elements equal to those of m
240  DenseMatrix &operator=(const DenseMatrix &m);
241
242  DenseMatrix &operator+=(const double *m);
244
246
247  DenseMatrix &operator*=(double c);
248
249  /// (*this) = -(*this)
250  void Neg();
251
252  /// Take the 2-norm of the columns of A and store in v
253  void Norm2(double *v) const;
254
255  /// Take the 2-norm of the columns of A and store in v
256  void Norm2(Vector &v) const
257  {
258  MFEM_ASSERT(v.Size() == Width(), "incompatible Vector size!");
259  Norm2(v.GetData());
260  }
261
262  /// Compute the norm ||A|| = max_{ij} |A_{ij}|
263  double MaxMaxNorm() const;
264
265  /// Compute the Frobenius norm of the matrix
266  double FNorm() const { double s, n2; FNorm(s, n2); return s*sqrt(n2); }
267
268  /// Compute the square of the Frobenius norm of the matrix
269  double FNorm2() const { double s, n2; FNorm(s, n2); return s*s*n2; }
270
271  /// Compute eigenvalues of A x = ev x where A = *this
272  void Eigenvalues(Vector &ev)
273  { Eigensystem(ev); }
274
275  /// Compute eigenvalues and eigenvectors of A x = ev x where A = *this
276  void Eigenvalues(Vector &ev, DenseMatrix &evect)
277  { Eigensystem(ev, &evect); }
278
279  /// Compute eigenvalues and eigenvectors of A x = ev x where A = *this
280  void Eigensystem(Vector &ev, DenseMatrix &evect)
281  { Eigensystem(ev, &evect); }
282
283  /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
284  where A = *this */
286  { Eigensystem(b, ev); }
287
288  /// Compute generalized eigenvalues of A x = ev B x, where A = *this
290  { Eigensystem(b, ev, &evect); }
291
292  /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
293  where A = *this */
295  { Eigensystem(b, ev, &evect); }
296
297  void SingularValues(Vector &sv) const;
298  int Rank(double tol) const;
299
300  /// Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
301  double CalcSingularvalue(const int i) const;
302
303  /** Return the eigenvalues (in increasing order) and eigenvectors of a
304  2x2 or 3x3 symmetric matrix. */
305  void CalcEigenvalues(double *lambda, double *vec) const;
306
307  void GetRow(int r, Vector &row) const;
308  void GetColumn(int c, Vector &col) const;
309  double *GetColumn(int col) { return data + col*height; }
310  const double *GetColumn(int col) const { return data + col*height; }
311
312  void GetColumnReference(int c, Vector &col)
313  { col.SetDataAndSize(data + c * height, height); }
314
315  void SetRow(int r, const double* row);
316  void SetRow(int r, const Vector &row);
317
318  void SetCol(int c, const double* col);
319  void SetCol(int c, const Vector &col);
320
321
322  /// Set all entries of a row to the specified value.
323  void SetRow(int row, double value);
324  /// Set all entries of a column to the specified value.
325  void SetCol(int col, double value);
326
327  /// Returns the diagonal of the matrix
328  void GetDiag(Vector &d) const;
329  /// Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|
330  void Getl1Diag(Vector &l) const;
331  /// Compute the row sums of the DenseMatrix
332  void GetRowSums(Vector &l) const;
333
334  /// Creates n x n diagonal matrix with diagonal elements c
335  void Diag(double c, int n);
336  /// Creates n x n diagonal matrix with diagonal given by diag
337  void Diag(double *diag, int n);
338
339  /// (*this) = (*this)^t
340  void Transpose();
341  /// (*this) = A^t
342  void Transpose(const DenseMatrix &A);
343  /// (*this) = 1/2 ((*this) + (*this)^t)
344  void Symmetrize();
345
346  void Lump();
347
348  /** Given a DShape matrix (from a scalar FE), stored in *this, returns the
349  CurlShape matrix. If *this is a N by D matrix, then curl is a D*N by
350  D*(D-1)/2 matrix. The size of curl must be set outside. The dimension D
351  can be either 2 or 3. In 2D this computes the scalar-valued curl of a
352  2D vector field */
354  /** Given a DShape matrix (from a scalar FE), stored in *this, returns the
355  CurlShape matrix. This computes the vector-valued curl of a scalar field.
356  *this is N by 2 matrix and curl is N by 2 matrix as well. */
358  /** Given a DShape matrix (from a scalar FE), stored in *this,
359  returns the DivShape vector. If *this is a N by dim matrix,
360  then div is a dim*N vector. The size of div must be set
361  outside. */
363
364  /// Copy rows row1 through row2 from A to *this
365  void CopyRows(const DenseMatrix &A, int row1, int row2);
366  /// Copy columns col1 through col2 from A to *this
367  void CopyCols(const DenseMatrix &A, int col1, int col2);
368  /// Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this
369  void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco);
370  /// Copy matrix A to the location in *this at row_offset, col_offset
371  void CopyMN(const DenseMatrix &A, int row_offset, int col_offset);
372  /// Copy matrix A^t to the location in *this at row_offset, col_offset
373  void CopyMNt(const DenseMatrix &A, int row_offset, int col_offset);
374  /** Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this at
375  row_offset, col_offset */
376  void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco,
377  int row_offset, int col_offset);
378  /// Copy c on the diagonal of size n to *this at row_offset, col_offset
379  void CopyMNDiag(double c, int n, int row_offset, int col_offset);
380  /// Copy diag on the diagonal of size n to *this at row_offset, col_offset
381  void CopyMNDiag(double *diag, int n, int row_offset, int col_offset);
382  /// Copy All rows and columns except m and n from A
383  void CopyExceptMN(const DenseMatrix &A, int m, int n);
384
385  /// Perform (ro+i,co+j)+=A(i,j) for 0<=i<A.Height, 0<=j<A.Width
386  void AddMatrix(DenseMatrix &A, int ro, int co);
387  /// Perform (ro+i,co+j)+=a*A(i,j) for 0<=i<A.Height, 0<=j<A.Width
388  void AddMatrix(double a, const DenseMatrix &A, int ro, int co);
389
390  /** Get the square submatrix which corresponds to the given indices @a idx.
391  Note: the @a A matrix will be resized to accommodate the data */
392  void GetSubMatrix(const Array<int> & idx, DenseMatrix & A) const;
393
394  /** Get the rectangular submatrix which corresponds to the given indices
395  @a idx_i and @a idx_j. Note: the @a A matrix will be resized to
396  accommodate the data */
397  void GetSubMatrix(const Array<int> & idx_i, const Array<int> & idx_j,
398  DenseMatrix & A) const;
399
400  /** Get the square submatrix which corresponds to the range
401  [ @a ibeg, @a iend ). Note: the @a A matrix will be resized
402  to accommodate the data */
403  void GetSubMatrix(int ibeg, int iend, DenseMatrix & A);
404
405  /** Get the square submatrix which corresponds to the range
406  i âˆˆ [ @a ibeg, @a iend ) and j âˆˆ [ @a jbeg, @a jend )
407  Note: the @a A matrix will be resized to accommodate the data */
408  void GetSubMatrix(int ibeg, int iend, int jbeg, int jend, DenseMatrix & A);
409
410  /// Set (*this)(idx[i],idx[j]) = A(i,j)
411  void SetSubMatrix(const Array<int> & idx, const DenseMatrix & A);
412
413  /// Set (*this)(idx_i[i],idx_j[j]) = A(i,j)
414  void SetSubMatrix(const Array<int> & idx_i, const Array<int> & idx_j,
415  const DenseMatrix & A);
416
417  /** Set a submatrix of (this) to the given matrix @a A
418  with row and column offset @a ibeg */
419  void SetSubMatrix(int ibeg, const DenseMatrix & A);
420
421  /** Set a submatrix of (this) to the given matrix @a A
422  with row and column offset @a ibeg and @a jbeg respectively */
423  void SetSubMatrix(int ibeg, int jbeg, const DenseMatrix & A);
424
425  /// (*this)(idx[i],idx[j]) += A(i,j)
426  void AddSubMatrix(const Array<int> & idx, const DenseMatrix & A);
427
428  /// (*this)(idx_i[i],idx_j[j]) += A(i,j)
429  void AddSubMatrix(const Array<int> & idx_i, const Array<int> & idx_j,
430  const DenseMatrix & A);
431
432  /** Add the submatrix @a A to this with row and column offset @a ibeg */
433  void AddSubMatrix(int ibeg, const DenseMatrix & A);
434
435  /** Add the submatrix @a A to this with row and column offsets
436  @a ibeg and @a jbeg respectively */
437  void AddSubMatrix(int ibeg, int jbeg, const DenseMatrix & A);
438
439  /// Add the matrix 'data' to the Vector 'v' at the given 'offset'
440  void AddToVector(int offset, Vector &v) const;
441  /// Get the matrix 'data' from the Vector 'v' at the given 'offset'
442  void GetFromVector(int offset, const Vector &v);
443  /** If (dofs[i] < 0 and dofs[j] >= 0) or (dofs[i] >= 0 and dofs[j] < 0)
444  then (*this)(i,j) = -(*this)(i,j). */
446
447  /// Replace small entries, abs(a_ij) <= eps, with zero.
448  void Threshold(double eps);
449
450  /** Count the number of entries in the matrix for which isfinite
451  is false, i.e. the entry is a NaN or +/-Inf. */
452  int CheckFinite() const { return mfem::CheckFinite(HostRead(), height*width); }
453
454  /// Prints matrix to stream out.
455  virtual void Print(std::ostream &out = mfem::out, int width_ = 4) const;
456  virtual void PrintMatlab(std::ostream &out = mfem::out) const;
457  /// Prints the transpose matrix to stream out.
458  virtual void PrintT(std::ostream &out = mfem::out, int width_ = 4) const;
459
460  /// Invert and print the numerical conditioning of the inversion.
461  void TestInversion();
462
463  std::size_t MemoryUsage() const { return data.Capacity() * sizeof(double); }
464
465  /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
466  const double *Read(bool on_dev = true) const
467  { return mfem::Read(data, Height()*Width(), on_dev); }
468
469  /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
471  { return mfem::Read(data, Height()*Width(), false); }
472
473  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
474  double *Write(bool on_dev = true)
475  { return mfem::Write(data, Height()*Width(), on_dev); }
476
477  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
478  double *HostWrite()
479  { return mfem::Write(data, Height()*Width(), false); }
480
481  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
482  double *ReadWrite(bool on_dev = true)
483  { return mfem::ReadWrite(data, Height()*Width(), on_dev); }
484
485  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
487  { return mfem::ReadWrite(data, Height()*Width(), false); }
488
489  void Swap(DenseMatrix &other);
490
491  /// Destroys dense matrix.
492  virtual ~DenseMatrix();
493 };
494
495 /// C = A + alpha*B
496 void Add(const DenseMatrix &A, const DenseMatrix &B,
497  double alpha, DenseMatrix &C);
498
499 /// C = alpha*A + beta*B
500 void Add(double alpha, const double *A,
501  double beta, const double *B, DenseMatrix &C);
502
503 /// C = alpha*A + beta*B
504 void Add(double alpha, const DenseMatrix &A,
505  double beta, const DenseMatrix &B, DenseMatrix &C);
506
507 /// @brief Solves the dense linear system, A * X = B for X
508 ///
509 /// @param [in,out] A the square matrix for the linear system
510 /// @param [in,out] X the rhs vector, B, on input, the solution, X, on output.
511 /// @param [in] TOL optional fuzzy comparison tolerance. Defaults to 1e-9.
512 ///
513 /// @return status set to true if successful, otherwise, false.
514 ///
515 /// @note This routine may replace the contents of the input Matrix, A, with the
516 /// corresponding LU factorization of the matrix. Matrices of size 1x1 and
517 /// 2x2 are handled explicitly.
518 ///
519 /// @pre A.IsSquare() == true
520 /// @pre X != nullptr
521 bool LinearSolve(DenseMatrix& A, double* X, double TOL = 1.e-9);
522
523 /// Matrix matrix multiplication. A = B * C.
524 void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);
525
526 /// Matrix matrix multiplication. A += B * C.
527 void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);
528
529 /// Matrix matrix multiplication. A += alpha * B * C.
530 void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c,
531  DenseMatrix &a);
532
533 /** Calculate the adjugate of a matrix (for NxN matrices, N=1,2,3) or the matrix
534  adj(A^t.A).A^t for rectangular matrices (2x1, 3x1, or 3x2). This operation
535  is well defined even when the matrix is not full rank. */
537
538 /// Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
540
541 /** Calculate the inverse of a matrix (for NxN matrices, N=1,2,3) or the
542  left inverse (A^t.A)^{-1}.A^t (for 2x1, 3x1, or 3x2 matrices) */
543 void CalcInverse(const DenseMatrix &a, DenseMatrix &inva);
544
545 /// Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
546 void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva);
547
548 /** For a given Nx(N-1) (N=2,3) matrix J, compute a vector n such that
549  n_k = (-1)^{k+1} det(J_k), k=1,..,N, where J_k is the matrix J with the
550  k-th row removed. Note: J^t.n = 0, det([n|J])=|n|^2=det(J^t.J). */
551 void CalcOrtho(const DenseMatrix &J, Vector &n);
552
553 /// Calculate the matrix A.At
554 void MultAAt(const DenseMatrix &a, DenseMatrix &aat);
555
556 /// ADAt = A D A^t, where D is diagonal
558
559 /// ADAt += A D A^t, where D is diagonal
561
562 /// Multiply a matrix A with the transpose of a matrix B: A*Bt
563 void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);
564
565 /// ADBt = A D B^t, where D is diagonal
566 void MultADBt(const DenseMatrix &A, const Vector &D,
567  const DenseMatrix &B, DenseMatrix &ADBt);
568
569 /// ABt += A * B^t
570 void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);
571
572 /// ADBt = A D B^t, where D is diagonal
574  const DenseMatrix &B, DenseMatrix &ADBt);
575
576 /// ABt += a * A * B^t
577 void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B,
578  DenseMatrix &ABt);
579
580 /// Multiply the transpose of a matrix A with a matrix B: At*B
581 void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB);
582
583 /// AAt += a * A * A^t
584 void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);
585
586 /// AAt = a * A * A^t
587 void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);
588
589 /// Make a matrix from a vector V.Vt
590 void MultVVt(const Vector &v, DenseMatrix &vvt);
591
592 void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);
593
594 /// VWt += v w^t
595 void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);
596
597 /// VVt += v v^t
598 void AddMultVVt(const Vector &v, DenseMatrix &VWt);
599
600 /// VWt += a * v w^t
601 void AddMult_a_VWt(const double a, const Vector &v, const Vector &w,
602  DenseMatrix &VWt);
603
604 /// VVt += a * v v^t
605 void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt);
606
607 /** Computes matrix P^t * A * P. Note: The @a RAP matrix will be resized
608  to accommodate the data */
609 void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix & RAP);
610
611 /** Computes the matrix Rt^t * A * P. Note: The @a RAP matrix will be resized
612  to accommodate the data */
613 void RAP(const DenseMatrix &Rt, const DenseMatrix &A,
614  const DenseMatrix &P, DenseMatrix & RAP);
615
616 /** Abstract class that can compute factorization of external data and perform various
617  operations with the factored data. */
618 class Factors
619 {
620 public:
621
622  double *data;
623
624  Factors() { }
625
626  Factors(double *data_) : data(data_) { }
627
628  virtual bool Factor(int m, double TOL = 0.0)
629  {
630  mfem_error("Factors::Factors(...)");
631  return false;
632  }
633
634  virtual double Det(int m) const
635  {
636  mfem_error("Factors::Det(...)");
637  return 0.;
638  }
639
640  virtual void Solve(int m, int n, double *X) const
641  {
642  mfem_error("Factors::Solve(...)");
643  }
644
645  virtual void GetInverseMatrix(int m, double *X) const
646  {
647  mfem_error("Factors::GetInverseMatrix(...)");
648  }
649
650  virtual ~Factors() {}
651 };
652
653
654 /** Class that can compute LU factorization of external data and perform various
655  operations with the factored data. */
656 class LUFactors : public Factors
657 {
658 public:
659  int *ipiv;
660 #ifdef MFEM_USE_LAPACK
661  static const int ipiv_base = 1;
662 #else
663  static const int ipiv_base = 0;
664 #endif
665
666  /** With this constructor, the (public) data and ipiv members should be set
667  explicitly before calling class methods. */
669
670  LUFactors(double *data_, int *ipiv_) : Factors(data_), ipiv(ipiv_) { }
671
672  /**
673  * @brief Compute the LU factorization of the current matrix
674  *
675  * Factorize the current matrix of size (m x m) overwriting it with the
676  * LU factors. The factorization is such that L.U = P.A, where A is the
677  * original matrix and P is a permutation matrix represented by ipiv.
678  *
679  * @param [in] m size of the square matrix
680  * @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0.
681  *
682  * @return status set to true if successful, otherwise, false.
683  */
684  virtual bool Factor(int m, double TOL = 0.0);
685
686  /** Assuming L.U = P.A factored data of size (m x m), compute |A|
687  from the diagonal values of U and the permutation information. */
688  virtual double Det(int m) const;
689
690  /** Assuming L.U = P.A factored data of size (m x m), compute X <- A X,
691  for a matrix X of size (m x n). */
692  void Mult(int m, int n, double *X) const;
693
694  /** Assuming L.U = P.A factored data of size (m x m), compute
695  X <- L^{-1} P X, for a matrix X of size (m x n). */
696  void LSolve(int m, int n, double *X) const;
697
698  /** Assuming L.U = P.A factored data of size (m x m), compute
699  X <- U^{-1} X, for a matrix X of size (m x n). */
700  void USolve(int m, int n, double *X) const;
701
702  /** Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1} X,
703  for a matrix X of size (m x n). */
704  virtual void Solve(int m, int n, double *X) const;
705
706  /** Assuming L.U = P.A factored data of size (m x m), compute X <- X A^{-1},
707  for a matrix X of size (n x m). */
708  void RightSolve(int m, int n, double *X) const;
709
710  /// Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1}.
711  virtual void GetInverseMatrix(int m, double *X) const;
712
713  /** Given an (n x m) matrix A21, compute X2 <- X2 - A21 X1, for matrices X1,
714  and X2 of size (m x r) and (n x r), respectively. */
715  static void SubMult(int m, int n, int r, const double *A21,
716  const double *X1, double *X2);
717
718  /** Assuming P.A = L.U factored data of size (m x m), compute the 2x2 block
719  decomposition:
720  | P 0 | | A A12 | = | L 0 | | U U12 |
721  | 0 I | | A21 A22 | | L21 I | | 0 S22 |
722  where A12, A21, and A22 are matrices of size (m x n), (n x m), and
723  (n x n), respectively. The blocks are overwritten as follows:
724  A12 <- U12 = L^{-1} P A12
725  A21 <- L21 = A21 U^{-1}
726  A22 <- S22 = A22 - L21 U12.
727  The block S22 is the Schur complement. */
728  void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const;
729
730  /** Given BlockFactor()'d data, perform the forward block solve for the
731  linear system:
732  | A A12 | | X1 | = | B1 |
733  | A21 A22 | | X2 | | B2 |
734  written in the factored form:
735  | L 0 | | U U12 | | X1 | = | P 0 | | B1 |
736  | L21 I | | 0 S22 | | X2 | | 0 I | | B2 |.
737  The resulting blocks Y1, Y2 solve the system:
738  | L 0 | | Y1 | = | P 0 | | B1 |
739  | L21 I | | Y2 | | 0 I | | B2 |
740  The blocks are overwritten as follows:
741  B1 <- Y1 = L^{-1} P B1
742  B2 <- Y2 = B2 - L21 Y1 = B2 - A21 A^{-1} B1
743  The blocks B1/Y1 and B2/Y2 are of size (m x r) and (n x r), respectively.
744  The Schur complement system is given by: S22 X2 = Y2. */
745  void BlockForwSolve(int m, int n, int r, const double *L21,
746  double *B1, double *B2) const;
747
748  /** Given BlockFactor()'d data, perform the backward block solve in
749  | U U12 | | X1 | = | Y1 |
750  | 0 S22 | | X2 | | Y2 |.
751  The input is the solution block X2 and the block Y1 resulting from
752  BlockForwSolve(). The result block X1 overwrites input block Y1:
753  Y1 <- X1 = U^{-1} (Y1 - U12 X2). */
754  void BlockBackSolve(int m, int n, int r, const double *U12,
755  const double *X2, double *Y1) const;
756 };
757
758
759 /** Class that can compute Cholesky factorizations of external data of an
760  SPD matrix and perform various operations with the factored data. */
761 class CholeskyFactors : public Factors
762 {
763 public:
764
765  /** With this constructor, the (public) data should be set
766  explicitly before calling class methods. */
768
769  CholeskyFactors(double *data_) : Factors(data_) { }
770
771  /**
772  * @brief Compute the Cholesky factorization of the current matrix
773  *
774  * Factorize the current matrix of size (m x m) overwriting it with the
775  * Cholesky factors. The factorization is such that LL^t = A, where A is the
776  * original matrix
777  *
778  * @param [in] m size of the square matrix
779  * @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0.
780  *
781  * @return status set to true if successful, otherwise, false.
782  */
783  virtual bool Factor(int m, double TOL = 0.0);
784
785  /** Assuming LL^t = A factored data of size (m x m), compute |A|
786  from the diagonal values of L */
787  virtual double Det(int m) const;
788
789  /** Assuming L.L^t = A factored data of size (m x m), compute X <- L X,
790  for a matrix X of size (m x n). */
791  void LMult(int m, int n, double *X) const;
792
793  /** Assuming L.L^t = A factored data of size (m x m), compute X <- L^t X,
794  for a matrix X of size (m x n). */
795  void UMult(int m, int n, double *X) const;
796
797  /** Assuming L L^t = A factored data of size (m x m), compute
798  X <- L^{-1} X, for a matrix X of size (m x n). */
799  void LSolve(int m, int n, double *X) const;
800
801  /** Assuming L L^t = A factored data of size (m x m), compute
802  X <- L^{-t} X, for a matrix X of size (m x n). */
803  void USolve(int m, int n, double *X) const;
804
805  /** Assuming L.L^t = A factored data of size (m x m), compute X <- A^{-1} X,
806  for a matrix X of size (m x n). */
807  virtual void Solve(int m, int n, double *X) const;
808
809  /** Assuming L.L^t = A factored data of size (m x m), compute X <- X A^{-1},
810  for a matrix X of size (n x m). */
811  void RightSolve(int m, int n, double *X) const;
812
813  /// Assuming L.L^t = A factored data of size (m x m), compute X <- A^{-1}.
814  virtual void GetInverseMatrix(int m, double *X) const;
815
816 };
817
818
819 /** Data type for inverse of square dense matrix.
820  Stores matrix factors, i.e., Cholesky factors if the matrix is SPD,
821  LU otherwise. */
823 {
824 private:
825  const DenseMatrix *a;
826  Factors * factors = nullptr;
827  bool spd = false;
828
829  void Init(int m);
830  bool own_data = false;
831 public:
832  /// Default constructor.
833  DenseMatrixInverse(bool spd_=false) : a(NULL), spd(spd_) { Init(0); }
834
835  /** Creates square dense matrix. Computes factorization of mat
836  and stores its factors. */
837  DenseMatrixInverse(const DenseMatrix &mat, bool spd_ = false);
838
839  /// Same as above but does not factorize the matrix.
840  DenseMatrixInverse(const DenseMatrix *mat, bool spd_ = false);
841
842  /// Get the size of the inverse matrix
843  int Size() const { return Width(); }
844
845  /// Factor the current DenseMatrix, *a
846  void Factor();
847
848  /// Factor a new DenseMatrix of the same size
849  void Factor(const DenseMatrix &mat);
850
851  virtual void SetOperator(const Operator &op);
852
853  /// Matrix vector multiplication with the inverse of dense matrix.
854  void Mult(const double *x, double *y) const;
855
856  /// Matrix vector multiplication with the inverse of dense matrix.
857  virtual void Mult(const Vector &x, Vector &y) const;
858
859  /// Multiply the inverse matrix by another matrix: X = A^{-1} B.
860  void Mult(const DenseMatrix &B, DenseMatrix &X) const;
861
862  /// Multiply the inverse matrix by another matrix: X <- A^{-1} X.
863  void Mult(DenseMatrix &X) const {factors->Solve(width, X.Width(), X.Data());}
864
865  using Operator::Mult;
866
867  /// Compute and return the inverse matrix in Ainv.
868  void GetInverseMatrix(DenseMatrix &Ainv) const;
869
870  /// Compute the determinant of the original DenseMatrix using the LU factors.
871  double Det() const { return factors->Det(width); }
872
873  /// Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
874  void TestInversion();
875
876  /// Destroys dense inverse matrix.
877  virtual ~DenseMatrixInverse();
878 };
879
880 #ifdef MFEM_USE_LAPACK
881
883 {
884  DenseMatrix &mat;
885  Vector EVal;
886  DenseMatrix EVect;
887  Vector ev;
888  int n;
889  double *work;
890  char jobz, uplo;
891  int lwork, info;
892 public:
893
896  void Eval();
897  Vector &Eigenvalues() { return EVal; }
898  DenseMatrix &Eigenvectors() { return EVect; }
899  double Eigenvalue(int i) { return EVal(i); }
900  const Vector &Eigenvector(int i)
901  {
902  ev.SetData(EVect.Data() + i * EVect.Height());
903  return ev;
904  }
906 };
907
909 {
910  DenseMatrix &A;
911  DenseMatrix &B;
912  DenseMatrix A_copy;
913  DenseMatrix B_copy;
914  Vector evalues_r;
915  Vector evalues_i;
916  DenseMatrix Vr;
917  DenseMatrix Vl;
918  int n;
919
920  double *alphar;
921  double *alphai;
922  double *beta;
923  double *work;
924  char jobvl, jobvr;
925  int lwork, info;
926
927 public:
928
930  bool left_eigen_vectors = false,
931  bool right_eigen_vectors = false);
932  void Eval();
933  Vector &EigenvaluesRealPart() { return evalues_r; }
934  Vector &EigenvaluesImagPart() { return evalues_i; }
935  double EigenvalueRealPart(int i) { return evalues_r(i); }
936  double EigenvalueImagPart(int i) { return evalues_i(i); }
937  DenseMatrix &LeftEigenvectors() { return Vl; }
938  DenseMatrix &RightEigenvectors() { return Vr; }
940 };
941
942
944 {
945  DenseMatrix Mc;
946  Vector sv;
947  DenseMatrix U,Vt;
948  int m, n;
949
950 #ifdef MFEM_USE_LAPACK
951  double *work;
952  char jobu, jobvt;
953  int lwork, info;
954 #endif
955
956  void Init();
957 public:
959  bool left_singular_vectors=false,
960  bool right_singlular_vectors=false);
961  DenseMatrixSVD(int h, int w,
962  bool left_singular_vectors=false,
963  bool right_singlular_vectors=false);
964  void Eval(DenseMatrix &M);
965  Vector &Singularvalues() { return sv; }
966  double Singularvalue(int i) { return sv(i); }
969  ~DenseMatrixSVD();
970 };
971
972 #endif // if MFEM_USE_LAPACK
973
974
975 class Table;
976
977 /// Rank 3 tensor (array of matrices)
979 {
980 private:
981  mutable DenseMatrix Mk;
982  Memory<double> tdata;
983  int nk;
984
985 public:
987  {
988  nk = 0;
989  }
990
991  DenseTensor(int i, int j, int k)
992  : Mk(NULL, i, j)
993  {
994  nk = k;
995  tdata.New(i*j*k);
996  }
997
998  DenseTensor(double *d, int i, int j, int k)
999  : Mk(NULL, i, j)
1000  {
1001  nk = k;
1002  tdata.Wrap(d, i*j*k, false);
1003  }
1004
1005  DenseTensor(int i, int j, int k, MemoryType mt)
1006  : Mk(NULL, i, j)
1007  {
1008  nk = k;
1009  tdata.New(i*j*k, mt);
1010  }
1011
1012  /// Copy constructor: deep copy
1013  DenseTensor(const DenseTensor &other)
1014  : Mk(NULL, other.Mk.height, other.Mk.width), nk(other.nk)
1015  {
1016  const int size = Mk.Height()*Mk.Width()*nk;
1017  if (size > 0)
1018  {
1019  tdata.New(size, other.tdata.GetMemoryType());
1020  tdata.CopyFrom(other.tdata, size);
1021  }
1022  }
1023
1024  int SizeI() const { return Mk.Height(); }
1025  int SizeJ() const { return Mk.Width(); }
1026  int SizeK() const { return nk; }
1027
1028  int TotalSize() const { return SizeI()*SizeJ()*SizeK(); }
1029
1030  void SetSize(int i, int j, int k, MemoryType mt_ = MemoryType::PRESERVE)
1031  {
1032  const MemoryType mt = mt_ == MemoryType::PRESERVE ? tdata.GetMemoryType() : mt_;
1033  tdata.Delete();
1034  Mk.UseExternalData(NULL, i, j);
1035  nk = k;
1036  tdata.New(i*j*k, mt);
1037  }
1038
1039  void UseExternalData(double *ext_data, int i, int j, int k)
1040  {
1041  tdata.Delete();
1042  Mk.UseExternalData(NULL, i, j);
1043  nk = k;
1044  tdata.Wrap(ext_data, i*j*k, false);
1045  }
1046
1047  /// Sets the tensor elements equal to constant c
1048  DenseTensor &operator=(double c);
1049
1050  /// Copy assignment operator (performs a deep copy)
1051  DenseTensor &operator=(const DenseTensor &other);
1052
1054  {
1055  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1056  Mk.data = Memory<double>(GetData(k), SizeI()*SizeJ(), false);
1057  return Mk;
1058  }
1059  const DenseMatrix &operator()(int k) const
1060  {
1061  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1062  Mk.data = Memory<double>(const_cast<double*>(GetData(k)), SizeI()*SizeJ(),
1063  false);
1064  return Mk;
1065  }
1066
1067  double &operator()(int i, int j, int k)
1068  {
1069  MFEM_ASSERT_INDEX_IN_RANGE(i, 0, SizeI());
1070  MFEM_ASSERT_INDEX_IN_RANGE(j, 0, SizeJ());
1071  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1072  return tdata[i+SizeI()*(j+SizeJ()*k)];
1073  }
1074
1075  const double &operator()(int i, int j, int k) const
1076  {
1077  MFEM_ASSERT_INDEX_IN_RANGE(i, 0, SizeI());
1078  MFEM_ASSERT_INDEX_IN_RANGE(j, 0, SizeJ());
1079  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1080  return tdata[i+SizeI()*(j+SizeJ()*k)];
1081  }
1082
1083  double *GetData(int k)
1084  {
1085  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1086  return tdata+k*Mk.Height()*Mk.Width();
1087  }
1088
1089  const double *GetData(int k) const
1090  {
1091  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1092  return tdata+k*Mk.Height()*Mk.Width();
1093  }
1094
1095  double *Data() { return tdata; }
1096
1097  const double *Data() const { return tdata; }
1098
1099  Memory<double> &GetMemory() { return tdata; }
1100  const Memory<double> &GetMemory() const { return tdata; }
1101
1102  /** Matrix-vector product from unassembled element matrices, assuming both
1103  'x' and 'y' use the same elem_dof table. */
1104  void AddMult(const Table &elem_dof, const Vector &x, Vector &y) const;
1105
1106  void Clear()
1107  { UseExternalData(NULL, 0, 0, 0); }
1108
1109  std::size_t MemoryUsage() const { return nk*Mk.MemoryUsage(); }
1110
1111  /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
1112  const double *Read(bool on_dev = true) const
1113  { return mfem::Read(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
1114
1115  /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
1117  { return mfem::Read(tdata, Mk.Height()*Mk.Width()*nk, false); }
1118
1119  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
1120  double *Write(bool on_dev = true)
1121  { return mfem::Write(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
1122
1123  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
1124  double *HostWrite()
1125  { return mfem::Write(tdata, Mk.Height()*Mk.Width()*nk, false); }
1126
1127  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
1128  double *ReadWrite(bool on_dev = true)
1129  { return mfem::ReadWrite(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
1130
1131  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
1133  { return mfem::ReadWrite(tdata, Mk.Height()*Mk.Width()*nk, false); }
1134
1136  {
1137  mfem::Swap(tdata, t.tdata);
1138  mfem::Swap(nk, t.nk);
1139  Mk.Swap(t.Mk);
1140  }
1141
1142  ~DenseTensor() { tdata.Delete(); }
1143 };
1144
1145 /** @brief Compute the LU factorization of a batch of matrices
1146
1147  Factorize n matrices of size (m x m) stored in a dense tensor overwriting it
1148  with the LU factors. The factorization is such that L.U = Piv.A, where A is
1149  the original matrix and Piv is a permutation matrix represented by P.
1150
1151  @param [in, out] Mlu batch of square matrices - dimension m x m x n.
1152  @param [out] P array storing pivot information - dimension m x n.
1153  @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0. */
1154 void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const double TOL = 0.0);
1155
1156 /** @brief Solve batch linear systems
1157
1158  Assuming L.U = P.A for n factored matrices (m x m), compute x <- A x, for n
1159  companion vectors.
1160
1161  @param [in] Mlu batch of LU factors for matrix M - dimension m x m x n.
1162  @param [in] P array storing pivot information - dimension m x n.
1163  @param [in, out] X vector storing right-hand side and then solution -
1164  dimension m x n. */
1165 void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X);
1166
1167
1168 // Inline methods
1169
1170 inline double &DenseMatrix::operator()(int i, int j)
1171 {
1172  MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
1173  return data[i+j*height];
1174 }
1175
1176 inline const double &DenseMatrix::operator()(int i, int j) const
1177 {
1178  MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
1179  return data[i+j*height];
1180 }
1181
1182 } // namespace mfem
1183
1184 #endif
const double * GetColumn(int col) const
Definition: densemat.hpp:310
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
Definition: densemat.cpp:1453
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:93
int CheckFinite() const
Definition: densemat.hpp:452
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
Definition: densemat.cpp:2781
DenseMatrix & operator-=(const DenseMatrix &m)
Definition: densemat.cpp:647
void SymmetricScaling(const Vector &s)
SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
Definition: densemat.cpp:409
void SquareRootInverse()
Replaces the current matrix with its square root inverse.
Definition: densemat.cpp:789
void RightSolve(int m, int n, double *X) const
Definition: densemat.cpp:3434
void RightSolve(int m, int n, double *X) const
Definition: densemat.cpp:3777
int CheckFinite(const double *v, const int n)
Definition: vector.hpp:492
void BatchLUFactor(DenseTensor &Mlu, Array< int > &P, const double TOL)
Compute the LU factorization of a batch of matrices.
Definition: densemat.cpp:4299
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:3166
void Mult(DenseMatrix &X) const
Multiply the inverse matrix by another matrix: X <- A^{-1} X.
Definition: densemat.hpp:863
DenseMatrix & operator*=(double c)
Definition: densemat.cpp:660
void SetCol(int c, const double *col)
Definition: densemat.cpp:2193
void UseExternalData(double *ext_data, int i, int j, int k)
Definition: densemat.hpp:1039
DenseTensor & operator=(double c)
Sets the tensor elements equal to constant c.
Definition: densemat.cpp:4282
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:3147
double * HostWrite()
Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:478
Memory< double > & GetMemory()
Definition: densemat.hpp:117
void SetRow(int r, const double *row)
Definition: densemat.cpp:2178
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
Definition: densemat.cpp:395
const DenseMatrix & operator()(int k) const
Definition: densemat.hpp:1059
virtual ~Factors()
Definition: densemat.hpp:650
void Mult(int m, int n, double *X) const
Definition: densemat.cpp:3345
void Eigenvalues(Vector &ev)
Compute eigenvalues of A x = ev x where A = *this.
Definition: densemat.hpp:272
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3963
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:3929
virtual double Det(int m) const
Definition: densemat.cpp:3328
void Delete()
Delete the owned pointers and reset the Memory object.
DenseMatrix & operator()(int k)
Definition: densemat.hpp:1053
void Eigensystem(DenseMatrix &b, Vector &ev, DenseMatrix &evect)
Definition: densemat.hpp:294
void AddMultTranspose_a(double a, const Vector &x, Vector &y) const
y += a * A^t x
Definition: densemat.cpp:317
Definition: densemat.hpp:1116
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
Definition: densemat.cpp:2479
DenseMatrixSVD(DenseMatrix &M, bool left_singular_vectors=false, bool right_singlular_vectors=false)
Definition: densemat.cpp:4154
double MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
Definition: densemat.cpp:847
double Det() const
Definition: densemat.cpp:488
const Memory< double > & GetMemory() const
Definition: densemat.hpp:1100
virtual void GetInverseMatrix(int m, double *X) const
Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1}.
Definition: densemat.cpp:3489
virtual void Solve(int m, int n, double *X) const
Definition: densemat.cpp:3420
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
T * Write(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for write access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true...
Definition: device.hpp:336
void TestInversion()
Invert and print the numerical conditioning of the inversion.
Definition: densemat.cpp:2293
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
DenseMatrix & LeftSingularvectors()
Definition: densemat.hpp:967
void CopyRows(const DenseMatrix &A, int row1, int row2)
Copy rows row1 through row2 from A to *this.
Definition: densemat.cpp:1566
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
void Swap(DenseTensor &t)
Definition: densemat.hpp:1135
void LSolve(int m, int n, double *X) const
Definition: densemat.cpp:3703
int SizeJ() const
Definition: densemat.hpp:1025
void Eval(DenseMatrix &M)
Definition: densemat.cpp:4188
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
Definition: densemat.cpp:1268
Abstract data type for matrix inverse.
Definition: matrix.hpp:62
double Det() const
Compute the determinant of the original DenseMatrix using the LU factors.
Definition: densemat.hpp:871
void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += a * A * B^t.
Definition: densemat.cpp:2979
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Definition: densemat.hpp:486
Factors(double *data_)
Definition: densemat.hpp:626
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3917
virtual void Solve(int m, int n, double *X) const
Definition: densemat.cpp:3763
void CopyFrom(const Memory &src, int size)
Copy size entries from src to *this.
Definition: densemat.cpp:1537
double Weight() const
Definition: densemat.cpp:545
const double * Data() const
Definition: densemat.hpp:1097
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2693
DenseMatrix & operator=(double c)
Sets the matrix elements equal to constant c.
Definition: densemat.cpp:601
void Set(double alpha, const double *A)
Set the matrix to alpha * A, assuming that A has the same dimensions as the matrix and uses column-ma...
Definition: densemat.cpp:572
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:464
void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt = a * A * A^t.
Definition: densemat.cpp:3120
void BlockForwSolve(int m, int n, int r, const double *L21, double *B1, double *B2) const
Definition: densemat.cpp:3595
void GetRowSums(Vector &l) const
Compute the row sums of the DenseMatrix.
Definition: densemat.cpp:1376
virtual double Det(int m) const
Definition: densemat.cpp:3656
const Memory< double > & GetMemory() const
Definition: densemat.hpp:118
void Getl1Diag(Vector &l) const
Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|.
Definition: densemat.cpp:1359
virtual double Det(int m) const
Definition: densemat.hpp:634
static void SubMult(int m, int n, int r, const double *A21, const double *X1, double *X2)
Definition: densemat.cpp:3552
virtual void GetInverseMatrix(int m, double *X) const
Definition: densemat.hpp:645
virtual void Solve(int m, int n, double *X) const
Definition: densemat.hpp:640
virtual void PrintMatlab(std::ostream &out=mfem::out) const
Prints operator in Matlab format.
Definition: densemat.cpp:2248
const Vector & Eigenvector(int i)
Definition: densemat.hpp:900
void USolve(int m, int n, double *X) const
Definition: densemat.cpp:3402
virtual void PrintT(std::ostream &out=mfem::out, int width_=4) const
Prints the transpose matrix to stream out.
Definition: densemat.cpp:2267
void Set(double alpha, const DenseMatrix &A)
Set the matrix to alpha * A.
Definition: densemat.hpp:220
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2321
double & operator()(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.hpp:1170
void Wrap(T *ptr, int size, bool own)
Wrap an externally allocated host pointer, ptr with the current host memory type returned by MemoryMa...
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:1112
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2551
double & operator()(int i, int j, int k)
Definition: densemat.hpp:1067
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2447
void Norm2(double *v) const
Take the 2-norm of the columns of A and store in v.
Definition: densemat.cpp:834
double Singularvalue(int i)
Definition: densemat.hpp:966
void Reset(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:90
int Size() const
For backward compatibility define Size to be synonym of Width()
Definition: densemat.hpp:102
DenseMatrix(const double(&values)[M][N])
Definition: densemat.hpp:64
Definition: densemat.hpp:470
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:581
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:215
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
void AddMult_a_VWt(const double a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
Definition: densemat.cpp:3211
double b
Definition: lissajous.cpp:42
void InvSymmetricScaling(const Vector &s)
InvSymmetricScaling this = diag(sqrt(1./s)) * this * diag(sqrt(1./s))
Definition: densemat.cpp:437
Abstract data type matrix.
Definition: matrix.hpp:27
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
Definition: densemat.cpp:336
virtual void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: densemat.cpp:2222
void MultADBt(const DenseMatrix &A, const Vector &D, const DenseMatrix &B, DenseMatrix &ADBt)
ADBt = A D B^t, where D is diagonal.
Definition: densemat.cpp:2846
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:679
bool LinearSolve(DenseMatrix &A, double *X, double TOL)
Solves the dense linear system, A * X = B for X
Definition: densemat.cpp:2343
virtual ~DenseMatrixInverse()
Destroys dense inverse matrix.
Definition: densemat.cpp:3995
void UMult(int m, int n, double *X) const
Definition: densemat.cpp:3685
void GetDiag(Vector &d) const
Returns the diagonal of the matrix.
Definition: densemat.cpp:1345
DenseTensor(int i, int j, int k, MemoryType mt)
Definition: densemat.hpp:1005
DenseTensor(double *d, int i, int j, int k)
Definition: densemat.hpp:998
DenseMatrixGeneralizedEigensystem(DenseMatrix &a, DenseMatrix &b, bool left_eigen_vectors=false, bool right_eigen_vectors=false)
Definition: densemat.cpp:4068
void LeftScaling(const Vector &s)
LeftScaling this = diag(s) * this.
Definition: densemat.cpp:354
void AddMult_a(double a, const Vector &x, Vector &y) const
y += a * A.x
Definition: densemat.cpp:299
void AddMultVVt(const Vector &v, DenseMatrix &VVt)
VVt += v v^t.
Definition: densemat.cpp:3187
void CopyMNDiag(double c, int n, int row_offset, int col_offset)
Copy c on the diagonal of size n to *this at row_offset, col_offset.
Definition: densemat.cpp:1662
double * GetData(int k)
Definition: densemat.hpp:1083
void SingularValues(Vector &sv) const
Definition: densemat.cpp:1213
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
Definition: densemat.cpp:3233
void Neg()
(*this) = -(*this)
Definition: densemat.cpp:670
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1331
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition: densemat.hpp:1030
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: densemat.cpp:3956
double FNorm2() const
Compute the square of the Frobenius norm of the matrix.
Definition: densemat.hpp:269
int Capacity() const
Return the size of the allocated memory.
double * data
Definition: densemat.hpp:622
void LSolve(int m, int n, double *X) const
Definition: densemat.cpp:3379
Definition: densemat.hpp:1128
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3252
void SetData(double *d)
Definition: vector.hpp:149
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Definition: device.hpp:319
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
DenseMatrix & Eigenvectors()
Definition: densemat.hpp:898
void Reset()
Reset the memory to be empty, ensuring that Delete() will be a no-op.
void LMult(int m, int n, double *X) const
Definition: densemat.cpp:3666
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
void Threshold(double eps)
Replace small entries, abs(a_ij) <= eps, with zero.
Definition: densemat.cpp:2208
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2587
void TestInversion()
Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
Definition: densemat.cpp:3984
int TotalSize() const
Definition: densemat.hpp:1028
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:635
void CalcEigenvalues(double *lambda, double *vec) const
Definition: densemat.cpp:1293
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1420
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
y += a * A.x
Definition: densemat.cpp:252
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
Definition: densemat.cpp:3136
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2886
void BatchLUSolve(const DenseTensor &Mlu, const Array< int > &P, Vector &X)
Solve batch linear systems.
Definition: densemat.cpp:4366
int Size() const
Get the size of the inverse matrix.
Definition: densemat.hpp:843
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1315
bool OwnsData() const
Return the DenseMatrix data (host pointer) ownership flag.
Definition: densemat.hpp:121
virtual bool Factor(int m, double TOL=0.0)
Compute the Cholesky factorization of the current matrix.
Definition: densemat.cpp:3614
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
double InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Definition: densemat.hpp:198
DenseTensor(int i, int j, int k)
Definition: densemat.hpp:991
void ClearExternalData()
Definition: densemat.hpp:95
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2719
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:312
void USolve(int m, int n, double *X) const
Definition: densemat.cpp:3733
void Swap(DenseMatrix &other)
Definition: densemat.cpp:2307
void AddMatrix(DenseMatrix &A, int ro, int co)
Perform (ro+i,co+j)+=A(i,j) for 0<=i<A.Height, 0<=j<A.Width.
Definition: densemat.cpp:1723
DenseMatrix(double *d, int h, int w)
Construct a DenseMatrix using an existing data array.
Definition: densemat.hpp:58
double * HostWrite()
Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:1124
void Clear()
Delete the matrix data array (if owned) and reset the matrix state.
Definition: densemat.hpp:98
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2654
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
Definition: densemat.hpp:1132
A class to initialize the size of a Tensor.
Definition: dtensor.hpp:54
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:156
void SetSubMatrix(const Array< int > &idx, const DenseMatrix &A)
Set (*this)(idx[i],idx[j]) = A(i,j)
Definition: densemat.cpp:1887
ADAt = A D A^t, where D is diagonal.
Definition: densemat.cpp:2765
void Eigenvalues(DenseMatrix &b, Vector &ev)
Definition: densemat.hpp:285
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
Definition: densemat.hpp:482
void CopyCols(const DenseMatrix &A, int col1, int col2)
Copy columns col1 through col2 from A to *this.
Definition: densemat.cpp:1579
double a
Definition: lissajous.cpp:41
double * GetColumn(int col)
Definition: densemat.hpp:309
ADBt = A D B^t, where D is diagonal.
Definition: densemat.cpp:2943
virtual ~DenseMatrix()
Destroys dense matrix.
Definition: densemat.cpp:2314
T * ReadWrite(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read+write access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true, or the mfem::Device&#39;s HostMemoryClass, otherwise.
Definition: device.hpp:353
void CopyMNt(const DenseMatrix &A, int row_offset, int col_offset)
Copy matrix A^t to the location in *this at row_offset, col_offset.
Definition: densemat.cpp:1618
void CopyExceptMN(const DenseMatrix &A, int m, int n)
Copy All rows and columns except m and n from A.
Definition: densemat.cpp:1697
void GetSubMatrix(const Array< int > &idx, DenseMatrix &A) const
Definition: densemat.cpp:1783
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
void Diag(double c, int n)
Creates n x n diagonal matrix with diagonal elements c.
Definition: densemat.cpp:1390
std::size_t MemoryUsage() const
Definition: densemat.hpp:1109
virtual bool Factor(int m, double TOL=0.0)
Compute the LU factorization of the current matrix.
Definition: densemat.cpp:3269
DenseMatrixInverse(bool spd_=false)
Default constructor.
Definition: densemat.hpp:833
const double * GetData(int k) const
Definition: densemat.hpp:1089
static const int ipiv_base
Definition: densemat.hpp:661
Definition: densemat.cpp:1478
int Rank(double tol) const
Definition: densemat.cpp:1253
virtual void GetInverseMatrix(int m, double *X) const
Assuming L.L^t = A factored data of size (m x m), compute X <- A^{-1}.
Definition: densemat.cpp:3824
int SizeI() const
Definition: densemat.hpp:1024
std::size_t MemoryUsage() const
Definition: densemat.hpp:463
virtual bool Factor(int m, double TOL=0.0)
Definition: densemat.hpp:628
virtual MatrixInverse * Inverse() const
Returns a pointer to the inverse matrix.
Definition: densemat.cpp:483
bool OwnsHostPtr() const
Return true if the host pointer is owned. Ownership indicates whether the pointer will be deleted by ...
void Norm2(Vector &v) const
Take the 2-norm of the columns of A and store in v.
Definition: densemat.hpp:256
DenseMatrix & RightSingularvectors()
Definition: densemat.hpp:968
RefCoord t[3]
void Eigenvalues(Vector &ev, DenseMatrix &evect)
Compute eigenvalues and eigenvectors of A x = ev x where A = *this.
Definition: densemat.hpp:276
DenseMatrixEigensystem(DenseMatrix &m)
Definition: densemat.cpp:4010
const double alpha
Definition: ex15.cpp:369
double operator*(const DenseMatrix &m) const
Matrix inner product: tr(A^t B)
Definition: densemat.cpp:200
LUFactors(double *data_, int *ipiv_)
Definition: densemat.hpp:670
void RightScaling(const Vector &s)
RightScaling: this = this * diag(s);.
Definition: densemat.cpp:380
DenseTensor(const DenseTensor &other)
Copy constructor: deep copy.
Definition: densemat.hpp:1013
Vector data type.
Definition: vector.hpp:60
void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const
Definition: densemat.cpp:3569
void BlockBackSolve(int m, int n, int r, const double *U12, const double *X2, double *Y1) const
Definition: densemat.cpp:3604
ADAt += A D A^t, where D is diagonal.
Definition: densemat.cpp:2737
void GetFromVector(int offset, const Vector &v)
Get the matrix &#39;data&#39; from the Vector &#39;v&#39; at the given &#39;offset&#39;.
Definition: densemat.cpp:2124
void AddMult(const Table &elem_dof, const Vector &x, Vector &y) const
Definition: densemat.cpp:4227
Memory< double > & GetMemory()
Definition: densemat.hpp:1099
int SizeK() const
Definition: densemat.hpp:1026
void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco)
Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this.
Definition: densemat.cpp:1592
double * Write(bool on_dev=true)
Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:1120
const double & operator()(int i, int j, int k) const
Definition: densemat.hpp:1075
RefCoord s[3]
double * Write(bool on_dev=true)
Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:474
double * Data()
Definition: densemat.hpp:1095
void InvLeftScaling(const Vector &s)
InvLeftScaling this = diag(1./s) * this.
Definition: densemat.cpp:367
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:80
void Eigensystem(Vector &ev, DenseMatrix &evect)
Compute eigenvalues and eigenvectors of A x = ev x where A = *this.
Definition: densemat.hpp:280
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
Abstract operator.
Definition: operator.hpp:24
void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += alpha * B * C.
Definition: densemat.cpp:2414
Vector & Singularvalues()
Definition: densemat.hpp:965
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:978
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.cpp:163
Definition: densemat.cpp:2135
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:466
Definition: densemat.cpp:1551
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:3095
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
virtual void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
y += a * A^t x
Definition: densemat.cpp:275
void AddSubMatrix(const Array< int > &idx, const DenseMatrix &A)
(*this)(idx[i],idx[j]) += A(i,j)
Definition: densemat.cpp:2000
MemoryType GetMemoryType() const
Return a MemoryType that is currently valid. If both the host and the device pointers are currently v...
void Eigenvalues(DenseMatrix &b, Vector &ev, DenseMatrix &evect)
Compute generalized eigenvalues of A x = ev B x, where A = *this.
Definition: densemat.hpp:289
void AddToVector(int offset, Vector &v) const
Add the matrix &#39;data&#39; to the Vector &#39;v&#39; at the given &#39;offset&#39;.
Definition: densemat.cpp:2113
double FNorm() const
Compute the Frobenius norm of the matrix.
Definition: densemat.hpp:266
DenseMatrix & operator+=(const double *m)
Definition: densemat.cpp:634
CholeskyFactors(double *data_)
Definition: densemat.hpp:769