MFEM  v4.5.2
Finite element discretization library
complex_densemat.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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_COMPLEX_DENSEMAT
13 #define MFEM_COMPLEX_DENSEMAT
14 
15 #include "complex_operator.hpp"
16 #include <complex>
17 
18 namespace 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
24  See ComplexOperator documentation for more information.
25  */
27 {
28 
29 public:
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. */
44  DenseMatrix * GetSystemMatrix() const;
45 
46  virtual Type GetType() const { return Complex_DenseMat; }
47 
49 
50 };
51 
52 /// Matrix matrix multiplication. A = B * C.
53 ComplexDenseMatrix * 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
57 ComplexDenseMatrix * 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 {
65 protected:
66 
67  // returns a new complex array
68  std::complex<double> * RealToComplex(int m, const double * x_r,
69  const double * x_i) const;
70  // copies the given complex array to real and imag arrays
71  void ComplexToReal(int m, const std::complex<double> * x, double * x_r,
72  double * x_i) const;
73 
74 public:
75 
76  double *data_r = nullptr;
77  double *data_i = nullptr;
78  std::complex<double> * data = nullptr;
79 
81 
82  ComplexFactors(double *data_r_, double *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;
90  SetComplexData(m);
91  }
92 
93  virtual bool Factor(int m, double TOL = 0.0)
94  {
95  mfem_error("ComplexFactors::ComplexFactors(...)");
96  return false;
97  }
98 
99  virtual std::complex<double> Det(int m) const
100  {
101  mfem_error("Factors::Det(...)");
102  return 0.;
103  }
104 
105  virtual void Solve(int m, int n, double *X_r, double * X_i) const
106  {
107  mfem_error("Factors::Solve(...)");
108  }
109 
110  virtual void GetInverseMatrix(int m, double *X_r, double * X_i) const
111  {
112  mfem_error("Factors::GetInverseMatrix(...)");
113  }
114 
115  virtual ~ComplexFactors()
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 {
125 public:
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(double *data_r_,double * 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, double 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<double> 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, double *X_r, double * X_i) const;
160 
161  void Mult(int m, int n, std::complex<double> *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, double *X_r, double *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, double *X_r, double *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, double *X_r, double *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, double *X_r, double *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, double *X_r, double * 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 {
188 public:
189 
190  /** With this constructor, the (public) data should be set
191  explicitly before calling class methods. */
193 
194  ComplexCholeskyFactors(double *data_r_, double * data_i_)
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, double 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<double> 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, double *X_r, double * 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, double *X_r, double *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, double *X_r, double * 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, double *X_r, double *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, double *X_r, double * 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, double *X_r, double *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, double *X_r, double * X_i) const;
241 
242 };
243 
244 } // namespace mfem
245 
246 #endif // MFEM_COMPLEX_DENSEMAT
ComplexDenseMatrix * ComputeInverse()
virtual DenseMatrix & real()
Real or imaginary part accessor methods.
void USolve(int m, int n, double *X_r, double *X_i) const
virtual std::complex< double > Det(int m) const
virtual bool Factor(int m, double TOL=0.0)
Compute the Cholesky factorization of the current matrix.
Mimic the action of a complex operator using two real operators.
virtual void GetInverseMatrix(int m, double *X_r, double *X_i) const
Assuming L.L^H = A factored data of size (m x m), compute X <- A^{-1}.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
virtual void GetInverseMatrix(int m, double *X_r, double *X_i) const
virtual bool Factor(int m, double TOL=0.0)
Compute the LU factorization of the current matrix.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual std::complex< double > Det(int m) const
std::complex< double > * data
virtual void Solve(int m, int n, double *X_r, double *X_i) const
ComplexLUFactors(double *data_r_, double *data_i, int *ipiv_)
void UMult(int m, int n, double *X_r, double *X_i) const
ID for class ComplexDenseMatrix.
Definition: operator.hpp:297
ComplexCholeskyFactors(double *data_r_, double *data_i_)
virtual void Solve(int m, int n, double *X_r, double *X_i) const
DenseMatrix * GetSystemMatrix() const
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
virtual std::complex< double > Det(int m) const
void USolve(int m, int n, double *X_r, double *X_i) const
void RightSolve(int m, int n, double *X_r, double *X_i) const
void LSolve(int m, int n, double *X_r, double *X_i) const
virtual bool Factor(int m, double TOL=0.0)
virtual Type GetType() const
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:283
ComplexDenseMatrix(DenseMatrix *A_Real, DenseMatrix *A_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
void LMult(int m, int n, double *X_r, double *X_i) const
ComplexFactors(double *data_r_, double *data_i_)
void ComplexToReal(int m, const std::complex< double > *x, double *x_r, double *x_i) const
void RightSolve(int m, int n, double *X_r, double *X_i) const
virtual void Solve(int m, int n, double *X_r, double *X_i) const
std::complex< double > * RealToComplex(int m, const double *x_r, const double *x_i) const
static const int ipiv_base
Native convention for Hermitian operators.
void LSolve(int m, int n, double *X_r, double *X_i) const
void Mult(int m, int n, double *X_r, double *X_i) const
virtual DenseMatrix & imag()
virtual void GetInverseMatrix(int m, double *X_r, double *X_i) const
Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1}.
Specialization of the ComplexOperator built from a pair of Dense Matrices. The purpose of this specia...