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