MFEM  v4.1.0 Finite element discretization library
densemat.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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  /// Change the data array and the size of the DenseMatrix.
62  /** The DenseMatrix does not assume ownership of the data array, i.e. it will
63  not delete the data array @a d. This method should not be used with
64  DenseMatrix that owns its current data array. */
65  void UseExternalData(double *d, int h, int w)
66  {
67  data.Wrap(d, h*w, false);
68  height = h; width = w;
69  }
70
71  /// Change the data array and the size of the DenseMatrix.
72  /** The DenseMatrix does not assume ownership of the data array, i.e. it will
73  not delete the new array @a d. This method will delete the current data
74  array, if owned. */
75  void Reset(double *d, int h, int w)
76  { if (OwnsData()) { data.Delete(); } UseExternalData(d, h, w); }
77
78  /** Clear the data array and the dimensions of the DenseMatrix. This method
79  should not be used with DenseMatrix that owns its current data array. */
80  void ClearExternalData() { data.Reset(); height = width = 0; }
81
82  /// Delete the matrix data array (if owned) and reset the matrix state.
83  void Clear()
84  { if (OwnsData()) { data.Delete(); } ClearExternalData(); }
85
86  /// For backward compatibility define Size to be synonym of Width()
87  int Size() const { return Width(); }
88
89  /// Change the size of the DenseMatrix to s x s.
90  void SetSize(int s) { SetSize(s, s); }
91
92  /// Change the size of the DenseMatrix to h x w.
93  void SetSize(int h, int w);
94
95  /// Returns the matrix data array.
96  inline double *Data() const
97  { return const_cast<double*>((const double*)data);}
98
99  /// Returns the matrix data array.
100  inline double *GetData() const { return Data(); }
101
102  Memory<double> &GetMemory() { return data; }
103  const Memory<double> &GetMemory() const { return data; }
104
105  /// Return the DenseMatrix data (host pointer) ownership flag.
106  inline bool OwnsData() const { return data.OwnsHostPtr(); }
107
108  /// Returns reference to a_{ij}.
109  inline double &operator()(int i, int j);
110
111  /// Returns constant reference to a_{ij}.
112  inline const double &operator()(int i, int j) const;
113
114  /// Matrix inner product: tr(A^t B)
115  double operator*(const DenseMatrix &m) const;
116
117  /// Trace of a square matrix
118  double Trace() const;
119
120  /// Returns reference to a_{ij}.
121  virtual double &Elem(int i, int j);
122
123  /// Returns constant reference to a_{ij}.
124  virtual const double &Elem(int i, int j) const;
125
126  /// Matrix vector multiplication.
127  void Mult(const double *x, double *y) const;
128
129  /// Matrix vector multiplication.
130  virtual void Mult(const Vector &x, Vector &y) const;
131
132  /// Multiply a vector with the transpose matrix.
133  void MultTranspose(const double *x, double *y) const;
134
135  /// Multiply a vector with the transpose matrix.
136  virtual void MultTranspose(const Vector &x, Vector &y) const;
137
138  /// y += A.x
139  void AddMult(const Vector &x, Vector &y) const;
140
141  /// y += A^t x
142  void AddMultTranspose(const Vector &x, Vector &y) const;
143
144  /// y += a * A.x
145  void AddMult_a(double a, const Vector &x, Vector &y) const;
146
147  /// y += a * A^t x
148  void AddMultTranspose_a(double a, const Vector &x, Vector &y) const;
149
150  /// Compute y^t A x
151  double InnerProduct(const double *x, const double *y) const;
152
153  /// LeftScaling this = diag(s) * this
154  void LeftScaling(const Vector & s);
155  /// InvLeftScaling this = diag(1./s) * this
156  void InvLeftScaling(const Vector & s);
157  /// RightScaling: this = this * diag(s);
158  void RightScaling(const Vector & s);
159  /// InvRightScaling: this = this * diag(1./s);
160  void InvRightScaling(const Vector & s);
161  /// SymmetricScaling this = diag(sqrt(s)) * this * diag(sqrt(s))
162  void SymmetricScaling(const Vector & s);
163  /// InvSymmetricScaling this = diag(sqrt(1./s)) * this * diag(sqrt(1./s))
164  void InvSymmetricScaling(const Vector & s);
165
166  /// Compute y^t A x
167  double InnerProduct(const Vector &x, const Vector &y) const
168  { return InnerProduct((const double *)x, (const double *)y); }
169
170  /// Returns a pointer to the inverse matrix.
171  virtual MatrixInverse *Inverse() const;
172
173  /// Replaces the current matrix with its inverse
174  void Invert();
175
176  /// Replaces the current matrix with its square root inverse
177  void SquareRootInverse();
178
179  /// Calculates the determinant of the matrix
180  /// (optimized for 2x2, 3x3, and 4x4 matrices)
181  double Det() const;
182
183  double Weight() const;
184
185  /** @brief Set the matrix to alpha * A, assuming that A has the same
186  dimensions as the matrix and uses column-major layout. */
187  void Set(double alpha, const double *A);
188  /// Set the matrix to alpha * A.
189  void Set(double alpha, const DenseMatrix &A)
190  {
191  SetSize(A.Height(), A.Width());
192  Set(alpha, A.GetData());
193  }
194
195  /// Adds the matrix A multiplied by the number c to the matrix
196  void Add(const double c, const DenseMatrix &A);
197
198  /// Sets the matrix elements equal to constant c
199  DenseMatrix &operator=(double c);
200
201  /// Copy the matrix entries from the given array
202  DenseMatrix &operator=(const double *d);
203
204  /// Sets the matrix size and elements equal to those of m
205  DenseMatrix &operator=(const DenseMatrix &m);
206
207  DenseMatrix &operator+=(const double *m);
209
211
212  DenseMatrix &operator*=(double c);
213
214  /// (*this) = -(*this)
215  void Neg();
216
217  /// Take the 2-norm of the columns of A and store in v
218  void Norm2(double *v) const;
219
220  /// Compute the norm ||A|| = max_{ij} |A_{ij}|
221  double MaxMaxNorm() const;
222
223  /// Compute the Frobenius norm of the matrix
224  double FNorm() const { double s, n2; FNorm(s, n2); return s*sqrt(n2); }
225
226  /// Compute the square of the Frobenius norm of the matrix
227  double FNorm2() const { double s, n2; FNorm(s, n2); return s*s*n2; }
228
229  /// Compute eigenvalues of A x = ev x where A = *this
230  void Eigenvalues(Vector &ev)
231  { Eigensystem(ev); }
232
233  /// Compute eigenvalues and eigenvectors of A x = ev x where A = *this
234  void Eigenvalues(Vector &ev, DenseMatrix &evect)
235  { Eigensystem(ev, &evect); }
236
237  /// Compute eigenvalues and eigenvectors of A x = ev x where A = *this
238  void Eigensystem(Vector &ev, DenseMatrix &evect)
239  { Eigensystem(ev, &evect); }
240
241  /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
242  where A = *this */
244  { Eigensystem(b, ev); }
245
246  /// Compute generalized eigenvalues of A x = ev B x, where A = *this
248  { Eigensystem(b, ev, &evect); }
249
250  /** Compute generalized eigenvalues and eigenvectors of A x = ev B x,
251  where A = *this */
253  { Eigensystem(b, ev, &evect); }
254
255  void SingularValues(Vector &sv) const;
256  int Rank(double tol) const;
257
258  /// Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
259  double CalcSingularvalue(const int i) const;
260
261  /** Return the eigenvalues (in increasing order) and eigenvectors of a
262  2x2 or 3x3 symmetric matrix. */
263  void CalcEigenvalues(double *lambda, double *vec) const;
264
265  void GetRow(int r, Vector &row) const;
266  void GetColumn(int c, Vector &col) const;
267  double *GetColumn(int col) { return data + col*height; }
268  const double *GetColumn(int col) const { return data + col*height; }
269
270  void GetColumnReference(int c, Vector &col)
271  { col.SetDataAndSize(data + c * height, height); }
272
273  void SetRow(int r, const double* row);
274  void SetRow(int r, const Vector &row);
275
276  void SetCol(int c, const double* col);
277  void SetCol(int c, const Vector &col);
278
279
280  /// Set all entries of a row to the specified value.
281  void SetRow(int row, double value);
282  /// Set all entries of a column to the specified value.
283  void SetCol(int col, double value);
284
285  /// Returns the diagonal of the matrix
286  void GetDiag(Vector &d) const;
287  /// Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|
288  void Getl1Diag(Vector &l) const;
289  /// Compute the row sums of the DenseMatrix
290  void GetRowSums(Vector &l) const;
291
292  /// Creates n x n diagonal matrix with diagonal elements c
293  void Diag(double c, int n);
294  /// Creates n x n diagonal matrix with diagonal given by diag
295  void Diag(double *diag, int n);
296
297  /// (*this) = (*this)^t
298  void Transpose();
299  /// (*this) = A^t
300  void Transpose(const DenseMatrix &A);
301  /// (*this) = 1/2 ((*this) + (*this)^t)
302  void Symmetrize();
303
304  void Lump();
305
306  /** Given a DShape matrix (from a scalar FE), stored in *this, returns the
307  CurlShape matrix. If *this is a N by D matrix, then curl is a D*N by
308  D*(D-1)/2 matrix. The size of curl must be set outside. The dimension D
309  can be either 2 or 3. */
311  /** Given a DShape matrix (from a scalar FE), stored in *this,
312  returns the DivShape vector. If *this is a N by dim matrix,
313  then div is a dim*N vector. The size of div must be set
314  outside. */
316
317  /// Copy rows row1 through row2 from A to *this
318  void CopyRows(const DenseMatrix &A, int row1, int row2);
319  /// Copy columns col1 through col2 from A to *this
320  void CopyCols(const DenseMatrix &A, int col1, int col2);
321  /// Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this
322  void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco);
323  /// Copy matrix A to the location in *this at row_offset, col_offset
324  void CopyMN(const DenseMatrix &A, int row_offset, int col_offset);
325  /// Copy matrix A^t to the location in *this at row_offset, col_offset
326  void CopyMNt(const DenseMatrix &A, int row_offset, int col_offset);
327  /** Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this at
328  row_offset, col_offset */
329  void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco,
330  int row_offset, int col_offset);
331  /// Copy c on the diagonal of size n to *this at row_offset, col_offset
332  void CopyMNDiag(double c, int n, int row_offset, int col_offset);
333  /// Copy diag on the diagonal of size n to *this at row_offset, col_offset
334  void CopyMNDiag(double *diag, int n, int row_offset, int col_offset);
335  /// Copy All rows and columns except m and n from A
336  void CopyExceptMN(const DenseMatrix &A, int m, int n);
337
338  /// Perform (ro+i,co+j)+=A(i,j) for 0<=i<A.Height, 0<=j<A.Width
339  void AddMatrix(DenseMatrix &A, int ro, int co);
340  /// Perform (ro+i,co+j)+=a*A(i,j) for 0<=i<A.Height, 0<=j<A.Width
341  void AddMatrix(double a, const DenseMatrix &A, int ro, int co);
342
343  /// Add the matrix 'data' to the Vector 'v' at the given 'offset'
344  void AddToVector(int offset, Vector &v) const;
345  /// Get the matrix 'data' from the Vector 'v' at the given 'offset'
346  void GetFromVector(int offset, const Vector &v);
347  /** If (dofs[i] < 0 and dofs[j] >= 0) or (dofs[i] >= 0 and dofs[j] < 0)
348  then (*this)(i,j) = -(*this)(i,j). */
350
351  /// Replace small entries, abs(a_ij) <= eps, with zero.
352  void Threshold(double eps);
353
354  /** Count the number of entries in the matrix for which isfinite
355  is false, i.e. the entry is a NaN or +/-Inf. */
356  int CheckFinite() const { return mfem::CheckFinite(data, height*width); }
357
358  /// Prints matrix to stream out.
359  virtual void Print(std::ostream &out = mfem::out, int width_ = 4) const;
360  virtual void PrintMatlab(std::ostream &out = mfem::out) const;
361  /// Prints the transpose matrix to stream out.
362  virtual void PrintT(std::ostream &out = mfem::out, int width_ = 4) const;
363
364  /// Invert and print the numerical conditioning of the inversion.
365  void TestInversion();
366
367  long MemoryUsage() const { return data.Capacity() * sizeof(double); }
368
369  /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
370  const double *Read(bool on_dev = true) const
371  { return mfem::Read(data, Height()*Width(), on_dev); }
372
373  /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
375  { return mfem::Read(data, Height()*Width(), false); }
376
377  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
378  double *Write(bool on_dev = true)
379  { return mfem::Write(data, Height()*Width(), on_dev); }
380
381  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
382  double *HostWrite()
383  { return mfem::Write(data, Height()*Width(), false); }
384
385  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
386  double *ReadWrite(bool on_dev = true)
387  { return mfem::ReadWrite(data, Height()*Width(), on_dev); }
388
389  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
391  { return mfem::ReadWrite(data, Height()*Width(), false); }
392
393  /// Destroys dense matrix.
394  virtual ~DenseMatrix();
395 };
396
397 /// C = A + alpha*B
398 void Add(const DenseMatrix &A, const DenseMatrix &B,
399  double alpha, DenseMatrix &C);
400
401 /// C = alpha*A + beta*B
402 void Add(double alpha, const double *A,
403  double beta, const double *B, DenseMatrix &C);
404
405 /// C = alpha*A + beta*B
406 void Add(double alpha, const DenseMatrix &A,
407  double beta, const DenseMatrix &B, DenseMatrix &C);
408
409 /// @brief Solves the dense linear system, A * X = B for X
410 ///
411 /// @param [in,out] A the square matrix for the linear system
412 /// @param [in,out] X the rhs vector, B, on input, the solution, X, on output.
413 /// @param [in] TOL optional fuzzy comparison tolerance. Defaults to 1e-9.
414 ///
415 /// @return status set to true if successful, otherwise, false.
416 ///
417 /// @note This routine may replace the contents of the input Matrix, A, with the
418 /// corresponding LU factorization of the matrix. Matrices of size 1x1 and
419 /// 2x2 are handled explicitly.
420 ///
421 /// @pre A.IsSquare() == true
422 /// @pre X != nullptr
423 bool LinearSolve(DenseMatrix& A, double* X, double TOL = 1.e-9);
424
425 /// Matrix matrix multiplication. A = B * C.
426 void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);
427
428 /// Matrix matrix multiplication. A += B * C.
429 void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a);
430
431 /// Matrix matrix multiplication. A += alpha * B * C.
432 void AddMult_a(double alpha, const DenseMatrix &b, const DenseMatrix &c,
433  DenseMatrix &a);
434
435 /** Calculate the adjugate of a matrix (for NxN matrices, N=1,2,3) or the matrix
436  adj(A^t.A).A^t for rectangular matrices (2x1, 3x1, or 3x2). This operation
437  is well defined even when the matrix is not full rank. */
439
440 /// Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
442
443 /** Calculate the inverse of a matrix (for NxN matrices, N=1,2,3) or the
444  left inverse (A^t.A)^{-1}.A^t (for 2x1, 3x1, or 3x2 matrices) */
445 void CalcInverse(const DenseMatrix &a, DenseMatrix &inva);
446
447 /// Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
448 void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva);
449
450 /** For a given Nx(N-1) (N=2,3) matrix J, compute a vector n such that
451  n_k = (-1)^{k+1} det(J_k), k=1,..,N, where J_k is the matrix J with the
452  k-th row removed. Note: J^t.n = 0, det([n|J])=|n|^2=det(J^t.J). */
453 void CalcOrtho(const DenseMatrix &J, Vector &n);
454
455 /// Calculate the matrix A.At
456 void MultAAt(const DenseMatrix &a, DenseMatrix &aat);
457
458 /// ADAt = A D A^t, where D is diagonal
460
461 /// ADAt += A D A^t, where D is diagonal
463
464 /// Multiply a matrix A with the transpose of a matrix B: A*Bt
465 void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);
466
467 /// ADBt = A D B^t, where D is diagonal
468 void MultADBt(const DenseMatrix &A, const Vector &D,
469  const DenseMatrix &B, DenseMatrix &ADBt);
470
471 /// ABt += A * B^t
472 void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt);
473
474 /// ADBt = A D B^t, where D is diagonal
476  const DenseMatrix &B, DenseMatrix &ADBt);
477
478 /// ABt += a * A * B^t
479 void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B,
480  DenseMatrix &ABt);
481
482 /// Multiply the transpose of a matrix A with a matrix B: At*B
483 void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB);
484
485 /// AAt += a * A * A^t
486 void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);
487
488 /// AAt = a * A * A^t
489 void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt);
490
491 /// Make a matrix from a vector V.Vt
492 void MultVVt(const Vector &v, DenseMatrix &vvt);
493
494 void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);
495
496 /// VWt += v w^t
497 void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt);
498
499 /// VVt += v v^t
500 void AddMultVVt(const Vector &v, DenseMatrix &VWt);
501
502 /// VWt += a * v w^t
503 void AddMult_a_VWt(const double a, const Vector &v, const Vector &w,
504  DenseMatrix &VWt);
505
506 /// VVt += a * v v^t
507 void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt);
508
509
510 /** Class that can compute LU factorization of external data and perform various
511  operations with the factored data. */
513 {
514 public:
515  double *data;
516  int *ipiv;
517 #ifdef MFEM_USE_LAPACK
518  static const int ipiv_base = 1;
519 #else
520  static const int ipiv_base = 0;
521 #endif
522
523  /** With this constructor, the (public) data and ipiv members should be set
524  explicitly before calling class methods. */
525  LUFactors() { }
526
527  LUFactors(double *data_, int *ipiv_) : data(data_), ipiv(ipiv_) { }
528
529  /**
530  * @brief Compute the LU factorization of the current matrix
531  *
532  * Factorize the current matrix of size (m x m) overwriting it with the
533  * LU factors. The factorization is such that L.U = P.A, where A is the
534  * original matrix and P is a permutation matrix represented by ipiv.
535  *
536  * @param [in] m size of the square matrix
537  * @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0.
538  *
539  * @return status set to true if successful, otherwise, false.
540  */
541  bool Factor(int m, double TOL = 0.0);
542
543  /** Assuming L.U = P.A factored data of size (m x m), compute |A|
544  from the diagonal values of U and the permutation information. */
545  double Det(int m) const;
546
547  /** Assuming L.U = P.A factored data of size (m x m), compute X <- A X,
548  for a matrix X of size (m x n). */
549  void Mult(int m, int n, double *X) const;
550
551  /** Assuming L.U = P.A factored data of size (m x m), compute
552  X <- L^{-1} P X, for a matrix X of size (m x n). */
553  void LSolve(int m, int n, double *X) const;
554
555  /** Assuming L.U = P.A factored data of size (m x m), compute
556  X <- U^{-1} X, for a matrix X of size (m x n). */
557  void USolve(int m, int n, double *X) const;
558
559  /** Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1} X,
560  for a matrix X of size (m x n). */
561  void Solve(int m, int n, double *X) const;
562
563  /** Assuming L.U = P.A factored data of size (m x m), compute X <- X A^{-1},
564  for a matrix X of size (n x m). */
565  void RightSolve(int m, int n, double *X) const;
566
567  /// Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1}.
568  void GetInverseMatrix(int m, double *X) const;
569
570  /** Given an (n x m) matrix A21, compute X2 <- X2 - A21 X1, for matrices X1,
571  and X2 of size (m x r) and (n x r), respectively. */
572  static void SubMult(int m, int n, int r, const double *A21,
573  const double *X1, double *X2);
574
575  /** Assuming P.A = L.U factored data of size (m x m), compute the 2x2 block
576  decomposition:
577  | P 0 | | A A12 | = | L 0 | | U U12 |
578  | 0 I | | A21 A22 | | L21 I | | 0 S22 |
579  where A12, A21, and A22 are matrices of size (m x n), (n x m), and
580  (n x n), respectively. The blocks are overwritten as follows:
581  A12 <- U12 = L^{-1} P A12
582  A21 <- L21 = A21 U^{-1}
583  A22 <- S22 = A22 - L21 U12.
584  The block S22 is the Schur complement. */
585  void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const;
586
587  /** Given BlockFactor()'d data, perform the forward block solve for the
588  linear system:
589  | A A12 | | X1 | = | B1 |
590  | A21 A22 | | X2 | | B2 |
591  written in the factored form:
592  | L 0 | | U U12 | | X1 | = | P 0 | | B1 |
593  | L21 I | | 0 S22 | | X2 | | 0 I | | B2 |.
594  The resulting blocks Y1, Y2 solve the system:
595  | L 0 | | Y1 | = | P 0 | | B1 |
596  | L21 I | | Y2 | | 0 I | | B2 |
597  The blocks are overwritten as follows:
598  B1 <- Y1 = L^{-1} P B1
599  B2 <- Y2 = B2 - L21 Y1 = B2 - A21 A^{-1} B1
600  The blocks B1/Y1 and B2/Y2 are of size (m x r) and (n x r), respectively.
601  The Schur complement system is given by: S22 X2 = Y2. */
602  void BlockForwSolve(int m, int n, int r, const double *L21,
603  double *B1, double *B2) const;
604
605  /** Given BlockFactor()'d data, perform the backward block solve in
606  | U U12 | | X1 | = | Y1 |
607  | 0 S22 | | X2 | | Y2 |.
608  The input is the solution block X2 and the block Y1 resulting from
609  BlockForwSolve(). The result block X1 overwrites input block Y1:
610  Y1 <- X1 = U^{-1} (Y1 - U12 X2). */
611  void BlockBackSolve(int m, int n, int r, const double *U12,
612  const double *X2, double *Y1) const;
613 };
614
615
616 /** Data type for inverse of square dense matrix.
617  Stores LU factors */
619 {
620 private:
621  const DenseMatrix *a;
622  LUFactors lu;
623
624 public:
625  /// Default constructor.
626  DenseMatrixInverse() : a(NULL), lu(NULL, NULL) { }
627
628  /** Creates square dense matrix. Computes factorization of mat
629  and stores LU factors. */
630  DenseMatrixInverse(const DenseMatrix &mat);
631
632  /// Same as above but does not factorize the matrix.
633  DenseMatrixInverse(const DenseMatrix *mat);
634
635  /// Get the size of the inverse matrix
636  int Size() const { return Width(); }
637
638  /// Factor the current DenseMatrix, *a
639  void Factor();
640
641  /// Factor a new DenseMatrix of the same size
642  void Factor(const DenseMatrix &mat);
643
644  virtual void SetOperator(const Operator &op);
645
646  /// Matrix vector multiplication with the inverse of dense matrix.
647  virtual void Mult(const Vector &x, Vector &y) const;
648
649  /// Multiply the inverse matrix by another matrix: X = A^{-1} B.
650  void Mult(const DenseMatrix &B, DenseMatrix &X) const;
651
652  /// Multiply the inverse matrix by another matrix: X <- A^{-1} X.
653  void Mult(DenseMatrix &X) const { lu.Solve(width, X.Width(), X.Data()); }
654
655  /// Compute and return the inverse matrix in Ainv.
656  void GetInverseMatrix(DenseMatrix &Ainv) const;
657
658  /// Compute the determinant of the original DenseMatrix using the LU factors.
659  double Det() const { return lu.Det(width); }
660
661  /// Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
662  void TestInversion();
663
664  /// Destroys dense inverse matrix.
665  virtual ~DenseMatrixInverse();
666 };
667
668
670 {
671  DenseMatrix &mat;
672  Vector EVal;
673  DenseMatrix EVect;
674  Vector ev;
675  int n;
676
677 #ifdef MFEM_USE_LAPACK
678  double *work;
679  char jobz, uplo;
680  int lwork, info;
681 #endif
682
683 public:
684
687  void Eval();
688  Vector &Eigenvalues() { return EVal; }
689  DenseMatrix &Eigenvectors() { return EVect; }
690  double Eigenvalue(int i) { return EVal(i); }
691  const Vector &Eigenvector(int i)
692  {
693  ev.SetData(EVect.Data() + i * EVect.Height());
694  return ev;
695  }
697 };
698
699
701 {
702  Vector sv;
703  int m, n;
704
705 #ifdef MFEM_USE_LAPACK
706  double *work;
707  char jobu, jobvt;
708  int lwork, info;
709 #endif
710
711  void Init();
712 public:
713
715  DenseMatrixSVD(int h, int w);
716  void Eval(DenseMatrix &M);
717  Vector &Singularvalues() { return sv; }
718  double Singularvalue(int i) { return sv(i); }
719  ~DenseMatrixSVD();
720 };
721
722 class Table;
723
724 /// Rank 3 tensor (array of matrices)
726 {
727 private:
728  DenseMatrix Mk;
729  Memory<double> tdata;
730  int nk;
731
732 public:
734  {
735  nk = 0;
736  tdata.Reset();
737  }
738
739  DenseTensor(int i, int j, int k)
740  : Mk(NULL, i, j)
741  {
742  nk = k;
743  tdata.New(i*j*k);
744  }
745
746  /// Copy constructor: deep copy
747  DenseTensor(const DenseTensor &other)
748  : Mk(NULL, other.Mk.height, other.Mk.width), nk(other.nk)
749  {
750  const int size = Mk.Height()*Mk.Width()*nk;
751  if (size > 0)
752  {
753  tdata.New(size, other.tdata.GetMemoryType());
754  tdata.CopyFrom(other.tdata, size);
755  }
756  else
757  {
758  tdata.Reset();
759  }
760  }
761
762  int SizeI() const { return Mk.Height(); }
763  int SizeJ() const { return Mk.Width(); }
764  int SizeK() const { return nk; }
765
766  int TotalSize() const { return SizeI()*SizeJ()*SizeK(); }
767
768  void SetSize(int i, int j, int k)
769  {
770  const MemoryType mt = tdata.GetMemoryType();
771  tdata.Delete();
772  Mk.UseExternalData(NULL, i, j);
773  nk = k;
774  tdata.New(i*j*k, mt);
775  }
776
777  void UseExternalData(double *ext_data, int i, int j, int k)
778  {
779  tdata.Delete();
780  Mk.UseExternalData(NULL, i, j);
781  nk = k;
782  tdata.Wrap(ext_data, i*j*k, false);
783  }
784
785  /// Sets the tensor elements equal to constant c
786  DenseTensor &operator=(double c);
787
789  {
790  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
791  Mk.data = Memory<double>(GetData(k), SizeI()*SizeJ(), false);
792  return Mk;
793  }
794  const DenseMatrix &operator()(int k) const
795  { return const_cast<DenseTensor&>(*this)(k); }
796
797  double &operator()(int i, int j, int k)
798  {
799  MFEM_ASSERT_INDEX_IN_RANGE(i, 0, SizeI());
800  MFEM_ASSERT_INDEX_IN_RANGE(j, 0, SizeJ());
801  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
802  return tdata[i+SizeI()*(j+SizeJ()*k)];
803  }
804
805  const double &operator()(int i, int j, int k) const
806  {
807  MFEM_ASSERT_INDEX_IN_RANGE(i, 0, SizeI());
808  MFEM_ASSERT_INDEX_IN_RANGE(j, 0, SizeJ());
809  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
810  return tdata[i+SizeI()*(j+SizeJ()*k)];
811  }
812
813  double *GetData(int k)
814  {
815  MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());
816  return tdata+k*Mk.Height()*Mk.Width();
817  }
818
819  double *Data() { return tdata; }
820
821  const double *Data() const { return tdata; }
822
823  Memory<double> &GetMemory() { return tdata; }
824  const Memory<double> &GetMemory() const { return tdata; }
825
826  /** Matrix-vector product from unassembled element matrices, assuming both
827  'x' and 'y' use the same elem_dof table. */
828  void AddMult(const Table &elem_dof, const Vector &x, Vector &y) const;
829
830  void Clear()
831  { UseExternalData(NULL, 0, 0, 0); }
832
833  long MemoryUsage() const { return nk*Mk.MemoryUsage(); }
834
835  /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
836  const double *Read(bool on_dev = true) const
837  { return mfem::Read(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
838
839  /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
841  { return mfem::Read(tdata, Mk.Height()*Mk.Width()*nk, false); }
842
843  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
844  double *Write(bool on_dev = true)
845  { return mfem::Write(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
846
847  /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
848  double *HostWrite()
849  { return mfem::Write(tdata, Mk.Height()*Mk.Width()*nk, false); }
850
851  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
852  double *ReadWrite(bool on_dev = true)
853  { return mfem::ReadWrite(tdata, Mk.Height()*Mk.Width()*nk, on_dev); }
854
855  /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
857  { return mfem::ReadWrite(tdata, Mk.Height()*Mk.Width()*nk, false); }
858
859  ~DenseTensor() { tdata.Delete(); }
860 };
861
862
863 // Inline methods
864
865 inline double &DenseMatrix::operator()(int i, int j)
866 {
867  MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
868  return data[i+j*height];
869 }
870
871 inline const double &DenseMatrix::operator()(int i, int j) const
872 {
873  MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
874  return data[i+j*height];
875 }
876
877 } // namespace mfem
878
879 #endif
int Size() const
For backward compatibility define Size to be synonym of Width()
Definition: densemat.hpp:87
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
Definition: densemat.cpp:1409
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:2391
DenseMatrix & operator-=(const DenseMatrix &m)
Definition: densemat.cpp:603
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:745
int CheckFinite(const double *v, const int n)
Definition: vector.hpp:383
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:2776
void Mult(DenseMatrix &X) const
Multiply the inverse matrix by another matrix: X &lt;- A^{-1} X.
Definition: densemat.hpp:653
DenseMatrix & operator*=(double c)
Definition: densemat.cpp:616
void GetDiag(Vector &d) const
Returns the diagonal of the matrix.
Definition: densemat.cpp:1301
void SetCol(int c, const double *col)
Definition: densemat.cpp:1805
void UseExternalData(double *ext_data, int i, int j, int k)
Definition: densemat.hpp:777
DenseTensor & operator=(double c)
Sets the tensor elements equal to constant c.
Definition: densemat.cpp:3496
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:2757
double * HostWrite()
Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:382
Memory< double > & GetMemory()
Definition: densemat.hpp:102
void SetRow(int r, const double *row)
Definition: densemat.cpp:1790
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
Definition: densemat.cpp:356
const DenseMatrix & operator()(int k) const
Definition: densemat.hpp:794
void Eigenvalues(Vector &ev)
Compute eigenvalues of A x = ev x where A = *this.
Definition: densemat.hpp:230
void SingularValues(Vector &sv) const
Definition: densemat.cpp:1169
void Delete()
Delete the owned pointers. The Memory is not reset by this method, i.e. it will, generally, not be Empty() after this call.
DenseMatrix & operator()(int k)
Definition: densemat.hpp:788
double Det() const
Definition: densemat.cpp:449
void Eigensystem(DenseMatrix &b, Vector &ev, DenseMatrix &evect)
Definition: densemat.hpp:252
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:63
int SizeK() const
Definition: densemat.hpp:764
long MemoryUsage() const
Definition: densemat.hpp:833
void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const
Definition: densemat.cpp:3173
void BlockBackSolve(int m, int n, int r, const double *U12, const double *X2, double *Y1) const
Definition: densemat.cpp:3209
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
Definition: densemat.cpp:297
Definition: densemat.cpp:2089
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:836
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:312
void AddMult(const Table &elem_dof, const Vector &x, Vector &y) const
Definition: densemat.cpp:3441
void TestInversion()
Invert and print the numerical conditioning of the inversion.
Definition: densemat.cpp:1905
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void CopyRows(const DenseMatrix &A, int row1, int row2)
Copy rows row1 through row2 from A to *this.
Definition: densemat.cpp:1508
double Det(int m) const
Definition: densemat.cpp:2922
void SetSize(int i, int j, int k)
Definition: densemat.hpp:768
void Eval(DenseMatrix &M)
Definition: densemat.cpp:3410
const double * Data() const
Definition: densemat.hpp:821
Abstract data type for matrix inverse.
Definition: matrix.hpp:69
void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += a * A * B^t.
Definition: densemat.cpp:2589
Definition: densemat.hpp:390
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:3250
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3238
void CopyFrom(const Memory &src, int size)
Copy size entries from src to *this.
void GetInverseMatrix(int m, double *X) const
Assuming L.U = P.A factored data of size (m x m), compute X &lt;- A^{-1}.
Definition: densemat.cpp:3091
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:100
const Memory< double > & GetMemory() const
Definition: densemat.hpp:824
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2303
DenseMatrix & operator=(double c)
Sets the matrix elements equal to constant c.
Definition: densemat.cpp:553
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:2730
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:3156
const Vector & Eigenvector(int i)
Definition: densemat.hpp:691
virtual void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: densemat.cpp:1834
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3278
const double & operator()(int i, int j, int k) const
Definition: densemat.hpp:805
void Set(double alpha, const DenseMatrix &A)
Set the matrix to alpha * A.
Definition: densemat.hpp:189
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1926
double & operator()(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.hpp:865
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:3000
double FNorm() const
Compute the Frobenius norm of the matrix.
Definition: densemat.hpp:224
Definition: densemat.hpp:374
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:2161
double & operator()(int i, int j, int k)
Definition: densemat.hpp:797
void RightSolve(int m, int n, double *X) const
Definition: densemat.cpp:3033
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2057
double operator*(const DenseMatrix &m) const
Matrix inner product: tr(A^t B)
Definition: densemat.cpp:186
double Singularvalue(int i)
Definition: densemat.hpp:718
void Reset(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:75
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
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 AddMult_a_VWt(const double a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
Definition: densemat.cpp:2821
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:1271
void BlockForwSolve(int m, int n, int r, const double *L21, double *B1, double *B2) const
Definition: densemat.cpp:3200
DenseMatrixSVD(DenseMatrix &M)
Definition: densemat.cpp:3376
Abstract data type matrix.
Definition: matrix.hpp:27
const double * GetColumn(int col) const
Definition: densemat.hpp:268
void Norm2(double *v) const
Take the 2-norm of the columns of A and store in v.
Definition: densemat.cpp:790
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:2456
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:635
bool LinearSolve(DenseMatrix &A, double *X, double TOL)
Solves the dense linear system, A * X = B for X
Definition: densemat.cpp:1953
virtual ~DenseMatrixInverse()
Destroys dense inverse matrix.
Definition: densemat.cpp:3301
void LSolve(int m, int n, double *X) const
Definition: densemat.cpp:2975
int TotalSize() const
Definition: densemat.hpp:766
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:659
void AddMultVVt(const Vector &v, DenseMatrix &VVt)
VVt += v v^t.
Definition: densemat.cpp:2797
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:1604
double * GetData(int k)
Definition: densemat.hpp:813
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:2843
void Neg()
(*this) = -(*this)
Definition: densemat.cpp:626
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: densemat.cpp:3271
void Solve(int m, int n, double *X) const
Definition: densemat.cpp:3019
Definition: densemat.hpp:852
void Getl1Diag(Vector &l) const
Returns the l1 norm of the rows of the matrix v_i = sum_j |a_ij|.
Definition: densemat.cpp:1315
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:1725
void SetData(double *d)
Definition: vector.hpp:118
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Definition: device.hpp:295
long MemoryUsage() const
Definition: densemat.hpp:367
DenseMatrix & Eigenvectors()
Definition: densemat.hpp:689
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:1287
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:1820
int SizeI() const
Definition: densemat.hpp:762
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2197
void TestInversion()
Print the numerical conditioning of the inversion: ||A^{-1} A - I||.
Definition: densemat.cpp:3290
double MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
Definition: densemat.cpp:803
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:96
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1376
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
Definition: densemat.cpp:2746
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:2496
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:27
DenseTensor(int i, int j, int k)
Definition: densemat.hpp:739
void ClearExternalData()
Definition: densemat.hpp:80
int CheckFinite() const
Definition: densemat.hpp:356
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2329
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:270
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:1665
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:848
void Clear()
Delete the matrix data array (if owned) and reset the matrix state.
Definition: densemat.hpp:83
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2264
Definition: densemat.hpp:856
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:125
ADAt = A D A^t, where D is diagonal.
Definition: densemat.cpp:2375
void Eigenvalues(DenseMatrix &b, Vector &ev)
Definition: densemat.hpp:243
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
Definition: densemat.hpp:386
virtual void PrintT(std::ostream &out=mfem::out, int width_=4) const
Prints the transpose matrix to stream out.
Definition: densemat.cpp:1879
void CopyCols(const DenseMatrix &A, int col1, int col2)
Copy columns col1 through col2 from A to *this.
Definition: densemat.cpp:1521
int SizeJ() const
Definition: densemat.hpp:763
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:267
bool OwnsData() const
Return the DenseMatrix data (host pointer) ownership flag.
Definition: densemat.hpp:106
ADBt = A D B^t, where D is diagonal.
Definition: densemat.cpp:2553
virtual ~DenseMatrix()
Destroys dense matrix.
Definition: densemat.cpp:1919
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:329
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:1560
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:1639
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:1346
void Mult(int m, int n, double *X) const
Definition: densemat.cpp:2939
bool Factor(int m, double TOL=0.0)
Compute the LU factorization of the current matrix.
Definition: densemat.cpp:2863
static const int ipiv_base
Definition: densemat.hpp:518
Definition: densemat.cpp:1434
DenseMatrixInverse()
Default constructor.
Definition: densemat.hpp:626
double CalcSingularvalue(const int i) const
Return the i-th singular value (decreasing order) of NxN matrix, N=1,2,3.
Definition: densemat.cpp:1224
void GetRowSums(Vector &l) const
Compute the row sums of the DenseMatrix.
Definition: densemat.cpp:1332
void CalcEigenvalues(double *lambda, double *vec) const
Definition: densemat.cpp:1249
const Memory< double > & GetMemory() const
Definition: densemat.hpp:103
void Eigenvalues(Vector &ev, DenseMatrix &evect)
Compute eigenvalues and eigenvectors of A x = ev x where A = *this.
Definition: densemat.hpp:234
DenseMatrixEigensystem(DenseMatrix &m)
Definition: densemat.cpp:3308
Definition: densemat.hpp:840
int Rank(double tol) const
Definition: densemat.cpp:1209
const double alpha
Definition: ex15.cpp:336
LUFactors(double *data_, int *ipiv_)
Definition: densemat.hpp:527
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:747
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:2648
Vector data type.
Definition: vector.hpp:48
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:2347
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:1736
Memory< double > & GetMemory()
Definition: densemat.hpp:823
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:1534
double * Write(bool on_dev=true)
Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:844
double * Write(bool on_dev=true)
Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:378
double * Data()
Definition: densemat.hpp:819
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:65
void Eigensystem(Vector &ev, DenseMatrix &evect)
Compute eigenvalues and eigenvectors of A x = ev x where A = *this.
Definition: densemat.hpp:238
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:90
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:2024
int Size() const
Get the size of the inverse matrix.
Definition: densemat.hpp:636
Vector & Singularvalues()
Definition: densemat.hpp:717
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:725
double InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Definition: densemat.hpp:167
double FNorm2() const
Compute the square of the Frobenius norm of the matrix.
Definition: densemat.hpp:227
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
Definition: densemat.cpp:163
Definition: densemat.cpp:1747
Definition: densemat.cpp:1493
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:2705
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void Eigenvalues(DenseMatrix &b, Vector &ev, DenseMatrix &evect)
Compute generalized eigenvalues of A x = ev B x, where A = *this.
Definition: densemat.hpp:247
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:370
virtual void PrintMatlab(std::ostream &out=mfem::out) const
Definition: densemat.cpp:1860
double * data
Definition: densemat.hpp:515
DenseMatrix & operator+=(const double *m)
Definition: densemat.cpp:586