MFEM  v4.6.0
Finite element discretization library
complex_operator.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_OPERATOR
13 #define MFEM_COMPLEX_OPERATOR
14 
15 #include "operator.hpp"
16 #include "sparsemat.hpp"
17 #ifdef MFEM_USE_MPI
18 #include "hypre.hpp"
19 #endif
20 
21 #ifdef MFEM_USE_SUITESPARSE
22 #include <umfpack.h>
23 #endif
24 
25 namespace mfem
26 {
27 
28 /** @brief Mimic the action of a complex operator using two real operators.
29 
30  This operator requires vectors that are twice the length of its internally
31  stored real operators, Op_Real and Op_Imag. It is assumed that these vectors
32  store the real part of the vector first followed by its imaginary part.
33 
34  ComplexOperator allows one to choose a convention upon construction, which
35  facilitates symmetry.
36 
37  If we let (y_r + i y_i) = (Op_r + i Op_i)(x_r + i x_i) then Matrix-vector
38  products are computed as:
39 
40  1. When Convention::HERMITIAN is used (default)
41  / y_r \ / Op_r -Op_i \ / x_r \
42  | | = | | | |
43  \ y_i / \ Op_i Op_r / \ x_i /
44 
45  2. When Convention::BLOCK_SYMMETRIC is used
46  / y_r \ / Op_r -Op_i \ / x_r \
47  | | = | | | |
48  \-y_i / \-Op_i -Op_r / \ x_i /
49  In other words, Matrix-vector products with Convention::BLOCK_SYMMETRIC
50  compute the complex conjugate of Op*x.
51 
52  Either convention can be used with a given complex operator, however, each
53  of them may be best suited for different classes of problems. For example:
54 
55  1. Convention::HERMITIAN, is well suited for Hermitian operators, i.e.,
56  operators where the real part is symmetric and the imaginary part of the
57  operator is anti-symmetric, hence the name. In such cases the resulting 2
58  x 2 operator will be symmetric.
59 
60  2. Convention::BLOCK_SYMMETRIC, is well suited for operators where both the
61  real and imaginary parts are symmetric. In this case the resulting 2 x 2
62  operator will also be symmetric. Such operators are common when studying
63  damped oscillations, for example.
64 
65  Note: this class cannot be used to represent a general nonlinear complex
66  operator.
67 */
68 class ComplexOperator : public Operator
69 {
70 public:
72  {
73  HERMITIAN, ///< Native convention for Hermitian operators
74  BLOCK_SYMMETRIC ///< Alternate convention for damping operators
75  };
76 
77  /** @brief Constructs complex operator object
78 
79  Note that either @p Op_Real or @p Op_Imag can be NULL, thus eliminating
80  their action (see documentation of the class for more details).
81 
82  In case ownership of the passed operator is transferred to this class
83  through @p ownReal and @p ownImag, the operators will be explicitly
84  destroyed at the end of the life of this object.
85  */
86  ComplexOperator(Operator * Op_Real, Operator * Op_Imag,
87  bool ownReal, bool ownImag,
88  Convention convention = HERMITIAN);
89 
90  virtual ~ComplexOperator();
91 
92  /** @brief Check for existence of real or imaginary part of the operator
93 
94  These methods do not check that the operators are non-zero but only that
95  the operators have been set.
96  */
97  bool hasRealPart() const { return Op_Real_ != NULL; }
98  bool hasImagPart() const { return Op_Imag_ != NULL; }
99 
100  /** @brief Real or imaginary part accessor methods
101 
102  The following accessor methods should only be called if the requested
103  part of the operator is known to exist. This can be checked with
104  hasRealPart() or hasImagPart().
105  */
106  virtual Operator & real();
107  virtual Operator & imag();
108  virtual const Operator & real() const;
109  virtual const Operator & imag() const;
110 
111  virtual void Mult(const Vector &x, Vector &y) const;
112  virtual void MultTranspose(const Vector &x, Vector &y) const;
113 
114  using Operator::Mult;
116 
117  virtual Type GetType() const { return Complex_Operator; }
118 
120 
121 protected:
122  // Let this be hidden from the public interface since the implementation
123  // depends on internal members
124  void Mult(const Vector &x_r, const Vector &x_i,
125  Vector &y_r, Vector &y_i) const;
126  void MultTranspose(const Vector &x_r, const Vector &x_i,
127  Vector &y_r, Vector &y_i) const;
128 
129 protected:
132 
133  bool ownReal_;
134  bool ownImag_;
135 
137 
138  mutable Vector x_r_, x_i_, y_r_, y_i_;
139  mutable Vector *u_, *v_;
140 };
141 
142 
143 /** @brief Specialization of the ComplexOperator built from a pair of Sparse
144  Matrices.
145 
146  The purpose of this specialization is to construct a single SparseMatrix
147  object which is equivalent to the 2x2 block system that the ComplexOperator
148  mimics. The resulting SparseMatrix can then be passed along to solvers which
149  require access to the CSR matrix data such as SuperLU, STRUMPACK, or similar
150  sparse linear solvers.
151 
152  See ComplexOperator documentation above for more information.
153  */
155 {
156 public:
158  bool ownReal, bool ownImag,
159  Convention convention = HERMITIAN)
160  : ComplexOperator(A_Real, A_Imag, ownReal, ownImag, convention)
161  {}
162 
163  virtual SparseMatrix & real();
164  virtual SparseMatrix & imag();
165 
166  virtual const SparseMatrix & real() const;
167  virtual const SparseMatrix & imag() const;
168 
169  /** Combine the blocks making up this complex operator into a single
170  SparseMatrix. The resulting matrix can be passed to solvers which require
171  access to the matrix entries themselves, such as sparse direct solvers,
172  rather than simply the action of the operator. Note that this combined
173  operator requires roughly twice the memory of the block structured
174  operator. */
175  SparseMatrix * GetSystemMatrix() const;
176 
177  virtual Type GetType() const { return MFEM_ComplexSparseMat; }
178 };
179 
180 #ifdef MFEM_USE_SUITESPARSE
181 /** @brief Interface with UMFPack solver specialized for ComplexSparseMatrix
182  This approach avoids forming a monolithic SparseMatrix which leads
183  to increased memory and flops
184  */
186 {
187 protected:
189  bool transa;
191 
192  void *Numeric;
193  SuiteSparse_long *AI, *AJ;
194 
195  void Init();
196 
197 public:
198  double Control[UMFPACK_CONTROL];
199  mutable double Info[UMFPACK_INFO];
200 
201  /** @brief For larger matrices, if the solver fails, set the parameter @a
202  use_long_ints_ = true. */
203  ComplexUMFPackSolver(bool use_long_ints_ = false, bool transa_ = false)
204  : use_long_ints(use_long_ints_), transa(transa_) { Init(); }
205  /** @brief Factorize the given ComplexSparseMatrix using the defaults.
206  For larger matrices, if the solver fails, set the parameter
207  @a use_long_ints_ = true. */
208  ComplexUMFPackSolver(ComplexSparseMatrix &A, bool use_long_ints_ = false,
209  bool transa_ = false)
210  : use_long_ints(use_long_ints_), transa(transa_) { Init(); SetOperator(A); }
211 
212  /** @brief Factorize the given Operator @a op which must be
213  a ComplexSparseMatrix.
214 
215  The factorization uses the parameters set in the #Control data member.
216  @note This method calls SparseMatrix::SortColumnIndices()
217  for real and imag parts of the ComplexSparseMatrix,
218  modifying the matrices if the column indices are not already sorted. */
219  virtual void SetOperator(const Operator &op);
220 
221  // Set the print level field in the #Control data member.
222  void SetPrintLevel(int print_lvl) { Control[UMFPACK_PRL] = print_lvl; }
223 
224  // This determines the action of MultTranspose (see below for details)
225  void SetTransposeSolve(bool transa_) { transa = transa_; }
226 
227  /** @brief This is solving the system A x = b */
228  virtual void Mult(const Vector &b, Vector &x) const;
229 
230  /** @brief
231  This is solving the system:
232  A^H x = b (when transa = false)
233  This is equivalent to solving the transpose block system for the
234  case of Convention = HERMITIAN
235  A^T x = b (when transa = true)
236  This is equivalent to solving the transpose block system for the
237  case of Convention = BLOCK_SYMMETRIC */
238  virtual void MultTranspose(const Vector &b, Vector &x) const;
239 
240  virtual ~ComplexUMFPackSolver();
241 };
242 
243 #endif
244 
245 #ifdef MFEM_USE_MPI
246 
247 /** @brief Specialization of the ComplexOperator built from a pair of
248  HypreParMatrices.
249 
250  The purpose of this specialization is to construct a single HypreParMatrix
251  object which is equivalent to the 2x2 block system that the ComplexOperator
252  mimics. The resulting HypreParMatrix can then be passed along to solvers
253  which require access to the CSR matrix data such as SuperLU, STRUMPACK, or
254  similar sparse linear solvers.
255 
256  See ComplexOperator documentation above for more information.
257  */
259 {
260 public:
262  bool ownReal, bool ownImag,
263  Convention convention = HERMITIAN);
264 
265  virtual HypreParMatrix & real();
266  virtual HypreParMatrix & imag();
267 
268  virtual const HypreParMatrix & real() const;
269  virtual const HypreParMatrix & imag() const;
270 
271  /** Combine the blocks making up this complex operator into a single
272  HypreParMatrix. The resulting matrix can be passed to solvers which
273  require access to the matrix entries themselves, such as sparse direct
274  solvers or Hypre preconditioners, rather than simply the action of the
275  operator. Note that this combined operator requires roughly twice the
276  memory of the block structured operator. */
278 
279  virtual Type GetType() const { return Complex_Hypre_ParCSR; }
280 
281 private:
282  void getColStartStop(const HypreParMatrix * A_r,
283  const HypreParMatrix * A_i,
284  int & num_recv_procs,
285  HYPRE_BigInt *& offd_col_start_stop) const;
286 
287  MPI_Comm comm_;
288  int myid_;
289  int nranks_;
290 };
291 
292 #endif // MFEM_USE_MPI
293 
294 }
295 
296 #endif // MFEM_COMPLEX_OPERATOR
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
double Control[UMFPACK_CONTROL]
virtual Operator & real()
Real or imaginary part accessor methods.
virtual void Mult(const Vector &b, Vector &x) const
This is solving the system A x = b.
void SetPrintLevel(int print_lvl)
virtual Type GetType() const
ComplexSparseMatrix * mat
Interface with UMFPack solver specialized for ComplexSparseMatrix This approach avoids forming a mono...
virtual Type GetType() const
ComplexHypreParMatrix(HypreParMatrix *A_Real, HypreParMatrix *A_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
Mimic the action of a complex operator using two real operators.
virtual void MultTranspose(const Vector &b, Vector &x) const
This is solving the system: A^H x = b (when transa = false) This is equivalent to solving the transpo...
ID for class ComplexHypreParMatrix.
Definition: operator.hpp:296
ComplexSparseMatrix(SparseMatrix *A_Real, SparseMatrix *A_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
bool hasRealPart() const
Check for existence of real or imaginary part of the operator.
Alternate convention for damping operators.
ComplexOperator(Operator *Op_Real, Operator *Op_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
Constructs complex operator object.
HypreParMatrix * GetSystemMatrix() const
Data type sparse matrix.
Definition: sparsemat.hpp:50
double b
Definition: lissajous.cpp:42
SparseMatrix * GetSystemMatrix() const
virtual SparseMatrix & real()
Real or imaginary part accessor methods.
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:283
void SetTransposeSolve(bool transa_)
Specialization of the ComplexOperator built from a pair of Sparse Matrices.
HYPRE_Int HYPRE_BigInt
virtual Operator & imag()
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a ComplexSparseMatrix.
Convention GetConvention() const
Specialization of the ComplexOperator built from a pair of HypreParMatrices.
virtual SparseMatrix & imag()
Native convention for Hermitian operators.
ID for class ComplexSparseMatrix.
Definition: operator.hpp:295
ComplexUMFPackSolver(ComplexSparseMatrix &A, bool use_long_ints_=false, bool transa_=false)
Factorize the given ComplexSparseMatrix using the defaults. For larger matrices, if the solver fails...
Vector data type.
Definition: vector.hpp:58
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 ...
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
virtual Type GetType() const
Base class for solvers.
Definition: operator.hpp:682
virtual HypreParMatrix & imag()
ID for class ComplexOperator.
Definition: operator.hpp:294
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
virtual HypreParMatrix & real()
Real or imaginary part accessor methods.
ComplexUMFPackSolver(bool use_long_ints_=false, bool transa_=false)
For larger matrices, if the solver fails, set the parameter use_long_ints_ = true.