MFEM  v4.4.0 Finite element discretization library
densemat.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #ifndef MFEM_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  virtual void Mult(const Vector &x, Vector &y) const;
146
147  /// Multiply a vector with the transpose matrix.
148  void MultTranspose(const double *x, double *y) const;
149
150  /// Multiply a vector with the transpose matrix.
151  virtual void MultTranspose(const Vector &x, Vector &y) const;
152
153  /// y += A.x
154  void AddMult(const Vector &x, Vector &y) const;
155
156  /// y += A^t x
157  void AddMultTranspose(const Vector &x, Vector &y) const;
158
159  /// y += a * A.x
160  void AddMult_a(double a, const Vector &x, Vector &y) const;
161
162  /// y += a * A^t x
163  void AddMultTranspose_a(double a, const Vector &x, Vector &y) const;
164
165  /// Compute y^t A x
166  double InnerProduct(const double *x, const double *y) const;
167
168  /// LeftScaling this = diag(s) * this
169  void LeftScaling(const Vector & s);
170  /// InvLeftScaling this = diag(1./s) * this
171  void InvLeftScaling(const Vector & s);
172  /// RightScaling: this = this * diag(s);
173  void RightScaling(const Vector & s);
174  /// InvRightScaling: this = this * diag(1./s);
175  void InvRightScaling(const Vector & s);
176  /// SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
177  void SymmetricScaling(const Vector & s);
178  /// InvSymmetricScaling this = diag(sqrt(1./s)) * this * diag(sqrt(1./s))
179  void InvSymmetricScaling(const Vector & s);
180
181  /// Compute y^t A x
182  double InnerProduct(const Vector &x, const Vector &y) const
183  { return InnerProduct((const double *)x, (const double *)y); }
184
185  /// Returns a pointer to the inverse matrix.
186  virtual MatrixInverse *Inverse() const;
187
188  /// Replaces the current matrix with its inverse
189  void Invert();
190
191  /// Replaces the current matrix with its square root inverse
192  void SquareRootInverse();
193
194  /// Calculates the determinant of the matrix
195  /// (optimized for 2x2, 3x3, and 4x4 matrices)
196  double Det() const;
197
198  double Weight() const;
199
200  /** @brief Set the matrix to alpha * A, assuming that A has the same
201  dimensions as the matrix and uses column-major layout. */
202  void Set(double alpha, const double *A);
203  /// Set the matrix to alpha * A.
204  void Set(double alpha, const DenseMatrix &A)
205  {
206  SetSize(A.Height(), A.Width());
207  Set(alpha, A.GetData());
208  }
209
210  /// Adds the matrix A multiplied by the number c to the matrix
211  void Add(const double c, const DenseMatrix &A);
212
213  /// Sets the matrix elements equal to constant c
214  DenseMatrix &operator=(double c);
215
216  /// Copy the matrix entries from the given array
217  DenseMatrix &operator=(const double *d);
218
219  /// Sets the matrix size and elements equal to those of m
220  DenseMatrix &operator=(const DenseMatrix &m);
221
222  DenseMatrix &operator+=(const double *m);
224
226
227  DenseMatrix &operator*=(double c);
228
229  /// (*this) = -(*this)
230  void Neg();
231
232  /// Take the 2-norm of the columns of A and store in v
233  void Norm2(double *v) const;
234
235  /// Compute the norm ||A|| = max_{ij} |A_{ij}|
236  double MaxMaxNorm() const;
237
238  /// Compute the Frobenius norm of the matrix
239  double FNorm() const { double s, n2; FNorm(s, n2); return s*sqrt(n2); }
240
241  /// Compute the square of the Frobenius norm of the matrix
242  double FNorm2() const { double s, n2; FNorm(s, n2); return s*s*n2; }
243
244  /// Compute eigenvalues of A x = ev x where A = *this
245  void Eigenvalues(Vector &ev)
246  { Eigensystem(ev); }
247
248  /// Compute eigenvalues and eigenvectors of A x = ev x where A = *this
249  void Eigenvalues(Vector &ev, DenseMatrix &evect)
250  { Eigensystem(ev, &evect); }
251
252  /// Compute eigenvalues and eigenvectors of A x = ev x where A = *this
253  void Eigensystem(Vector &ev, DenseMatrix &evect)
254  { Eigensystem(ev, &evect); }
255
256  /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
257  where A = *this */
259  { Eigensystem(b, ev); }
260
261  /// Compute generalized eigenvalues of A x = ev B x, where A = *this
263  { Eigensystem(b, ev, &evect); }
264
265  /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
266  where A = *this */
268  { Eigensystem(b, ev, &evect); }
269
270  void SingularValues(Vector &sv) const;
271  int Rank(double tol) const;
272
273  /// Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
274  double CalcSingularvalue(const int i) const;
275
276  /** Return the eigenvalues (in increasing order) and eigenvectors of a
277  2x2 or 3x3 symmetric matrix. */
278  void CalcEigenvalues(double *lambda, double *vec) const;
279
280  void GetRow(int r, Vector &row) const;
281  void GetColumn(int c, Vector &col) const;
282  double *GetColumn(int col) { return data + col*height; }
283  const double *GetColumn(int col) const { return data + col*height; }
284
285  void GetColumnReference(int c, Vector &col)
286  { col.SetDataAndSize(data + c * height, height); }
287
288  void SetRow(int r, const double* row);
289  void SetRow(int r, const Vector &row);
290
291  void SetCol(int c, const double* col);
292  void SetCol(int c, const Vector &col);
293
294
295  /// Set all entries of a row to the specified value.
296  void SetRow(int row, double value);
297  /// Set all entries of a column to the specified value.
298  void SetCol(int col, double value);
299
300  /// Returns the diagonal of the matrix
301  void GetDiag(Vector &d) const;
302  /// Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|
303  void Getl1Diag(Vector &l) const;
304  /// Compute the row sums of the DenseMatrix
305  void GetRowSums(Vector &l) const;
306
307  /// Creates n x n diagonal matrix with diagonal elements c
308  void Diag(double c, int n);
309  /// Creates n x n diagonal matrix with diagonal given by diag
310  void Diag(double *diag, int n);
311
312  /// (*this) = (*this)^t
313  void Transpose();
314  /// (*this) = A^t
315  void Transpose(const DenseMatrix &A);
316  /// (*this) = 1/2 ((*this) + (*this)^t)
317  void Symmetrize();
318
319  void Lump();
320
321  /** Given a DShape matrix (from a scalar FE), stored in *this, returns the
322  CurlShape matrix. If *this is a N by D matrix, then curl is a D*N by
323  D*(D-1)/2 matrix. The size of curl must be set outside. The dimension D
324  can be either 2 or 3. */
325  void GradToCurl(DenseMatrix &curl);
326  /** Given a DShape matrix (from a scalar FE), stored in *this,
327  returns the DivShape vector. If *this is a N by dim matrix,
328  then div is a dim*N vector. The size of div must be set
329  outside. */
330  void GradToDiv(Vector &div);
331
332  /// Copy rows row1 through row2 from A to *this
333  void CopyRows(const DenseMatrix &A, int row1, int row2);
334  /// Copy columns col1 through col2 from A to *this
335  void CopyCols(const DenseMatrix &A, int col1, int col2);
336  /// Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this
337  void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco);
338  /// Copy matrix A to the location in *this at row_offset, col_offset
339  void CopyMN(const DenseMatrix &A, int row_offset, int col_offset);
340  /// Copy matrix A^t to the location in *this at row_offset, col_offset
341  void CopyMNt(const DenseMatrix &A, int row_offset, int col_offset);
342  /** Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this at
343  row_offset, col_offset */
344  void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco,
345  int row_offset, int col_offset);
346  /// Copy c on the diagonal of size n to *this at row_offset, col_offset
347  void CopyMNDiag(double c, int n, int row_offset, int col_offset);
348  /// Copy diag on the diagonal of size n to *this at row_offset, col_offset
349  void CopyMNDiag(double *diag, int n, int row_offset, int col_offset);
350  /// Copy All rows and columns except m and n from A
351  void CopyExceptMN(const DenseMatrix &A, int m, int n);
352
353  /// Perform (ro+i,co+j)+=A(i,j) for 0<=i<A.Height, 0<=j<A.Width
354  void AddMatrix(DenseMatrix &A, int ro, int co);
355  /// Perform (ro+i,co+j)+=a*A(i,j) for 0<=i<A.Height, 0<=j<A.Width
356  void AddMatrix(double a, const DenseMatrix &A, int ro, int co);
357
358  /// Add the matrix 'data' to the Vector 'v' at the given 'offset'
359  void AddToVector(int offset, Vector &v) const;
360  /// Get the matrix 'data' from the Vector 'v' at the given 'offset'
361  void GetFromVector(int offset, const Vector &v);
362  /** If (dofs[i] < 0 and dofs[j] >= 0) or (dofs[i] >= 0 and dofs[j] < 0)
363  then (*this)(i,j) = -(*this)(i,j). */
364  void AdjustDofDirection(Array<int> &dofs);
365
366  /// Replace small entries, abs(a_ij) <= eps, with zero.
367  void Threshold(double eps);
368
369  /** Count the number of entries in the matrix for which isfinite
370  is false, i.e. the entry is a NaN or +/-Inf. */
371  int CheckFinite() const { return mfem::CheckFinite(data, height*width); }
372
373  /// Prints matrix to stream out.
374  virtual void Print(std::ostream &out = mfem::out, int width_ = 4) const;
375  virtual void PrintMatlab(std::ostream &out = mfem::out) const;
376  /// Prints the transpose matrix to stream out.
377  virtual void PrintT(std::ostream &out = mfem::out, int width_ = 4) const;
378
379  /// Invert and print the numerical conditioning of the inversion.
380  void TestInversion();
381
382  long MemoryUsage() const { return data.Capacity() * sizeof(double); }
383
384  /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
385  const double *Read(bool on_dev = true) const
386  { return mfem::Read(data, Height()*Width(), on_dev); }
387
388  /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
389  const double *HostRead() const
390  { return mfem::Read(data, Height()*Width(), false); }
391
392  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
393  double *Write(bool on_dev = true)
394  { return mfem::Write(data, Height()*Width(), on_dev); }
395
396  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
397  double *HostWrite()
398  { return mfem::Write(data, Height()*Width(), false); }
399
400  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
401  double *ReadWrite(bool on_dev = true)
402  { return mfem::ReadWrite(data, Height()*Width(), on_dev); }
403
404  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
405  double *HostReadWrite()
406  { return mfem::ReadWrite(data, Height()*Width(), false); }
407
408  void Swap(DenseMatrix &other);
409
410  /// Destroys dense matrix.
411  virtual ~DenseMatrix();
412 };
413
414 /// C = A + alpha*B
415 void Add(const DenseMatrix &A, const DenseMatrix &B,
416  double alpha, DenseMatrix &C);
417
418 /// C = alpha*A + beta*B
419 void Add(double alpha, const double *A,
420  double beta, const double *B, DenseMatrix &C);
421
422 /// C = alpha*A + beta*B
423 void Add(double alpha, const DenseMatrix &A,
424  double beta, const DenseMatrix &B, DenseMatrix &C);
425
426 /// @brief Solves the dense linear system, A * X = B for X
427 ///
428 /// @param [in,out] A the square matrix for the linear system
429 /// @param [in,out] X the rhs vector, B, on input, the solution, X, on output.
430 /// @param [in] TOL optional fuzzy comparison tolerance. Defaults to 1e-9.
431 ///
432 /// @return status set to true if successful, otherwise, false.
433 ///
434 /// @note This routine may replace the contents of the input Matrix, A, with the
435 /// corresponding LU factorization of the matrix. Matrices of size 1x1 and
436 /// 2x2 are handled explicitly.
437 ///
438 /// @pre A.IsSquare() == true
439 /// @pre X != nullptr
440 bool LinearSolve(DenseMatrix& A, double* X, double TOL = 1.e-9);
441
442 /// Matrix matrix multiplication. A = B * C.
443 void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);
444
445 /// Matrix matrix multiplication. A += B * C.
446 void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);
447
448 /// Matrix matrix multiplication. A += alpha * B * C.
449 void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c,
450  DenseMatrix &a);
451
452 /** Calculate the adjugate of a matrix (for NxN matrices, N=1,2,3) or the matrix
453  adj(A^t.A).A^t for rectangular matrices (2x1, 3x1, or 3x2). This operation
454  is well defined even when the matrix is not full rank. */
455 void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja);
456
457 /// Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
458 void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat);
459
460 /** Calculate the inverse of a matrix (for NxN matrices, N=1,2,3) or the
461  left inverse (A^t.A)^{-1}.A^t (for 2x1, 3x1, or 3x2 matrices) */
462 void CalcInverse(const DenseMatrix &a, DenseMatrix &inva);
463
464 /// Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
465 void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva);
466
467 /** For a given Nx(N-1) (N=2,3) matrix J, compute a vector n such that
468  n_k = (-1)^{k+1} det(J_k), k=1,..,N, where J_k is the matrix J with the
469  k-th row removed. Note: J^t.n = 0, det([n|J])=|n|^2=det(J^t.J). */
470 void CalcOrtho(const DenseMatrix &J, Vector &n);
471
472 /// Calculate the matrix A.At
473 void MultAAt(const DenseMatrix &a, DenseMatrix &aat);
474
475 /// ADAt = A D A^t, where D is diagonal
476 void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt);
477
478 /// ADAt += A D A^t, where D is diagonal
479 void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt);
480
481 /// Multiply a matrix A with the transpose of a matrix B: A*Bt
482 void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);
483
484 /// ADBt = A D B^t, where D is diagonal
485 void MultADBt(const DenseMatrix &A, const Vector &D,
486  const DenseMatrix &B, DenseMatrix &ADBt);
487
488 /// ABt += A * B^t
489 void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);
490
491 /// ADBt = A D B^t, where D is diagonal
492 void AddMultADBt(const DenseMatrix &A, const Vector &D,
493  const DenseMatrix &B, DenseMatrix &ADBt);
494
495 /// ABt += a * A * B^t
496 void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B,
497  DenseMatrix &ABt);
498
499 /// Multiply the transpose of a matrix A with a matrix B: At*B
500 void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB);
501
502 /// AAt += a * A * A^t
503 void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);
504
505 /// AAt = a * A * A^t
506 void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);
507
508 /// Make a matrix from a vector V.Vt
509 void MultVVt(const Vector &v, DenseMatrix &vvt);
510
511 void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);
512
513 /// VWt += v w^t
514 void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);
515
516 /// VVt += v v^t
517 void AddMultVVt(const Vector &v, DenseMatrix &VWt);
518
519 /// VWt += a * v w^t
520 void AddMult_a_VWt(const double a, const Vector &v, const Vector &w,
521  DenseMatrix &VWt);
522
523 /// VVt += a * v v^t
524 void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt);
525
526
527 /** Class that can compute LU factorization of external data and perform various
528  operations with the factored data. */
530 {
531 public:
532  double *data;
533  int *ipiv;
534 #ifdef MFEM_USE_LAPACK
535  static const int ipiv_base = 1;
536 #else
537  static const int ipiv_base = 0;
538 #endif
539
540  /** With this constructor, the (public) data and ipiv members should be set
541  explicitly before calling class methods. */
542  LUFactors() { }
543
544  LUFactors(double *data_, int *ipiv_) : data(data_), ipiv(ipiv_) { }
545
546  /**
547  * @brief Compute the LU factorization of the current matrix
548  *
549  * Factorize the current matrix of size (m x m) overwriting it with the
550  * LU factors. The factorization is such that L.U = P.A, where A is the
551  * original matrix and P is a permutation matrix represented by ipiv.
552  *
553  * @param [in] m size of the square matrix
554  * @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0.
555  *
556  * @return status set to true if successful, otherwise, false.
557  */
558  bool Factor(int m, double TOL = 0.0);
559
560  /** Assuming L.U = P.A factored data of size (m x m), compute |A|
561  from the diagonal values of U and the permutation information. */
562  double Det(int m) const;
563
564  /** Assuming L.U = P.A factored data of size (m x m), compute X <- A X,
565  for a matrix X of size (m x n). */
566  void Mult(int m, int n, double *X) const;
567
568  /** Assuming L.U = P.A factored data of size (m x m), compute
569  X <- L^{-1} P X, for a matrix X of size (m x n). */
570  void LSolve(int m, int n, double *X) const;
571
572  /** Assuming L.U = P.A factored data of size (m x m), compute
573  X <- U^{-1} X, for a matrix X of size (m x n). */
574  void USolve(int m, int n, double *X) const;
575
576  /** Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1} X,
577  for a matrix X of size (m x n). */
578  void Solve(int m, int n, double *X) const;
579
580  /** Assuming L.U = P.A factored data of size (m x m), compute X <- X A^{-1},
581  for a matrix X of size (n x m). */
582  void RightSolve(int m, int n, double *X) const;
583
584  /// Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1}.
585  void GetInverseMatrix(int m, double *X) const;
586
587  /** Given an (n x m) matrix A21, compute X2 <- X2 - A21 X1, for matrices X1,
588  and X2 of size (m x r) and (n x r), respectively. */
589  static void SubMult(int m, int n, int r, const double *A21,
590  const double *X1, double *X2);
591
592  /** Assuming P.A = L.U factored data of size (m x m), compute the 2x2 block
593  decomposition:
594  | P 0 | | A A12 | = | L 0 | | U U12 |
595  | 0 I | | A21 A22 | | L21 I | | 0 S22 |
596  where A12, A21, and A22 are matrices of size (m x n), (n x m), and
597  (n x n), respectively. The blocks are overwritten as follows:
598  A12 <- U12 = L^{-1} P A12
599  A21 <- L21 = A21 U^{-1}
600  A22 <- S22 = A22 - L21 U12.
601  The block S22 is the Schur complement. */
602  void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const;
603
604  /** Given BlockFactor()'d data, perform the forward block solve for the
605  linear system:
606  | A A12 | | X1 | = | B1 |
607  | A21 A22 | | X2 | | B2 |
608  written in the factored form:
609  | L 0 | | U U12 | | X1 | = | P 0 | | B1 |
610  | L21 I | | 0 S22 | | X2 | | 0 I | | B2 |.
611  The resulting blocks Y1, Y2 solve the system:
612  | L 0 | | Y1 | = | P 0 | | B1 |
613  | L21 I | | Y2 | | 0 I | | B2 |
614  The blocks are overwritten as follows:
615  B1 <- Y1 = L^{-1} P B1
616  B2 <- Y2 = B2 - L21 Y1 = B2 - A21 A^{-1} B1
617  The blocks B1/Y1 and B2/Y2 are of size (m x r) and (n x r), respectively.
618  The Schur complement system is given by: S22 X2 = Y2. */
619  void BlockForwSolve(int m, int n, int r, const double *L21,
620  double *B1, double *B2) const;
621
622  /** Given BlockFactor()'d data, perform the backward block solve in
623  | U U12 | | X1 | = | Y1 |
624  | 0 S22 | | X2 | | Y2 |.
625  The input is the solution block X2 and the block Y1 resulting from
626  BlockForwSolve(). The result block X1 overwrites input block Y1:
627  Y1 <- X1 = U^{-1} (Y1 - U12 X2). */
628  void BlockBackSolve(int m, int n, int r, const double *U12,
629  const double *X2, double *Y1) const;
630 };
631
632
633 /** Data type for inverse of square dense matrix.
634  Stores LU factors */
636 {
637 private:
638  const DenseMatrix *a;
639  LUFactors lu;
640
641 public:
642  /// Default constructor.
643  DenseMatrixInverse() : a(NULL), lu(NULL, NULL) { }
644
645  /** Creates square dense matrix. Computes factorization of mat
646  and stores LU factors. */
647  DenseMatrixInverse(const DenseMatrix &mat);
648
649  /// Same as above but does not factorize the matrix.
650  DenseMatrixInverse(const DenseMatrix *mat);
651
652  /// Get the size of the inverse matrix
653  int Size() const { return Width(); }
654
655  /// Factor the current DenseMatrix, *a
656  void Factor();
657
658  /// Factor a new DenseMatrix of the same size
659  void Factor(const DenseMatrix &mat);
660
661  virtual void SetOperator(const Operator &op);
662
663  /// Matrix vector multiplication with the inverse of dense matrix.
664  void Mult(const double *x, double *y) const;
665
666  /// Matrix vector multiplication with the inverse of dense matrix.
667  virtual void Mult(const Vector &x, Vector &y) const;
668
669  /// Multiply the inverse matrix by another matrix: X = A^{-1} B.
670  void Mult(const DenseMatrix &B, DenseMatrix &X) const;
671
672  /// Multiply the inverse matrix by another matrix: X <- A^{-1} X.
673  void Mult(DenseMatrix &X) const { lu.Solve(width, X.Width(), X.Data()); }
674
675  /// Compute and return the inverse matrix in Ainv.
676  void GetInverseMatrix(DenseMatrix &Ainv) const;
677
678  /// Compute the determinant of the original DenseMatrix using the LU factors.
679  double Det() const { return lu.Det(width); }
680
681  /// Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
682  void TestInversion();
683
684  /// Destroys dense inverse matrix.
685  virtual ~DenseMatrixInverse();
686 };
687
688 #ifdef MFEM_USE_LAPACK
689
691 {
692  DenseMatrix &mat;
693  Vector EVal;
694  DenseMatrix EVect;
695  Vector ev;
696  int n;
697  double *work;
698  char jobz, uplo;
699  int lwork, info;
700 public:
701
704  void Eval();
705  Vector &Eigenvalues() { return EVal; }
706  DenseMatrix &Eigenvectors() { return EVect; }
707  double Eigenvalue(int i) { return EVal(i); }
708  const Vector &Eigenvector(int i)
709  {
710  ev.SetData(EVect.Data() + i * EVect.Height());
711  return ev;
712  }
714 };
715
717 {
718  DenseMatrix &A;
719  DenseMatrix &B;
720  DenseMatrix A_copy;
721  DenseMatrix B_copy;
722  Vector evalues_r;
723  Vector evalues_i;
724  DenseMatrix Vr;
725  DenseMatrix Vl;
726  int n;
727
728  double *alphar;
729  double *alphai;
730  double *beta;
731  double *work;
732  char jobvl, jobvr;
733  int lwork, info;
734
735 public:
736
738  bool left_eigen_vectors = false,
739  bool right_eigen_vectors = false);
740  void Eval();
741  Vector &EigenvaluesRealPart() { return evalues_r; }
742  Vector &EigenvaluesImagPart() { return evalues_i; }
743  double EigenvalueRealPart(int i) { return evalues_r(i); }
744  double EigenvalueImagPart(int i) { return evalues_i(i); }
745  DenseMatrix &LeftEigenvectors() { return Vl; }
746  DenseMatrix &RightEigenvectors() { return Vr; }
748 };
749
750
752 {
753  DenseMatrix Mc;
754  Vector sv;
755  DenseMatrix U,Vt;
756  int m, n;
757
758 #ifdef MFEM_USE_LAPACK
759  double *work;
760  char jobu, jobvt;
761  int lwork, info;
762 #endif
763
764  void Init();
765 public:
767  bool left_singular_vectors=false,
768  bool right_singlular_vectors=false);
769  DenseMatrixSVD(int h, int w,
770  bool left_singular_vectors=false,
771  bool right_singlular_vectors=false);
772  void Eval(DenseMatrix &M);
773  Vector &Singularvalues() { return sv; }
774  double Singularvalue(int i) { return sv(i); }
777  ~DenseMatrixSVD();
778 };
779
780 #endif // if MFEM_USE_LAPACK
781
782
783 class Table;
784
785 /// Rank 3 tensor (array of matrices)
787 {
788 private:
789  DenseMatrix Mk;
790  Memory<double> tdata;
791  int nk;
792
793 public:
795  {
796  nk = 0;
797  }
798
799  DenseTensor(int i, int j, int k)
800  : Mk(NULL, i, j)
801  {
802  nk = k;
803  tdata.New(i*j*k);
804  }
805
806  DenseTensor(double *d, int i, int j, int k)
807  : Mk(NULL, i, j)
808  {
809  nk = k;
810  tdata.Wrap(d, i*j*k, false);
811  }
812
813  DenseTensor(int i, int j, int k, MemoryType mt)
814  : Mk(NULL, i, j)
815  {
816  nk = k;
817  tdata.New(i*j*k, mt);
818  }
819
820  /// Copy constructor: deep copy
821  DenseTensor(const DenseTensor &other)
822  : Mk(NULL, other.Mk.height, other.Mk.width), nk(other.nk)
823  {
824  const int size = Mk.Height()*Mk.Width()*nk;
825  if (size > 0)
826  {
827  tdata.New(size, other.tdata.GetMemoryType());
828  tdata.CopyFrom(other.tdata, size);
829  }
830  }
831
832  int SizeI() const { return Mk.Height(); }
833  int SizeJ() const { return Mk.Width(); }
834  int SizeK() const { return nk; }
835
836  int TotalSize() const { return SizeI()*SizeJ()*SizeK(); }
837
838  void SetSize(int i, int j, int k, MemoryType mt_ = MemoryType::PRESERVE)
839  {
840  const MemoryType mt = mt_ == MemoryType::PRESERVE ? tdata.GetMemoryType() : mt_;
841  tdata.Delete();
842  Mk.UseExternalData(NULL, i, j);
843  nk = k;
844  tdata.New(i*j*k, mt);
845  }
846
847  void UseExternalData(double *ext_data, int i, int j, int k)
848  {
849  tdata.Delete();
850  Mk.UseExternalData(NULL, i, j);
851  nk = k;
852  tdata.Wrap(ext_data, i*j*k, false);
853  }
854
855  /// Sets the tensor elements equal to constant c
856  DenseTensor &operator=(double c);
857
858  /// Copy assignment operator (performs a deep copy)
859  DenseTensor &operator=(const DenseTensor &other);
860
862  {
863  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
864  Mk.data = Memory<double>(GetData(k), SizeI()*SizeJ(), false);
865  return Mk;
866  }
867  const DenseMatrix &operator()(int k) const
868  { return const_cast<DenseTensor&>(*this)(k); }
869
870  double &operator()(int i, int j, int k)
871  {
872  MFEM_ASSERT_INDEX_IN_RANGE(i, 0, SizeI());
873  MFEM_ASSERT_INDEX_IN_RANGE(j, 0, SizeJ());
874  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
875  return tdata[i+SizeI()*(j+SizeJ()*k)];
876  }
877
878  const double &operator()(int i, int j, int k) const
879  {
880  MFEM_ASSERT_INDEX_IN_RANGE(i, 0, SizeI());
881  MFEM_ASSERT_INDEX_IN_RANGE(j, 0, SizeJ());
882  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
883  return tdata[i+SizeI()*(j+SizeJ()*k)];
884  }
885
886  double *GetData(int k)
887  {
888  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
889  return tdata+k*Mk.Height()*Mk.Width();
890  }
891
892  double *Data() { return tdata; }
893
894  const double *Data() const { return tdata; }
895
896  Memory<double> &GetMemory() { return tdata; }
897  const Memory<double> &GetMemory() const { return tdata; }
898
899  /** Matrix-vector product from unassembled element matrices, assuming both
900  'x' and 'y' use the same elem_dof table. */
901  void AddMult(const Table &elem_dof, const Vector &x, Vector &y) const;
902
903  void Clear()
904  { UseExternalData(NULL, 0, 0, 0); }
905
906  long MemoryUsage() const { return nk*Mk.MemoryUsage(); }
907
908  /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
909  const double *Read(bool on_dev = true) const
910  { return mfem::Read(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
911
912  /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
913  const double *HostRead() const
914  { return mfem::Read(tdata, Mk.Height()*Mk.Width()*nk, false); }
915
916  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
917  double *Write(bool on_dev = true)
918  { return mfem::Write(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
919
920  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
921  double *HostWrite()
922  { return mfem::Write(tdata, Mk.Height()*Mk.Width()*nk, false); }
923
924  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
925  double *ReadWrite(bool on_dev = true)
926  { return mfem::ReadWrite(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
927
928  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
929  double *HostReadWrite()
930  { return mfem::ReadWrite(tdata, Mk.Height()*Mk.Width()*nk, false); }
931
933  {
934  mfem::Swap(tdata, t.tdata);
935  mfem::Swap(nk, t.nk);
936  Mk.Swap(t.Mk);
937  }
938
939  ~DenseTensor() { tdata.Delete(); }
940 };
941
942 /** @brief Compute the LU factorization of a batch of matrices
943
944  Factorize n matrices of size (m x m) stored in a dense tensor overwriting it
945  with the LU factors. The factorization is such that L.U = Piv.A, where A is
946  the original matrix and Piv is a permutation matrix represented by P.
947
948  @param [in, out] Mlu batch of square matrices - dimension m x m x n.
949  @param [out] P array storing pivot information - dimension m x n.
950  @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0. */
951 void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const double TOL = 0.0);
952
953 /** @brief Solve batch linear systems
954
955  Assuming L.U = P.A for n factored matrices (m x m), compute x <- A x, for n
956  companion vectors.
957
958  @param [in] Mlu batch of LU factors for matrix M - dimension m x m x n.
959  @param [in] P array storing pivot information - dimension m x n.
960  @param [in, out] X vector storing right-hand side and then solution -
961  dimension m x n. */
962 void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X);
963
964
965 // Inline methods
966
967 inline double &DenseMatrix::operator()(int i, int j)
968 {
969  MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
970  return data[i+j*height];
971 }
972
973 inline const double &DenseMatrix::operator()(int i, int j) const
974 {
975  MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
976  return data[i+j*height];
977 }
978
979 } // namespace mfem
980
981 #endif
int Size() const
For backward compatibility define Size to be synonym of Width()
Definition: densemat.hpp:102
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
Definition: densemat.cpp:1392
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:2376
DenseMatrix & operator-=(const DenseMatrix &m)
Definition: densemat.cpp:586
void SymmetricScaling(const Vector &s)
SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
Definition: densemat.cpp:357
void SquareRootInverse()
Replaces the current matrix with its square root inverse.
Definition: densemat.cpp:728
int CheckFinite(const double *v, const int n)
Definition: vector.hpp:501
void BatchLUFactor(DenseTensor &Mlu, Array< int > &P, const double TOL)
Compute the LU factorization of a batch of matrices.
Definition: densemat.cpp:3582
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:2761
void Mult(DenseMatrix &X) const
Multiply the inverse matrix by another matrix: X &lt;- A^{-1} X.
Definition: densemat.hpp:673
DenseMatrix & operator*=(double c)
Definition: densemat.cpp:599
void GetDiag(Vector &d) const
Returns the diagonal of the matrix.
Definition: densemat.cpp:1284
void SetCol(int c, const double *col)
Definition: densemat.cpp:1788
void UseExternalData(double *ext_data, int i, int j, int k)
Definition: densemat.hpp:847
DenseTensor & operator=(double c)
Sets the tensor elements equal to constant c.
Definition: densemat.cpp:3565
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:2742
double * HostWrite()
Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:397
Memory< double > & GetMemory()
Definition: densemat.hpp:117
void SetRow(int r, const double *row)
Definition: densemat.cpp:1773
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
Definition: densemat.cpp:343
const DenseMatrix & operator()(int k) const
Definition: densemat.hpp:867
void Eigenvalues(Vector &ev)
Compute eigenvalues of A x = ev x where A = *this.
Definition: densemat.hpp:245
void SingularValues(Vector &sv) const
Definition: densemat.cpp:1152
void Delete()
Delete the owned pointers and reset the Memory object.
DenseMatrix & operator()(int k)
Definition: densemat.hpp:861
double Det() const
Definition: densemat.cpp:436
void Eigensystem(DenseMatrix &b, Vector &ev, DenseMatrix &evect)
Definition: densemat.hpp:267
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:472
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
int SizeK() const
Definition: densemat.hpp:834
long MemoryUsage() const
Definition: densemat.hpp:906
void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const
Definition: densemat.cpp:3148
void BlockBackSolve(int m, int n, int r, const double *U12, const double *X2, double *Y1) const
Definition: densemat.cpp:3183
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
Definition: densemat.cpp:284
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:2074
DenseMatrixSVD(DenseMatrix &M, bool left_singular_vectors=false, bool right_singlular_vectors=false)
Definition: densemat.cpp:3437
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:909
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 AddMult(const Table &elem_dof, const Vector &x, Vector &y) const
Definition: densemat.cpp:3510
void TestInversion()
Invert and print the numerical conditioning of the inversion.
Definition: densemat.cpp:1888
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
DenseMatrix & LeftSingularvectors()
Definition: densemat.hpp:775
void CopyRows(const DenseMatrix &A, int row1, int row2)
Copy rows row1 through row2 from A to *this.
Definition: densemat.cpp:1491
double Det(int m) const
Definition: densemat.cpp:2907
void Swap(DenseTensor &t)
Definition: densemat.hpp:932
void Eval(DenseMatrix &M)
Definition: densemat.cpp:3471
const double * Data() const
Definition: densemat.hpp:894
Abstract data type for matrix inverse.
Definition: matrix.hpp:62
void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += a * A * B^t.
Definition: densemat.cpp:2574
double * HostReadWrite()
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:405
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:3224
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3212
void CopyFrom(const Memory &src, int size)
Copy size entries from src to *this.
void GetInverseMatrix(int m, double *X) const
Assuming L.U = P.A factored data of size (m x m), compute X &lt;- A^{-1}.
Definition: densemat.cpp:3068
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
const Memory< double > & GetMemory() const
Definition: densemat.hpp:897
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2288
DenseMatrix & operator=(double c)
Sets the matrix elements equal to constant c.
Definition: densemat.cpp:540
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:520
int Capacity() const
Return the size of the allocated memory.
void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt = a * A * A^t.
Definition: densemat.cpp:2715
MemoryType GetMemoryType() const
Return a MemoryType that is currently valid. If both the host and the device pointers are currently v...
static void SubMult(int m, int n, int r, const double *A21, const double *X1, double *X2)
Definition: densemat.cpp:3131
const Vector & Eigenvector(int i)
Definition: densemat.hpp:708
virtual void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: densemat.cpp:1817
const double & operator()(int i, int j, int k) const
Definition: densemat.hpp:878
void Set(double alpha, const DenseMatrix &A)
Set the matrix to alpha * A.
Definition: densemat.hpp:204
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1916
double & operator()(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.hpp:967
double Weight() const
Definition: densemat.cpp:493
void Wrap(T *ptr, int size, bool own)
Wrap an externally allocated host pointer, ptr with the current host memory type returned by MemoryMa...
void USolve(int m, int n, double *X) const
Definition: densemat.cpp:2981
double FNorm() const
Compute the Frobenius norm of the matrix.
Definition: densemat.hpp:239
const double * HostRead() const
Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:389
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:188
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2146
double & operator()(int i, int j, int k)
Definition: densemat.hpp:870
void RightSolve(int m, int n, double *X) const
Definition: densemat.cpp:3013
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2042
double operator*(const DenseMatrix &m) const
Matrix inner product: tr(A^t B)
Definition: densemat.cpp:173
double Singularvalue(int i)
Definition: densemat.hpp:774
void Reset(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:90
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
DenseMatrix(const double(&values)[M][N])
Definition: densemat.hpp:64
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:529
void AddMult_a_VWt(const double a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
Definition: densemat.cpp:2806
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:385
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1254
void BlockForwSolve(int m, int n, int r, const double *L21, double *B1, double *B2) const
Definition: densemat.cpp:3174
Abstract data type matrix.
Definition: matrix.hpp:27
const double * GetColumn(int col) const
Definition: densemat.hpp:283
void Norm2(double *v) const
Take the 2-norm of the columns of A and store in v.
Definition: densemat.cpp:773
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:2441
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:618
bool LinearSolve(DenseMatrix &A, double *X, double TOL)
Solves the dense linear system, A * X = B for X
Definition: densemat.cpp:1938
virtual ~DenseMatrixInverse()
Destroys dense inverse matrix.
Definition: densemat.cpp:3284
void LSolve(int m, int n, double *X) const
Definition: densemat.cpp:2958
DenseTensor(int i, int j, int k, MemoryType mt)
Definition: densemat.hpp:813
int TotalSize() const
Definition: densemat.hpp:836
DenseTensor(double *d, int i, int j, int k)
Definition: densemat.hpp:806
DenseMatrixGeneralizedEigensystem(DenseMatrix &a, DenseMatrix &b, bool left_eigen_vectors=false, bool right_eigen_vectors=false)
Definition: densemat.cpp:3350
void LeftScaling(const Vector &s)
LeftScaling this = diag(s) * this.
Definition: densemat.cpp:302
double Det() const
Compute the determinant of the original DenseMatrix using the LU factors.
Definition: densemat.hpp:679
void AddMultVVt(const Vector &v, DenseMatrix &VVt)
VVt += v v^t.
Definition: densemat.cpp:2782
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:1587
double * GetData(int k)
Definition: densemat.hpp:886
bool OwnsHostPtr() const
Return true if the host pointer is owned. Ownership indicates whether the pointer will be deleted by ...
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
Definition: densemat.cpp:2828
void Neg()
(*this) = -(*this)
Definition: densemat.cpp:609
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition: densemat.hpp:838
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: densemat.cpp:3245
void Solve(int m, int n, double *X) const
Definition: densemat.cpp:2999
double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:925
void Getl1Diag(Vector &l) const
Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|.
Definition: densemat.cpp:1298
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:1708
void SetData(double *d)
Definition: vector.hpp:149
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true...
Definition: device.hpp:319
long MemoryUsage() const
Definition: densemat.hpp:382
DenseMatrix & Eigenvectors()
Definition: densemat.hpp:706
void Reset()
Reset the memory to be empty, ensuring that Delete() will be a no-op.
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1270
void AddMult(const Vector &x, Vector &y) const
y += A.x
Definition: densemat.cpp:211
void Threshold(double eps)
Replace small entries, abs(a_ij) &lt;= eps, with zero.
Definition: densemat.cpp:1803
int SizeI() const
Definition: densemat.hpp:832
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2182
void TestInversion()
Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
Definition: densemat.cpp:3273
double MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
Definition: densemat.cpp:786
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:630
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1359
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
Definition: densemat.cpp:2731
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:412
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2481
void BatchLUSolve(const DenseTensor &Mlu, const Array< int > &P, Vector &X)
Solve batch linear systems.
Definition: densemat.cpp:3649
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
DenseTensor(int i, int j, int k)
Definition: densemat.hpp:799
void ClearExternalData()
Definition: densemat.hpp:95
int CheckFinite() const
Definition: densemat.hpp:371
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2314
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:285
void Swap(DenseMatrix &other)
Definition: densemat.cpp:1902
void AddMatrix(DenseMatrix &A, int ro, int co)
Perform (ro+i,co+j)+=A(i,j) for 0&lt;=i&lt;A.Height, 0&lt;=j&lt;A.Width.
Definition: densemat.cpp:1648
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:921
void Clear()
Delete the matrix data array (if owned) and reset the matrix state.
Definition: densemat.hpp:98
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3252
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2249
double * HostReadWrite()
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:929
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 MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt = A D A^t, where D is diagonal.
Definition: densemat.cpp:2360
void Eigenvalues(DenseMatrix &b, Vector &ev)
Definition: densemat.hpp:258
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:401
virtual void PrintT(std::ostream &out=mfem::out, int width_=4) const
Prints the transpose matrix to stream out.
Definition: densemat.cpp:1862
void CopyCols(const DenseMatrix &A, int col1, int col2)
Copy columns col1 through col2 from A to *this.
Definition: densemat.cpp:1504
int SizeJ() const
Definition: densemat.hpp:833
virtual MatrixInverse * Inverse() const
Returns a pointer to the inverse matrix.
Definition: densemat.cpp:431
double a
Definition: lissajous.cpp:41
double * GetColumn(int col)
Definition: densemat.hpp:282
bool OwnsData() const
Return the DenseMatrix data (host pointer) ownership flag.
Definition: densemat.hpp:121
void AddMultADBt(const DenseMatrix &A, const Vector &D, const DenseMatrix &B, DenseMatrix &ADBt)
ADBt = A D B^t, where D is diagonal.
Definition: densemat.cpp:2538
virtual ~DenseMatrix()
Destroys dense matrix.
Definition: densemat.cpp:1909
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:1543
void AddMultTranspose(const Vector &x, Vector &y) const
y += A^t x
Definition: densemat.cpp:229
void CopyExceptMN(const DenseMatrix &A, int m, int n)
Copy All rows and columns except m and n from A.
Definition: densemat.cpp:1622
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:1329
void Mult(int m, int n, double *X) const
Definition: densemat.cpp:2924
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
bool Factor(int m, double TOL=0.0)
Compute the LU factorization of the current matrix.
Definition: densemat.cpp:2848
static const int ipiv_base
Definition: densemat.hpp:535
void GradToCurl(DenseMatrix &curl)
Definition: densemat.cpp:1417
DenseMatrixInverse()
Default constructor.
Definition: densemat.hpp:643
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
Definition: densemat.cpp:1207
void GetRowSums(Vector &l) const
Compute the row sums of the DenseMatrix.
Definition: densemat.cpp:1315
void CalcEigenvalues(double *lambda, double *vec) const
Definition: densemat.cpp:1232
const Memory< double > & GetMemory() const
Definition: densemat.hpp:118
DenseMatrix & RightSingularvectors()
Definition: densemat.hpp:776
RefCoord t[3]
void Eigenvalues(Vector &ev, DenseMatrix &evect)
Compute eigenvalues and eigenvectors of A x = ev x where A = *this.
Definition: densemat.hpp:249
DenseMatrixEigensystem(DenseMatrix &m)
Definition: densemat.cpp:3292
const double * HostRead() const
Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:913
int Rank(double tol) const
Definition: densemat.cpp:1192
const double alpha
Definition: ex15.cpp:369
LUFactors(double *data_, int *ipiv_)
Definition: densemat.hpp:544
void AddMult_a(double a, const Vector &x, Vector &y) const
y += a * A.x
Definition: densemat.cpp:247
void RightScaling(const Vector &s)
RightScaling: this = this * diag(s);.
Definition: densemat.cpp:328
DenseTensor(const DenseTensor &other)
Copy constructor: deep copy.
Definition: densemat.hpp:821
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:2633
Vector data type.
Definition: vector.hpp:60
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:160
void AddMultTranspose_a(double a, const Vector &x, Vector &y) const
y += a * A^t x
Definition: densemat.cpp:265
void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt += A D A^t, where D is diagonal.
Definition: densemat.cpp:2332
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:1719
Memory< double > & GetMemory()
Definition: densemat.hpp:896
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:1517
double * Write(bool on_dev=true)
Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:917
RefCoord s[3]
double * Write(bool on_dev=true)
Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:393
double * Data()
Definition: densemat.hpp:892
void InvLeftScaling(const Vector &s)
InvLeftScaling this = diag(1./s) * this.
Definition: densemat.cpp:315
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:253
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
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
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:2009
int Size() const
Get the size of the inverse matrix.
Definition: densemat.hpp:653
Vector & Singularvalues()
Definition: densemat.hpp:773
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:786
double InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Definition: densemat.hpp:182
double FNorm2() const
Compute the square of the Frobenius norm of the matrix.
Definition: densemat.hpp:242
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.cpp:150
void AdjustDofDirection(Array< int > &dofs)
Definition: densemat.cpp:1730
void GradToDiv(Vector &div)
Definition: densemat.cpp:1476
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:2690
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void Eigenvalues(DenseMatrix &b, Vector &ev, DenseMatrix &evect)
Compute generalized eigenvalues of A x = ev B x, where A = *this.
Definition: densemat.hpp:262
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:385
virtual void PrintMatlab(std::ostream &out=mfem::out) const
Prints operator in Matlab format.
Definition: densemat.cpp:1843
double * data
Definition: densemat.hpp:532
DenseMatrix & operator+=(const double *m)
Definition: densemat.cpp:573