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