MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
complex_densemat.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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_COMPLEX_DENSEMAT
13#define MFEM_COMPLEX_DENSEMAT
14
15#include "complex_operator.hpp"
16#include <complex>
17
18namespace mfem
19{
20
21/** @brief Specialization of the ComplexOperator built from a pair of
22 Dense Matrices. The purpose of this specialization is to support
23 the inverse of a ComplexDenseMatrix and various MatMat operations
25 */
27{
28
29public:
31 bool ownReal, bool ownImag, Convention convention = HERMITIAN)
32 : ComplexOperator(A_Real, A_Imag, ownReal, ownImag, convention)
33 { }
34
35 virtual DenseMatrix & real();
36 virtual DenseMatrix & imag();
37
38 virtual const DenseMatrix & real() const;
39 virtual const DenseMatrix & imag() const;
40
41 /** Combine the blocks making up this complex operator into a single
42 DenseMatrix. Note that this combined operator requires roughly
43 twice the memory of the block structured operator. */
45
46 virtual Type GetType() const { return Complex_DenseMat; }
47
49
50};
51
52/// Matrix matrix multiplication. A = B * C.
53ComplexDenseMatrix * Mult(const ComplexDenseMatrix &B,
54 const ComplexDenseMatrix &C);
55
56/// Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B
57ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A,
58 const ComplexDenseMatrix &B);
59
60
61/** Abstract class that can compute factorization of external data and perform various
62 operations with the factored data. */
64{
65protected:
66
67 // returns a new complex array
68 std::complex<real_t> * RealToComplex(int m, const real_t * x_r,
69 const real_t * x_i) const;
70 // copies the given complex array to real and imag arrays
71 void ComplexToReal(int m, const std::complex<real_t> * x, real_t * x_r,
72 real_t * x_i) const;
73
74public:
75
76 real_t *data_r = nullptr;
77 real_t *data_i = nullptr;
78 std::complex<real_t> * data = nullptr;
79
81
82 ComplexFactors(real_t *data_r_, real_t *data_i_)
83 : data_r(data_r_), data_i(data_i_) { }
84
85 void SetComplexData(int m);
86
87 void ResetComplexData(int m)
88 {
89 delete [] data; data = nullptr;
91 }
92
93 virtual bool Factor(int m, real_t TOL = 0.0)
94 {
95 mfem_error("ComplexFactors::ComplexFactors(...)");
96 return false;
97 }
98
99 virtual std::complex<real_t> Det(int m) const
100 {
101 mfem_error("Factors::Det(...)");
102 return 0.;
103 }
104
105 virtual void Solve(int m, int n, real_t *X_r, real_t * X_i) const
106 {
107 mfem_error("Factors::Solve(...)");
108 }
109
110 virtual void GetInverseMatrix(int m, real_t *X_r, real_t * X_i) const
111 {
112 mfem_error("Factors::GetInverseMatrix(...)");
113 }
114
116 {
117 delete [] data; data = nullptr;
118 }
119};
120
121/** Class that computes factorization of external data and perform various
122 operations with the factored data. */
124{
125public:
126 int *ipiv;
127#ifdef MFEM_USE_LAPACK
128 static const int ipiv_base = 1;
129#else
130 static const int ipiv_base = 0;
131#endif
132
133 /** With this constructor, the (public) data and ipiv members should be set
134 explicitly before calling class methods. */
136
137 ComplexLUFactors(real_t *data_r_,real_t * data_i, int *ipiv_)
138 : ComplexFactors(data_r_, data_i), ipiv(ipiv_) { }
139 /**
140 * @brief Compute the LU factorization of the current matrix
141 *
142 * Factorize the current matrix of size (m x m) overwriting it with the
143 * LU factors. The factorization is such that L.U = P.A, where A is the
144 * original matrix and P is a permutation matrix represented by ipiv.
145 *
146 * @param [in] m size of the square matrix
147 * @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0.
148 *
149 * @return status set to true if successful, otherwise, false.
150 */
151 virtual bool Factor(int m, real_t TOL = 0.0);
152
153 /** Assuming L.U = P.A factored data of size (m x m), compute |A|
154 from the diagonal values of U and the permutation information. */
155 virtual std::complex<real_t> Det(int m) const;
156
157 /** Assuming L.U = P.A factored data of size (m x m), compute X <- A X,
158 for a matrix X of size (m x n). */
159 void Mult(int m, int n, real_t *X_r, real_t * X_i) const;
160
161 void Mult(int m, int n, std::complex<real_t> *X) const;
162
163 /** Assuming L.U = P.A factored data of size (m x m), compute
164 X <- L^{-1} P X, for a matrix X of size (m x n). */
165 void LSolve(int m, int n, real_t *X_r, real_t *X_i) const;
166
167 /** Assuming L.U = P.A factored data of size (m x m), compute
168 X <- U^{-1} X, for a matrix X of size (m x n). */
169 void USolve(int m, int n, real_t *X_r, real_t *X_i) const;
170
171 /** Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1} X,
172 for a matrix X of size (m x n). */
173 virtual void Solve(int m, int n, real_t *X_r, real_t *X_i) const;
174
175 /** Assuming L.U = P.A factored data of size (m x m), compute X <- X A^{-1},
176 for a matrix X of size (n x m). */
177 void RightSolve(int m, int n, real_t *X_r, real_t *X_i) const;
178
179 /// Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1}.
180 virtual void GetInverseMatrix(int m, real_t *X_r, real_t * X_i) const;
181};
182
183
184/** Class that can compute Cholesky factorizations of external data of an
185 Hermitian positive matrix and perform various operations with the factored data. */
187{
188public:
189
190 /** With this constructor, the (public) data should be set
191 explicitly before calling class methods. */
193
195 : ComplexFactors(data_r_, data_i_) { }
196
197 /**
198 * @brief Compute the Cholesky factorization of the current matrix
199 *
200 * Factorize the current matrix of size (m x m) overwriting it with the
201 * Cholesky factors. The factorization is such that LL^H = A, where A is the
202 * original matrix
203 *
204 * @param [in] m size of the square matrix
205 * @param [in] TOL optional fuzzy comparison tolerance. Defaults to 0.0.
206 *
207 * @return status set to true if successful, otherwise, false.
208 */
209 virtual bool Factor(int m, real_t TOL = 0.0);
210
211 /** Assuming LL^H = A factored data of size (m x m), compute |A|
212 from the diagonal values of L */
213 virtual std::complex<real_t> Det(int m) const;
214
215 /** Assuming L.L^H = A factored data of size (m x m), compute X <- L X,
216 for a matrix X of size (m x n). */
217 void LMult(int m, int n, real_t *X_r, real_t * X_i) const;
218
219 /** Assuming L.L^H = A factored data of size (m x m), compute X <- L^t X,
220 for a matrix X of size (m x n). */
221 void UMult(int m, int n, real_t *X_r, real_t *X_i) const;
222
223 /** Assuming L L^H = A factored data of size (m x m), compute
224 X <- L^{-1} X, for a matrix X of size (m x n). */
225 void LSolve(int m, int n, real_t *X_r, real_t * X_i) const;
226
227 /** Assuming L L^H = A factored data of size (m x m), compute
228 X <- L^{-t} X, for a matrix X of size (m x n). */
229 void USolve(int m, int n, real_t *X_r, real_t *X_i) const;
230
231 /** Assuming L.L^H = A factored data of size (m x m), compute X <- A^{-1} X,
232 for a matrix X of size (m x n). */
233 virtual void Solve(int m, int n, real_t *X_r, real_t * X_i) const;
234
235 /** Assuming L.L^H = A factored data of size (m x m), compute X <- X A^{-1},
236 for a matrix X of size (n x m). */
237 void RightSolve(int m, int n, real_t *X_r, real_t *X_i) const;
238
239 /// Assuming L.L^H = A factored data of size (m x m), compute X <- A^{-1}.
240 virtual void GetInverseMatrix(int m, real_t *X_r, real_t * X_i) const;
241
242};
243
244} // namespace mfem
245
246#endif // MFEM_COMPLEX_DENSEMAT
void UMult(int m, int n, real_t *X_r, real_t *X_i) const
void LSolve(int m, int n, real_t *X_r, real_t *X_i) const
ComplexCholeskyFactors(real_t *data_r_, real_t *data_i_)
virtual void GetInverseMatrix(int m, real_t *X_r, real_t *X_i) const
Assuming L.L^H = A factored data of size (m x m), compute X <- A^{-1}.
virtual void Solve(int m, int n, real_t *X_r, real_t *X_i) const
void USolve(int m, int n, real_t *X_r, real_t *X_i) const
virtual std::complex< real_t > Det(int m) const
void RightSolve(int m, int n, real_t *X_r, real_t *X_i) const
virtual bool Factor(int m, real_t TOL=0.0)
Compute the Cholesky factorization of the current matrix.
void LMult(int m, int n, real_t *X_r, real_t *X_i) const
Specialization of the ComplexOperator built from a pair of Dense Matrices. The purpose of this specia...
virtual Type GetType() const
ComplexDenseMatrix(DenseMatrix *A_Real, DenseMatrix *A_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
DenseMatrix * GetSystemMatrix() const
ComplexDenseMatrix * ComputeInverse()
virtual DenseMatrix & real()
Real or imaginary part accessor methods.
virtual DenseMatrix & imag()
virtual std::complex< real_t > Det(int m) const
virtual void Solve(int m, int n, real_t *X_r, real_t *X_i) const
virtual bool Factor(int m, real_t TOL=0.0)
std::complex< real_t > * data
virtual void GetInverseMatrix(int m, real_t *X_r, real_t *X_i) const
void ComplexToReal(int m, const std::complex< real_t > *x, real_t *x_r, real_t *x_i) const
ComplexFactors(real_t *data_r_, real_t *data_i_)
std::complex< real_t > * RealToComplex(int m, const real_t *x_r, const real_t *x_i) const
void Mult(int m, int n, std::complex< real_t > *X) const
virtual void GetInverseMatrix(int m, real_t *X_r, real_t *X_i) const
Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1}.
ComplexLUFactors(real_t *data_r_, real_t *data_i, int *ipiv_)
void LSolve(int m, int n, real_t *X_r, real_t *X_i) const
void USolve(int m, int n, real_t *X_r, real_t *X_i) const
void RightSolve(int m, int n, real_t *X_r, real_t *X_i) const
virtual bool Factor(int m, real_t TOL=0.0)
Compute the LU factorization of the current matrix.
virtual void Solve(int m, int n, real_t *X_r, real_t *X_i) const
virtual std::complex< real_t > Det(int m) const
void Mult(int m, int n, real_t *X_r, real_t *X_i) const
Mimic the action of a complex operator using two real operators.
@ HERMITIAN
Native convention for Hermitian operators.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Type
Enumeration defining IDs for some classes derived from Operator.
Definition operator.hpp:284
@ Complex_DenseMat
ID for class ComplexDenseMatrix.
Definition operator.hpp:297
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:476
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