MFEM  v4.5.1 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  /// Adds the matrix A multiplied by the number c to the matrix,
214  /// assuming A has the same dimensions and uses column-major layout.
215  void Add(const double c, const double *A);
216
217  /// Sets the matrix elements equal to constant c
218  DenseMatrix &operator=(double c);
219
220  /// Copy the matrix entries from the given array
221  DenseMatrix &operator=(const double *d);
222
223  /// Sets the matrix size and elements equal to those of m
224  DenseMatrix &operator=(const DenseMatrix &m);
225
226  DenseMatrix &operator+=(const double *m);
228
230
231  DenseMatrix &operator*=(double c);
232
233  /// (*this) = -(*this)
234  void Neg();
235
236  /// Take the 2-norm of the columns of A and store in v
237  void Norm2(double *v) const;
238
239  /// Compute the norm ||A|| = max_{ij} |A_{ij}|
240  double MaxMaxNorm() const;
241
242  /// Compute the Frobenius norm of the matrix
243  double FNorm() const { double s, n2; FNorm(s, n2); return s*sqrt(n2); }
244
245  /// Compute the square of the Frobenius norm of the matrix
246  double FNorm2() const { double s, n2; FNorm(s, n2); return s*s*n2; }
247
248  /// Compute eigenvalues of A x = ev x where A = *this
249  void Eigenvalues(Vector &ev)
250  { Eigensystem(ev); }
251
252  /// Compute eigenvalues and eigenvectors of A x = ev x where A = *this
253  void Eigenvalues(Vector &ev, DenseMatrix &evect)
254  { Eigensystem(ev, &evect); }
255
256  /// Compute eigenvalues and eigenvectors of A x = ev x where A = *this
257  void Eigensystem(Vector &ev, DenseMatrix &evect)
258  { Eigensystem(ev, &evect); }
259
260  /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
261  where A = *this */
263  { Eigensystem(b, ev); }
264
265  /// Compute generalized eigenvalues of A x = ev B x, where A = *this
267  { Eigensystem(b, ev, &evect); }
268
269  /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
270  where A = *this */
272  { Eigensystem(b, ev, &evect); }
273
274  void SingularValues(Vector &sv) const;
275  int Rank(double tol) const;
276
277  /// Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
278  double CalcSingularvalue(const int i) const;
279
280  /** Return the eigenvalues (in increasing order) and eigenvectors of a
281  2x2 or 3x3 symmetric matrix. */
282  void CalcEigenvalues(double *lambda, double *vec) const;
283
284  void GetRow(int r, Vector &row) const;
285  void GetColumn(int c, Vector &col) const;
286  double *GetColumn(int col) { return data + col*height; }
287  const double *GetColumn(int col) const { return data + col*height; }
288
289  void GetColumnReference(int c, Vector &col)
290  { col.SetDataAndSize(data + c * height, height); }
291
292  void SetRow(int r, const double* row);
293  void SetRow(int r, const Vector &row);
294
295  void SetCol(int c, const double* col);
296  void SetCol(int c, const Vector &col);
297
298
299  /// Set all entries of a row to the specified value.
300  void SetRow(int row, double value);
301  /// Set all entries of a column to the specified value.
302  void SetCol(int col, double value);
303
304  /// Returns the diagonal of the matrix
305  void GetDiag(Vector &d) const;
306  /// Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|
307  void Getl1Diag(Vector &l) const;
308  /// Compute the row sums of the DenseMatrix
309  void GetRowSums(Vector &l) const;
310
311  /// Creates n x n diagonal matrix with diagonal elements c
312  void Diag(double c, int n);
313  /// Creates n x n diagonal matrix with diagonal given by diag
314  void Diag(double *diag, int n);
315
316  /// (*this) = (*this)^t
317  void Transpose();
318  /// (*this) = A^t
319  void Transpose(const DenseMatrix &A);
320  /// (*this) = 1/2 ((*this) + (*this)^t)
321  void Symmetrize();
322
323  void Lump();
324
325  /** Given a DShape matrix (from a scalar FE), stored in *this, returns the
326  CurlShape matrix. If *this is a N by D matrix, then curl is a D*N by
327  D*(D-1)/2 matrix. The size of curl must be set outside. The dimension D
328  can be either 2 or 3. In 2D this computes the scalar-valued curl of a
329  2D vector field */
330  void GradToCurl(DenseMatrix &curl);
331  /** Given a DShape matrix (from a scalar FE), stored in *this, returns the
332  CurlShape matrix. This computes the vector-valued curl of a scalar field.
333  *this is N by 2 matrix and curl is N by 2 matrix as well. */
334  void GradToVectorCurl2D(DenseMatrix &curl);
335  /** Given a DShape matrix (from a scalar FE), stored in *this,
336  returns the DivShape vector. If *this is a N by dim matrix,
337  then div is a dim*N vector. The size of div must be set
338  outside. */
339  void GradToDiv(Vector &div);
340
341  /// Copy rows row1 through row2 from A to *this
342  void CopyRows(const DenseMatrix &A, int row1, int row2);
343  /// Copy columns col1 through col2 from A to *this
344  void CopyCols(const DenseMatrix &A, int col1, int col2);
345  /// Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this
346  void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco);
347  /// Copy matrix A to the location in *this at row_offset, col_offset
348  void CopyMN(const DenseMatrix &A, int row_offset, int col_offset);
349  /// Copy matrix A^t to the location in *this at row_offset, col_offset
350  void CopyMNt(const DenseMatrix &A, int row_offset, int col_offset);
351  /** Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this at
352  row_offset, col_offset */
353  void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco,
354  int row_offset, int col_offset);
355  /// Copy c on the diagonal of size n to *this at row_offset, col_offset
356  void CopyMNDiag(double c, int n, int row_offset, int col_offset);
357  /// Copy diag on the diagonal of size n to *this at row_offset, col_offset
358  void CopyMNDiag(double *diag, int n, int row_offset, int col_offset);
359  /// Copy All rows and columns except m and n from A
360  void CopyExceptMN(const DenseMatrix &A, int m, int n);
361
362  /// Perform (ro+i,co+j)+=A(i,j) for 0<=i<A.Height, 0<=j<A.Width
363  void AddMatrix(DenseMatrix &A, int ro, int co);
364  /// Perform (ro+i,co+j)+=a*A(i,j) for 0<=i<A.Height, 0<=j<A.Width
365  void AddMatrix(double a, const DenseMatrix &A, int ro, int co);
366
367  /** Get the square submatrix which corresponds to the given indices @a idx.
368  Note: the @a A matrix will be resized to accommodate the data */
369  void GetSubMatrix(const Array<int> & idx, DenseMatrix & A) const;
370
371  /** Get the rectangular submatrix which corresponds to the given indices
372  @a idx_i and @a idx_j. Note: the @a A matrix will be resized to
373  accommodate the data */
374  void GetSubMatrix(const Array<int> & idx_i, const Array<int> & idx_j,
375  DenseMatrix & A) const;
376
377  /** Get the square submatrix which corresponds to the range
378  [ @a ibeg, @a iend ). Note: the @a A matrix will be resized
379  to accommodate the data */
380  void GetSubMatrix(int ibeg, int iend, DenseMatrix & A);
381
382  /** Get the square submatrix which corresponds to the range
383  i âˆˆ [ @a ibeg, @a iend ) and j âˆˆ [ @a jbeg, @a jend )
384  Note: the @a A matrix will be resized to accommodate the data */
385  void GetSubMatrix(int ibeg, int iend, int jbeg, int jend, DenseMatrix & A);
386
387  /// Set (*this)(idx[i],idx[j]) = A(i,j)
388  void SetSubMatrix(const Array<int> & idx, const DenseMatrix & A);
389
390  /// Set (*this)(idx_i[i],idx_j[j]) = A(i,j)
391  void SetSubMatrix(const Array<int> & idx_i, const Array<int> & idx_j,
392  const DenseMatrix & A);
393
394  /** Set a submatrix of (this) to the given matrix @a A
395  with row and column offset @a ibeg */
396  void SetSubMatrix(int ibeg, const DenseMatrix & A);
397
398  /** Set a submatrix of (this) to the given matrix @a A
399  with row and column offset @a ibeg and @a jbeg respectively */
400  void SetSubMatrix(int ibeg, int jbeg, const DenseMatrix & A);
401
402  /// (*this)(idx[i],idx[j]) += A(i,j)
403  void AddSubMatrix(const Array<int> & idx, const DenseMatrix & A);
404
405  /// (*this)(idx_i[i],idx_j[j]) += A(i,j)
406  void AddSubMatrix(const Array<int> & idx_i, const Array<int> & idx_j,
407  const DenseMatrix & A);
408
409  /** Add the submatrix @a A to this with row and column offset @a ibeg */
410  void AddSubMatrix(int ibeg, const DenseMatrix & A);
411
412  /** Add the submatrix @a A to this with row and column offsets
413  @a ibeg and @a jbeg respectively */
414  void AddSubMatrix(int ibeg, int jbeg, const DenseMatrix & A);
415
416  /// Add the matrix 'data' to the Vector 'v' at the given 'offset'
417  void AddToVector(int offset, Vector &v) const;
418  /// Get the matrix 'data' from the Vector 'v' at the given 'offset'
419  void GetFromVector(int offset, const Vector &v);
420  /** If (dofs[i] < 0 and dofs[j] >= 0) or (dofs[i] >= 0 and dofs[j] < 0)
421  then (*this)(i,j) = -(*this)(i,j). */
422  void AdjustDofDirection(Array<int> &dofs);
423
424  /// Replace small entries, abs(a_ij) <= eps, with zero.
425  void Threshold(double eps);
426
427  /** Count the number of entries in the matrix for which isfinite
428  is false, i.e. the entry is a NaN or +/-Inf. */
429  int CheckFinite() const { return mfem::CheckFinite(HostRead(), height*width); }
430
431  /// Prints matrix to stream out.
432  virtual void Print(std::ostream &out = mfem::out, int width_ = 4) const;
433  virtual void PrintMatlab(std::ostream &out = mfem::out) const;
434  /// Prints the transpose matrix to stream out.
435  virtual void PrintT(std::ostream &out = mfem::out, int width_ = 4) const;
436
437  /// Invert and print the numerical conditioning of the inversion.
438  void TestInversion();
439
440  std::size_t MemoryUsage() const { return data.Capacity() * sizeof(double); }
441
442  /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
443  const double *Read(bool on_dev = true) const
444  { return mfem::Read(data, Height()*Width(), on_dev); }
445
446  /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
447  const double *HostRead() const
448  { return mfem::Read(data, Height()*Width(), false); }
449
450  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
451  double *Write(bool on_dev = true)
452  { return mfem::Write(data, Height()*Width(), on_dev); }
453
454  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
455  double *HostWrite()
456  { return mfem::Write(data, Height()*Width(), false); }
457
458  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
459  double *ReadWrite(bool on_dev = true)
460  { return mfem::ReadWrite(data, Height()*Width(), on_dev); }
461
462  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
464  { return mfem::ReadWrite(data, Height()*Width(), false); }
465
466  void Swap(DenseMatrix &other);
467
468  /// Destroys dense matrix.
469  virtual ~DenseMatrix();
470 };
471
472 /// C = A + alpha*B
473 void Add(const DenseMatrix &A, const DenseMatrix &B,
474  double alpha, DenseMatrix &C);
475
476 /// C = alpha*A + beta*B
477 void Add(double alpha, const double *A,
478  double beta, const double *B, DenseMatrix &C);
479
480 /// C = alpha*A + beta*B
481 void Add(double alpha, const DenseMatrix &A,
482  double beta, const DenseMatrix &B, DenseMatrix &C);
483
484 /// @brief Solves the dense linear system, A * X = B for X
485 ///
486 /// @param [in,out] A the square matrix for the linear system
487 /// @param [in,out] X the rhs vector, B, on input, the solution, X, on output.
488 /// @param [in] TOL optional fuzzy comparison tolerance. Defaults to 1e-9.
489 ///
490 /// @return status set to true if successful, otherwise, false.
491 ///
492 /// @note This routine may replace the contents of the input Matrix, A, with the
493 /// corresponding LU factorization of the matrix. Matrices of size 1x1 and
494 /// 2x2 are handled explicitly.
495 ///
496 /// @pre A.IsSquare() == true
497 /// @pre X != nullptr
498 bool LinearSolve(DenseMatrix& A, double* X, double TOL = 1.e-9);
499
500 /// Matrix matrix multiplication. A = B * C.
501 void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);
502
503 /// Matrix matrix multiplication. A += B * C.
504 void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);
505
506 /// Matrix matrix multiplication. A += alpha * B * C.
507 void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c,
508  DenseMatrix &a);
509
510 /** Calculate the adjugate of a matrix (for NxN matrices, N=1,2,3) or the matrix
511  adj(A^t.A).A^t for rectangular matrices (2x1, 3x1, or 3x2). This operation
512  is well defined even when the matrix is not full rank. */
513 void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja);
514
515 /// Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
516 void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat);
517
518 /** Calculate the inverse of a matrix (for NxN matrices, N=1,2,3) or the
519  left inverse (A^t.A)^{-1}.A^t (for 2x1, 3x1, or 3x2 matrices) */
520 void CalcInverse(const DenseMatrix &a, DenseMatrix &inva);
521
522 /// Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
523 void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva);
524
525 /** For a given Nx(N-1) (N=2,3) matrix J, compute a vector n such that
526  n_k = (-1)^{k+1} det(J_k), k=1,..,N, where J_k is the matrix J with the
527  k-th row removed. Note: J^t.n = 0, det([n|J])=|n|^2=det(J^t.J). */
528 void CalcOrtho(const DenseMatrix &J, Vector &n);
529
530 /// Calculate the matrix A.At
531 void MultAAt(const DenseMatrix &a, DenseMatrix &aat);
532
533 /// ADAt = A D A^t, where D is diagonal
534 void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt);
535
536 /// ADAt += A D A^t, where D is diagonal
537 void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt);
538
539 /// Multiply a matrix A with the transpose of a matrix B: A*Bt
540 void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);
541
542 /// ADBt = A D B^t, where D is diagonal
543 void MultADBt(const DenseMatrix &A, const Vector &D,
544  const DenseMatrix &B, DenseMatrix &ADBt);
545
546 /// ABt += A * B^t
547 void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);
548
549 /// ADBt = A D B^t, where D is diagonal
550 void AddMultADBt(const DenseMatrix &A, const Vector &D,
551  const DenseMatrix &B, DenseMatrix &ADBt);
552
553 /// ABt += a * A * B^t
554 void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B,
555  DenseMatrix &ABt);
556
557 /// Multiply the transpose of a matrix A with a matrix B: At*B
558 void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB);
559
560 /// AAt += a * A * A^t
561 void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);
562
563 /// AAt = a * A * A^t
564 void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);
565
566 /// Make a matrix from a vector V.Vt
567 void MultVVt(const Vector &v, DenseMatrix &vvt);
568
569 void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);
570
571 /// VWt += v w^t
572 void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);
573
574 /// VVt += v v^t
575 void AddMultVVt(const Vector &v, DenseMatrix &VWt);
576
577 /// VWt += a * v w^t
578 void AddMult_a_VWt(const double a, const Vector &v, const Vector &w,
579  DenseMatrix &VWt);
580
581 /// VVt += a * v v^t
582 void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt);
583
584 /** Computes matrix P^t * A * P. Note: The @a RAP matrix will be resized
585  to accommodate the data */
586 void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix & RAP);
587
588 /** Computes the matrix Rt^t * A * P. Note: The @a RAP matrix will be resized
589  to accommodate the data */
590 void RAP(const DenseMatrix &Rt, const DenseMatrix &A,
591  const DenseMatrix &P, DenseMatrix & RAP);
592
593 /** Abstract class that can compute factorization of external data and perform various
594  operations with the factored data. */
595 class Factors
596 {
597 public:
598
599  double *data;
600
601  Factors() { }
602
603  Factors(double *data_) : data(data_) { }
604
605  virtual bool Factor(int m, double TOL = 0.0)
606  {
607  mfem_error("Factors::Factors(...)");
608  return false;
609  }
610
611  virtual double Det(int m) const
612  {
613  mfem_error("Factors::Det(...)");
614  return 0.;
615  }
616
617  virtual void Solve(int m, int n, double *X) const
618  {
619  mfem_error("Factors::Solve(...)");
620  }
621
622  virtual void GetInverseMatrix(int m, double *X) const
623  {
624  mfem_error("Factors::GetInverseMatrix(...)");
625  }
626
627  virtual ~Factors() {}
628 };
629
630
631 /** Class that can compute LU factorization of external data and perform various
632  operations with the factored data. */
633 class LUFactors : public Factors
634 {
635 public:
636  int *ipiv;
637 #ifdef MFEM_USE_LAPACK
638  static const int ipiv_base = 1;
639 #else
640  static const int ipiv_base = 0;
641 #endif
642
643  /** With this constructor, the (public) data and ipiv members should be set
644  explicitly before calling class methods. */
646
647  LUFactors(double *data_, int *ipiv_) : Factors(data_), ipiv(ipiv_) { }
648
649  /**
650  * @brief Compute the LU factorization of the current matrix
651  *
652  * Factorize the current matrix of size (m x m) overwriting it with the
653  * LU factors. The factorization is such that L.U = P.A, where A is the
654  * original matrix and P is a permutation matrix represented by ipiv.
655  *
656  * @param [in] m size of the square matrix
657  * @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0.
658  *
659  * @return status set to true if successful, otherwise, false.
660  */
661  virtual bool Factor(int m, double TOL = 0.0);
662
663  /** Assuming L.U = P.A factored data of size (m x m), compute |A|
664  from the diagonal values of U and the permutation information. */
665  virtual double Det(int m) const;
666
667  /** Assuming L.U = P.A factored data of size (m x m), compute X <- A X,
668  for a matrix X of size (m x n). */
669  void Mult(int m, int n, double *X) const;
670
671  /** Assuming L.U = P.A factored data of size (m x m), compute
672  X <- L^{-1} P X, for a matrix X of size (m x n). */
673  void LSolve(int m, int n, double *X) const;
674
675  /** Assuming L.U = P.A factored data of size (m x m), compute
676  X <- U^{-1} X, for a matrix X of size (m x n). */
677  void USolve(int m, int n, double *X) const;
678
679  /** Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1} X,
680  for a matrix X of size (m x n). */
681  virtual void Solve(int m, int n, double *X) const;
682
683  /** Assuming L.U = P.A factored data of size (m x m), compute X <- X A^{-1},
684  for a matrix X of size (n x m). */
685  void RightSolve(int m, int n, double *X) const;
686
687  /// Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1}.
688  virtual void GetInverseMatrix(int m, double *X) const;
689
690  /** Given an (n x m) matrix A21, compute X2 <- X2 - A21 X1, for matrices X1,
691  and X2 of size (m x r) and (n x r), respectively. */
692  static void SubMult(int m, int n, int r, const double *A21,
693  const double *X1, double *X2);
694
695  /** Assuming P.A = L.U factored data of size (m x m), compute the 2x2 block
696  decomposition:
697  | P 0 | | A A12 | = | L 0 | | U U12 |
698  | 0 I | | A21 A22 | | L21 I | | 0 S22 |
699  where A12, A21, and A22 are matrices of size (m x n), (n x m), and
700  (n x n), respectively. The blocks are overwritten as follows:
701  A12 <- U12 = L^{-1} P A12
702  A21 <- L21 = A21 U^{-1}
703  A22 <- S22 = A22 - L21 U12.
704  The block S22 is the Schur complement. */
705  void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const;
706
707  /** Given BlockFactor()'d data, perform the forward block solve for the
708  linear system:
709  | A A12 | | X1 | = | B1 |
710  | A21 A22 | | X2 | | B2 |
711  written in the factored form:
712  | L 0 | | U U12 | | X1 | = | P 0 | | B1 |
713  | L21 I | | 0 S22 | | X2 | | 0 I | | B2 |.
714  The resulting blocks Y1, Y2 solve the system:
715  | L 0 | | Y1 | = | P 0 | | B1 |
716  | L21 I | | Y2 | | 0 I | | B2 |
717  The blocks are overwritten as follows:
718  B1 <- Y1 = L^{-1} P B1
719  B2 <- Y2 = B2 - L21 Y1 = B2 - A21 A^{-1} B1
720  The blocks B1/Y1 and B2/Y2 are of size (m x r) and (n x r), respectively.
721  The Schur complement system is given by: S22 X2 = Y2. */
722  void BlockForwSolve(int m, int n, int r, const double *L21,
723  double *B1, double *B2) const;
724
725  /** Given BlockFactor()'d data, perform the backward block solve in
726  | U U12 | | X1 | = | Y1 |
727  | 0 S22 | | X2 | | Y2 |.
728  The input is the solution block X2 and the block Y1 resulting from
729  BlockForwSolve(). The result block X1 overwrites input block Y1:
730  Y1 <- X1 = U^{-1} (Y1 - U12 X2). */
731  void BlockBackSolve(int m, int n, int r, const double *U12,
732  const double *X2, double *Y1) const;
733 };
734
735
736 /** Class that can compute Cholesky factorizations of external data of an
737  SPD matrix and perform various operations with the factored data. */
738 class CholeskyFactors : public Factors
739 {
740 public:
741
742  /** With this constructor, the (public) data should be set
743  explicitly before calling class methods. */
745
746  CholeskyFactors(double *data_) : Factors(data_) { }
747
748  /**
749  * @brief Compute the Cholesky factorization of the current matrix
750  *
751  * Factorize the current matrix of size (m x m) overwriting it with the
752  * Cholesky factors. The factorization is such that LL^t = A, where A is the
753  * original matrix
754  *
755  * @param [in] m size of the square matrix
756  * @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0.
757  *
758  * @return status set to true if successful, otherwise, false.
759  */
760  virtual bool Factor(int m, double TOL = 0.0);
761
762  /** Assuming LL^t = A factored data of size (m x m), compute |A|
763  from the diagonal values of L */
764  virtual double Det(int m) const;
765
766  /** Assuming L.L^t = A factored data of size (m x m), compute X <- L X,
767  for a matrix X of size (m x n). */
768  void LMult(int m, int n, double *X) const;
769
770  /** Assuming L.L^t = A factored data of size (m x m), compute X <- L^t X,
771  for a matrix X of size (m x n). */
772  void UMult(int m, int n, double *X) const;
773
774  /** Assuming L L^t = A factored data of size (m x m), compute
775  X <- L^{-1} X, for a matrix X of size (m x n). */
776  void LSolve(int m, int n, double *X) const;
777
778  /** Assuming L L^t = A factored data of size (m x m), compute
779  X <- L^{-t} X, for a matrix X of size (m x n). */
780  void USolve(int m, int n, double *X) const;
781
782  /** Assuming L.L^t = A factored data of size (m x m), compute X <- A^{-1} X,
783  for a matrix X of size (m x n). */
784  virtual void Solve(int m, int n, double *X) const;
785
786  /** Assuming L.L^t = A factored data of size (m x m), compute X <- X A^{-1},
787  for a matrix X of size (n x m). */
788  void RightSolve(int m, int n, double *X) const;
789
790  /// Assuming L.L^t = A factored data of size (m x m), compute X <- A^{-1}.
791  virtual void GetInverseMatrix(int m, double *X) const;
792
793 };
794
795
796 /** Data type for inverse of square dense matrix.
797  Stores matrix factors, i.e., Cholesky factors if the matrix is SPD,
798  LU otherwise. */
800 {
801 private:
802  const DenseMatrix *a;
803  Factors * factors = nullptr;
804  bool spd = false;
805
806  void Init(int m);
807  bool own_data = false;
808 public:
809  /// Default constructor.
810  DenseMatrixInverse(bool spd_=false) : a(NULL), spd(spd_) { Init(0); }
811
812  /** Creates square dense matrix. Computes factorization of mat
813  and stores its factors. */
814  DenseMatrixInverse(const DenseMatrix &mat, bool spd_ = false);
815
816  /// Same as above but does not factorize the matrix.
817  DenseMatrixInverse(const DenseMatrix *mat, bool spd_ = false);
818
819  /// Get the size of the inverse matrix
820  int Size() const { return Width(); }
821
822  /// Factor the current DenseMatrix, *a
823  void Factor();
824
825  /// Factor a new DenseMatrix of the same size
826  void Factor(const DenseMatrix &mat);
827
828  virtual void SetOperator(const Operator &op);
829
830  /// Matrix vector multiplication with the inverse of dense matrix.
831  void Mult(const double *x, double *y) const;
832
833  /// Matrix vector multiplication with the inverse of dense matrix.
834  virtual void Mult(const Vector &x, Vector &y) const;
835
836  /// Multiply the inverse matrix by another matrix: X = A^{-1} B.
837  void Mult(const DenseMatrix &B, DenseMatrix &X) const;
838
839  /// Multiply the inverse matrix by another matrix: X <- A^{-1} X.
840  void Mult(DenseMatrix &X) const {factors->Solve(width, X.Width(), X.Data());}
841
842  /// Compute and return the inverse matrix in Ainv.
843  void GetInverseMatrix(DenseMatrix &Ainv) const;
844
845  /// Compute the determinant of the original DenseMatrix using the LU factors.
846  double Det() const { return factors->Det(width); }
847
848  /// Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
849  void TestInversion();
850
851  /// Destroys dense inverse matrix.
852  virtual ~DenseMatrixInverse();
853 };
854
855 #ifdef MFEM_USE_LAPACK
856
858 {
859  DenseMatrix &mat;
860  Vector EVal;
861  DenseMatrix EVect;
862  Vector ev;
863  int n;
864  double *work;
865  char jobz, uplo;
866  int lwork, info;
867 public:
868
871  void Eval();
872  Vector &Eigenvalues() { return EVal; }
873  DenseMatrix &Eigenvectors() { return EVect; }
874  double Eigenvalue(int i) { return EVal(i); }
875  const Vector &Eigenvector(int i)
876  {
877  ev.SetData(EVect.Data() + i * EVect.Height());
878  return ev;
879  }
881 };
882
884 {
885  DenseMatrix &A;
886  DenseMatrix &B;
887  DenseMatrix A_copy;
888  DenseMatrix B_copy;
889  Vector evalues_r;
890  Vector evalues_i;
891  DenseMatrix Vr;
892  DenseMatrix Vl;
893  int n;
894
895  double *alphar;
896  double *alphai;
897  double *beta;
898  double *work;
899  char jobvl, jobvr;
900  int lwork, info;
901
902 public:
903
905  bool left_eigen_vectors = false,
906  bool right_eigen_vectors = false);
907  void Eval();
908  Vector &EigenvaluesRealPart() { return evalues_r; }
909  Vector &EigenvaluesImagPart() { return evalues_i; }
910  double EigenvalueRealPart(int i) { return evalues_r(i); }
911  double EigenvalueImagPart(int i) { return evalues_i(i); }
912  DenseMatrix &LeftEigenvectors() { return Vl; }
913  DenseMatrix &RightEigenvectors() { return Vr; }
915 };
916
917
919 {
920  DenseMatrix Mc;
921  Vector sv;
922  DenseMatrix U,Vt;
923  int m, n;
924
925 #ifdef MFEM_USE_LAPACK
926  double *work;
927  char jobu, jobvt;
928  int lwork, info;
929 #endif
930
931  void Init();
932 public:
934  bool left_singular_vectors=false,
935  bool right_singlular_vectors=false);
936  DenseMatrixSVD(int h, int w,
937  bool left_singular_vectors=false,
938  bool right_singlular_vectors=false);
939  void Eval(DenseMatrix &M);
940  Vector &Singularvalues() { return sv; }
941  double Singularvalue(int i) { return sv(i); }
944  ~DenseMatrixSVD();
945 };
946
947 #endif // if MFEM_USE_LAPACK
948
949
950 class Table;
951
952 /// Rank 3 tensor (array of matrices)
954 {
955 private:
956  mutable DenseMatrix Mk;
957  Memory<double> tdata;
958  int nk;
959
960 public:
962  {
963  nk = 0;
964  }
965
966  DenseTensor(int i, int j, int k)
967  : Mk(NULL, i, j)
968  {
969  nk = k;
970  tdata.New(i*j*k);
971  }
972
973  DenseTensor(double *d, int i, int j, int k)
974  : Mk(NULL, i, j)
975  {
976  nk = k;
977  tdata.Wrap(d, i*j*k, false);
978  }
979
980  DenseTensor(int i, int j, int k, MemoryType mt)
981  : Mk(NULL, i, j)
982  {
983  nk = k;
984  tdata.New(i*j*k, mt);
985  }
986
987  /// Copy constructor: deep copy
988  DenseTensor(const DenseTensor &other)
989  : Mk(NULL, other.Mk.height, other.Mk.width), nk(other.nk)
990  {
991  const int size = Mk.Height()*Mk.Width()*nk;
992  if (size > 0)
993  {
994  tdata.New(size, other.tdata.GetMemoryType());
995  tdata.CopyFrom(other.tdata, size);
996  }
997  }
998
999  int SizeI() const { return Mk.Height(); }
1000  int SizeJ() const { return Mk.Width(); }
1001  int SizeK() const { return nk; }
1002
1003  int TotalSize() const { return SizeI()*SizeJ()*SizeK(); }
1004
1005  void SetSize(int i, int j, int k, MemoryType mt_ = MemoryType::PRESERVE)
1006  {
1007  const MemoryType mt = mt_ == MemoryType::PRESERVE ? tdata.GetMemoryType() : mt_;
1008  tdata.Delete();
1009  Mk.UseExternalData(NULL, i, j);
1010  nk = k;
1011  tdata.New(i*j*k, mt);
1012  }
1013
1014  void UseExternalData(double *ext_data, int i, int j, int k)
1015  {
1016  tdata.Delete();
1017  Mk.UseExternalData(NULL, i, j);
1018  nk = k;
1019  tdata.Wrap(ext_data, i*j*k, false);
1020  }
1021
1022  /// Sets the tensor elements equal to constant c
1023  DenseTensor &operator=(double c);
1024
1025  /// Copy assignment operator (performs a deep copy)
1026  DenseTensor &operator=(const DenseTensor &other);
1027
1029  {
1030  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1031  Mk.data = Memory<double>(GetData(k), SizeI()*SizeJ(), false);
1032  return Mk;
1033  }
1034  const DenseMatrix &operator()(int k) const
1035  {
1036  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1037  Mk.data = Memory<double>(const_cast<double*>(GetData(k)), SizeI()*SizeJ(),
1038  false);
1039  return Mk;
1040  }
1041
1042  double &operator()(int i, int j, int k)
1043  {
1044  MFEM_ASSERT_INDEX_IN_RANGE(i, 0, SizeI());
1045  MFEM_ASSERT_INDEX_IN_RANGE(j, 0, SizeJ());
1046  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1047  return tdata[i+SizeI()*(j+SizeJ()*k)];
1048  }
1049
1050  const double &operator()(int i, int j, int k) const
1051  {
1052  MFEM_ASSERT_INDEX_IN_RANGE(i, 0, SizeI());
1053  MFEM_ASSERT_INDEX_IN_RANGE(j, 0, SizeJ());
1054  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1055  return tdata[i+SizeI()*(j+SizeJ()*k)];
1056  }
1057
1058  double *GetData(int k)
1059  {
1060  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1061  return tdata+k*Mk.Height()*Mk.Width();
1062  }
1063
1064  const double *GetData(int k) const
1065  {
1066  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
1067  return tdata+k*Mk.Height()*Mk.Width();
1068  }
1069
1070  double *Data() { return tdata; }
1071
1072  const double *Data() const { return tdata; }
1073
1074  Memory<double> &GetMemory() { return tdata; }
1075  const Memory<double> &GetMemory() const { return tdata; }
1076
1077  /** Matrix-vector product from unassembled element matrices, assuming both
1078  'x' and 'y' use the same elem_dof table. */
1079  void AddMult(const Table &elem_dof, const Vector &x, Vector &y) const;
1080
1081  void Clear()
1082  { UseExternalData(NULL, 0, 0, 0); }
1083
1084  std::size_t MemoryUsage() const { return nk*Mk.MemoryUsage(); }
1085
1086  /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
1087  const double *Read(bool on_dev = true) const
1088  { return mfem::Read(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
1089
1090  /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
1091  const double *HostRead() const
1092  { return mfem::Read(tdata, Mk.Height()*Mk.Width()*nk, false); }
1093
1094  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
1095  double *Write(bool on_dev = true)
1096  { return mfem::Write(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
1097
1098  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
1099  double *HostWrite()
1100  { return mfem::Write(tdata, Mk.Height()*Mk.Width()*nk, false); }
1101
1102  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
1103  double *ReadWrite(bool on_dev = true)
1104  { return mfem::ReadWrite(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
1105
1106  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
1108  { return mfem::ReadWrite(tdata, Mk.Height()*Mk.Width()*nk, false); }
1109
1111  {
1112  mfem::Swap(tdata, t.tdata);
1113  mfem::Swap(nk, t.nk);
1114  Mk.Swap(t.Mk);
1115  }
1116
1117  ~DenseTensor() { tdata.Delete(); }
1118 };
1119
1120 /** @brief Compute the LU factorization of a batch of matrices
1121
1122  Factorize n matrices of size (m x m) stored in a dense tensor overwriting it
1123  with the LU factors. The factorization is such that L.U = Piv.A, where A is
1124  the original matrix and Piv is a permutation matrix represented by P.
1125
1126  @param [in, out] Mlu batch of square matrices - dimension m x m x n.
1127  @param [out] P array storing pivot information - dimension m x n.
1128  @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0. */
1129 void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const double TOL = 0.0);
1130
1131 /** @brief Solve batch linear systems
1132
1133  Assuming L.U = P.A for n factored matrices (m x m), compute x <- A x, for n
1134  companion vectors.
1135
1136  @param [in] Mlu batch of LU factors for matrix M - dimension m x m x n.
1137  @param [in] P array storing pivot information - dimension m x n.
1138  @param [in, out] X vector storing right-hand side and then solution -
1139  dimension m x n. */
1140 void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X);
1141
1142
1143 // Inline methods
1144
1145 inline double &DenseMatrix::operator()(int i, int j)
1146 {
1147  MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
1148  return data[i+j*height];
1149 }
1150
1151 inline const double &DenseMatrix::operator()(int i, int j) const
1152 {
1153  MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
1154  return data[i+j*height];
1155 }
1156
1157 } // namespace mfem
1158
1159 #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:1414
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:2742
DenseMatrix & operator-=(const DenseMatrix &m)
Definition: densemat.cpp:608
void SymmetricScaling(const Vector &s)
SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
Definition: densemat.cpp:370
void SquareRootInverse()
Replaces the current matrix with its square root inverse.
Definition: densemat.cpp:750
int CheckFinite(const double *v, const int n)
Definition: vector.hpp:493
void BatchLUFactor(DenseTensor &Mlu, Array< int > &P, const double TOL)
Compute the LU factorization of a batch of matrices.
Definition: densemat.cpp:4260
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:3127
void Mult(DenseMatrix &X) const
Multiply the inverse matrix by another matrix: X &lt;- A^{-1} X.
Definition: densemat.hpp:840
DenseMatrix & operator*=(double c)
Definition: densemat.cpp:621
void GetDiag(Vector &d) const
Returns the diagonal of the matrix.
Definition: densemat.cpp:1306
void SetCol(int c, const double *col)
Definition: densemat.cpp:2154
void UseExternalData(double *ext_data, int i, int j, int k)
Definition: densemat.hpp:1014
DenseTensor & operator=(double c)
Sets the tensor elements equal to constant c.
Definition: densemat.cpp:4243
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:3108
double * HostWrite()
Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:455
Memory< double > & GetMemory()
Definition: densemat.hpp:117
void SetRow(int r, const double *row)
Definition: densemat.cpp:2139
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
Definition: densemat.cpp:356
const DenseMatrix & operator()(int k) const
Definition: densemat.hpp:1034
virtual ~Factors()
Definition: densemat.hpp:627
void Eigenvalues(Vector &ev)
Compute eigenvalues of A x = ev x where A = *this.
Definition: densemat.hpp:249
void SingularValues(Vector &sv) const
Definition: densemat.cpp:1174
void Delete()
Delete the owned pointers and reset the Memory object.
DenseMatrix & operator()(int k)
Definition: densemat.hpp:1028
double Det() const
Definition: densemat.cpp:449
void Eigensystem(DenseMatrix &b, Vector &ev, DenseMatrix &evect)
Definition: densemat.hpp:271
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:73
int SizeK() const
Definition: densemat.hpp:1001
virtual void GetInverseMatrix(int m, double *X) const
Assuming L.L^t = A factored data of size (m x m), compute X &lt;- A^{-1}.
Definition: densemat.cpp:3785
void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const
Definition: densemat.cpp:3530
void BlockBackSolve(int m, int n, int r, const double *U12, const double *X2, double *Y1) const
Definition: densemat.cpp:3565
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
Definition: densemat.cpp:297
std::size_t MemoryUsage() const
Definition: densemat.hpp:440
Definition: densemat.cpp:2440
DenseMatrixSVD(DenseMatrix &M, bool left_singular_vectors=false, bool right_singlular_vectors=false)
Definition: densemat.cpp:4115
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:1087
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:4188
void TestInversion()
Invert and print the numerical conditioning of the inversion.
Definition: densemat.cpp:2254
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
DenseMatrix & LeftSingularvectors()
Definition: densemat.hpp:942
void CopyRows(const DenseMatrix &A, int row1, int row2)
Copy rows row1 through row2 from A to *this.
Definition: densemat.cpp:1527
virtual double Det(int m) const
Definition: densemat.cpp:3289
void Swap(DenseTensor &t)
Definition: densemat.hpp:1110
virtual void GetInverseMatrix(int m, double *X) const
Definition: densemat.hpp:622
void Eval(DenseMatrix &M)
Definition: densemat.cpp:4149
const double * Data() const
Definition: densemat.hpp:1072
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:2940
void LMult(int m, int n, double *X) const
Definition: densemat.cpp:3627
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:463
Factors(double *data_)
Definition: densemat.hpp:603
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:3890
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3878
void CopyFrom(const Memory &src, int size)
Copy size entries from src to *this.
Definition: densemat.cpp:1498
virtual 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:3450
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
const Memory< double > & GetMemory() const
Definition: densemat.hpp:1075
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2654
DenseMatrix & operator=(double c)
Sets the matrix elements equal to constant c.
Definition: densemat.cpp:562
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:533
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:3081
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:3513
void USolve(int m, int n, double *X) const
Definition: densemat.cpp:3694
const Vector & Eigenvector(int i)
Definition: densemat.hpp:875
virtual void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: densemat.cpp:2183
const double & operator()(int i, int j, int k) const
Definition: densemat.hpp:1050
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:2282
double & operator()(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.hpp:1145
double Weight() const
Definition: densemat.cpp:506
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:3363
double FNorm() const
Compute the Frobenius norm of the matrix.
Definition: densemat.hpp:243
const double * HostRead() const
Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:447
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:201
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2512
double & operator()(int i, int j, int k)
Definition: densemat.hpp:1042
void RightSolve(int m, int n, double *X) const
Definition: densemat.cpp:3395
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2408
double operator*(const DenseMatrix &m) const
Matrix inner product: tr(A^t B)
Definition: densemat.cpp:186
double Singularvalue(int i)
Definition: densemat.hpp:941
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:67
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:542
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:3172
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:398
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1276
void BlockForwSolve(int m, int n, int r, const double *L21, double *B1, double *B2) const
Definition: densemat.cpp:3556
Abstract data type matrix.
Definition: matrix.hpp:27
const double * GetColumn(int col) const
Definition: densemat.hpp:287
void Norm2(double *v) const
Take the 2-norm of the columns of A and store in v.
Definition: densemat.cpp:795
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:2807
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:640
bool LinearSolve(DenseMatrix &A, double *X, double TOL)
Solves the dense linear system, A * X = B for X
Definition: densemat.cpp:2304
virtual ~DenseMatrixInverse()
Destroys dense inverse matrix.
Definition: densemat.cpp:3956
void LSolve(int m, int n, double *X) const
Definition: densemat.cpp:3340
DenseTensor(int i, int j, int k, MemoryType mt)
Definition: densemat.hpp:980
int TotalSize() const
Definition: densemat.hpp:1003
DenseTensor(double *d, int i, int j, int k)
Definition: densemat.hpp:973
DenseMatrixGeneralizedEigensystem(DenseMatrix &a, DenseMatrix &b, bool left_eigen_vectors=false, bool right_eigen_vectors=false)
Definition: densemat.cpp:4029
void LeftScaling(const Vector &s)
LeftScaling this = diag(s) * this.
Definition: densemat.cpp:315
double Det() const
Compute the determinant of the original DenseMatrix using the LU factors.
Definition: densemat.hpp:846
void AddMultVVt(const Vector &v, DenseMatrix &VVt)
VVt += v v^t.
Definition: densemat.cpp:3148
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:1623
double * GetData(int k)
Definition: densemat.hpp:1058
virtual double Det(int m) const
Definition: densemat.cpp:3617
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:3194
void Neg()
(*this) = -(*this)
Definition: densemat.cpp:631
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition: densemat.hpp:1005
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: densemat.cpp:3917
virtual void Solve(int m, int n, double *X) const
Definition: densemat.cpp:3381
double * data
Definition: densemat.hpp:599
double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:1103
void Getl1Diag(Vector &l) const
Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|.
Definition: densemat.cpp:1320
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3213
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:2074
void SetData(double *d)
Definition: vector.hpp:150
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
DenseMatrix & Eigenvectors()
Definition: densemat.hpp:873
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:1292
void AddMult(const Vector &x, Vector &y) const
y += A.x
Definition: densemat.cpp:224
void Threshold(double eps)
Replace small entries, abs(a_ij) &lt;= eps, with zero.
Definition: densemat.cpp:2169
int SizeI() const
Definition: densemat.hpp:999
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2548
void TestInversion()
Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
Definition: densemat.cpp:3945
double MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
Definition: densemat.cpp:808
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:630
void LSolve(int m, int n, double *X) const
Definition: densemat.cpp:3664
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1381
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
Definition: densemat.cpp:3097
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:425
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2847
void BatchLUSolve(const DenseTensor &Mlu, const Array< int > &P, Vector &X)
Solve batch linear systems.
Definition: densemat.cpp:4327
virtual bool Factor(int m, double TOL=0.0)
Compute the Cholesky factorization of the current matrix.
Definition: densemat.cpp:3575
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
DenseTensor(int i, int j, int k)
Definition: densemat.hpp:966
void ClearExternalData()
Definition: densemat.hpp:95
int CheckFinite() const
Definition: densemat.hpp:429
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2680
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:289
void Swap(DenseMatrix &other)
Definition: densemat.cpp:2268
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:1684
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:1099
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:3924
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2615
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:1107
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:157
void SetSubMatrix(const Array< int > &idx, const DenseMatrix &A)
Set (*this)(idx[i],idx[j]) = A(i,j)
Definition: densemat.cpp:1848
void MultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt = A D A^t, where D is diagonal.
Definition: densemat.cpp:2726
void Eigenvalues(DenseMatrix &b, Vector &ev)
Definition: densemat.hpp:262
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:459
virtual void PrintT(std::ostream &out=mfem::out, int width_=4) const
Prints the transpose matrix to stream out.
Definition: densemat.cpp:2228
void CopyCols(const DenseMatrix &A, int col1, int col2)
Copy columns col1 through col2 from A to *this.
Definition: densemat.cpp:1540
int SizeJ() const
Definition: densemat.hpp:1000
virtual MatrixInverse * Inverse() const
Returns a pointer to the inverse matrix.
Definition: densemat.cpp:444
double a
Definition: lissajous.cpp:41
double * GetColumn(int col)
Definition: densemat.hpp:286
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:2904
virtual ~DenseMatrix()
Destroys dense matrix.
Definition: densemat.cpp:2275
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:1579
void AddMultTranspose(const Vector &x, Vector &y) const
y += A^t x
Definition: densemat.cpp:242
void CopyExceptMN(const DenseMatrix &A, int m, int n)
Copy All rows and columns except m and n from A.
Definition: densemat.cpp:1658
virtual double Det(int m) const
Definition: densemat.hpp:611
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:1351
void Mult(int m, int n, double *X) const
Definition: densemat.cpp:3306
virtual bool Factor(int m, double TOL=0.0)
Compute the LU factorization of the current matrix.
Definition: densemat.cpp:3230
DenseMatrixInverse(bool spd_=false)
Default constructor.
Definition: densemat.hpp:810
static const int ipiv_base
Definition: densemat.hpp:638
Definition: densemat.cpp:1439
virtual void Solve(int m, int n, double *X) const
Definition: densemat.hpp:617
const double * GetData(int k) const
Definition: densemat.hpp:1064
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
Definition: densemat.cpp:1229
void GetRowSums(Vector &l) const
Compute the row sums of the DenseMatrix.
Definition: densemat.cpp:1337
void CalcEigenvalues(double *lambda, double *vec) const
Definition: densemat.cpp:1254
virtual bool Factor(int m, double TOL=0.0)
Definition: densemat.hpp:605
const Memory< double > & GetMemory() const
Definition: densemat.hpp:118
DenseMatrix & RightSingularvectors()
Definition: densemat.hpp:943
void RightSolve(int m, int n, double *X) const
Definition: densemat.cpp:3738
RefCoord t[3]
void Eigenvalues(Vector &ev, DenseMatrix &evect)
Compute eigenvalues and eigenvectors of A x = ev x where A = *this.
Definition: densemat.hpp:253
DenseMatrixEigensystem(DenseMatrix &m)
Definition: densemat.cpp:3971
const double * HostRead() const
Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:1091
int Rank(double tol) const
Definition: densemat.cpp:1214
const double alpha
Definition: ex15.cpp:369
LUFactors(double *data_, int *ipiv_)
Definition: densemat.hpp:647
void AddMult_a(double a, const Vector &x, Vector &y) const
y += a * A.x
Definition: densemat.cpp:260
void RightScaling(const Vector &s)
RightScaling: this = this * diag(s);.
Definition: densemat.cpp:341
DenseTensor(const DenseTensor &other)
Copy constructor: deep copy.
Definition: densemat.hpp:988
Vector data type.
Definition: vector.hpp:60
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
void AddMultTranspose_a(double a, const Vector &x, Vector &y) const
y += a * A^t x
Definition: densemat.cpp:278
ADAt += A D A^t, where D is diagonal.
Definition: densemat.cpp:2698
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:2085
Memory< double > & GetMemory()
Definition: densemat.hpp:1074
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:1553
double * Write(bool on_dev=true)
Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:1095
void GetSubMatrix(const Array< int > &idx, DenseMatrix &A) const
Definition: densemat.cpp:1744
RefCoord s[3]
double * Write(bool on_dev=true)
Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:451
double * Data()
Definition: densemat.hpp:1070
void InvLeftScaling(const Vector &s)
InvLeftScaling this = diag(1./s) * this.
Definition: densemat.cpp:328
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:257
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
std::size_t MemoryUsage() const
Definition: densemat.hpp:1084
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:2375
int Size() const
Get the size of the inverse matrix.
Definition: densemat.hpp:820
Vector & Singularvalues()
Definition: densemat.hpp:940
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:953
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:246
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.cpp:163
void AdjustDofDirection(Array< int > &dofs)
Definition: densemat.cpp:2096
Definition: densemat.cpp:1512
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:3056
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void AddSubMatrix(const Array< int > &idx, const DenseMatrix &A)
(*this)(idx[i],idx[j]) += A(i,j)
Definition: densemat.cpp:1961
void Eigenvalues(DenseMatrix &b, Vector &ev, DenseMatrix &evect)
Compute generalized eigenvalues of A x = ev B x, where A = *this.
Definition: densemat.hpp:266
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:443
virtual void PrintMatlab(std::ostream &out=mfem::out) const
Prints operator in Matlab format.
Definition: densemat.cpp:2209
virtual void Solve(int m, int n, double *X) const
Definition: densemat.cpp:3724
DenseMatrix & operator+=(const double *m)
Definition: densemat.cpp:595
CholeskyFactors(double *data_)
Definition: densemat.hpp:746
void UMult(int m, int n, double *X) const
Definition: densemat.cpp:3646