MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
complex_operator.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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  virtual Type GetType() const { return Complex_Operator; }
115 
117 
118 protected:
119  // Let this be hidden from the public interface since the implementation
120  // depends on internal members
121  void Mult(const Vector &x_r, const Vector &x_i,
122  Vector &y_r, Vector &y_i) const;
123  void MultTranspose(const Vector &x_r, const Vector &x_i,
124  Vector &y_r, Vector &y_i) const;
125 
126 protected:
129 
130  bool ownReal_;
131  bool ownImag_;
132 
134 
135  mutable Vector x_r_, x_i_, y_r_, y_i_;
136  mutable Vector *u_, *v_;
137 };
138 
139 
140 /** @brief Specialization of the ComplexOperator built from a pair of Sparse
141  Matrices.
142 
143  The purpose of this specialization is to construct a single SparseMatrix
144  object which is equivalent to the 2x2 block system that the ComplexOperator
145  mimics. The resulting SparseMatrix can then be passed along to solvers which
146  require access to the CSR matrix data such as SuperLU, STRUMPACK, or similar
147  sparse linear solvers.
148 
149  See ComplexOperator documentation above for more information.
150  */
152 {
153 public:
155  bool ownReal, bool ownImag,
156  Convention convention = HERMITIAN)
157  : ComplexOperator(A_Real, A_Imag, ownReal, ownImag, convention)
158  {}
159 
160  virtual SparseMatrix & real();
161  virtual SparseMatrix & imag();
162 
163  virtual const SparseMatrix & real() const;
164  virtual const SparseMatrix & imag() const;
165 
166  /** Combine the blocks making up this complex operator into a single
167  SparseMatrix. The resulting matrix can be passed to solvers which require
168  access to the matrix entries themselves, such as sparse direct solvers,
169  rather than simply the action of the operator. Note that this combined
170  operator requires roughly twice the memory of the block structured
171  operator. */
172  SparseMatrix * GetSystemMatrix() const;
173 
174  virtual Type GetType() const { return MFEM_ComplexSparseMat; }
175 };
176 
177 #ifdef MFEM_USE_SUITESPARSE
178 /** @brief Interface with UMFPack solver specialized for ComplexSparseMatrix
179  This approach avoids forming a monolithic SparseMatrix which leads
180  to increased memory and flops
181  */
183 {
184 protected:
186  bool transa;
188 
189  void *Numeric;
190  SuiteSparse_long *AI, *AJ;
191 
192  void Init();
193 
194 public:
195  double Control[UMFPACK_CONTROL];
196  mutable double Info[UMFPACK_INFO];
197 
198  /** @brief For larger matrices, if the solver fails, set the parameter @a
199  _use_long_ints = true. */
200  ComplexUMFPackSolver(bool _use_long_ints = false, bool transa_ = false)
201  : use_long_ints(_use_long_ints), transa(transa_) { Init(); }
202  /** @brief Factorize the given ComplexSparseMatrix using the defaults.
203  For larger matrices, if the solver fails, set the parameter
204  @a _use_long_ints = true. */
205  ComplexUMFPackSolver(ComplexSparseMatrix &A, bool _use_long_ints = false,
206  bool transa_ = false)
207  : use_long_ints(_use_long_ints), transa(transa_) { Init(); SetOperator(A); }
208 
209  /** @brief Factorize the given Operator @a op which must be
210  a ComplexSparseMatrix.
211 
212  The factorization uses the parameters set in the #Control data member.
213  @note This method calls SparseMatrix::SortColumnIndices()
214  for real and imag parts of the ComplexSparseMatrix,
215  modifying the matrices if the column indices are not already sorted. */
216  virtual void SetOperator(const Operator &op);
217 
218  // Set the print level field in the #Control data member.
219  void SetPrintLevel(int print_lvl) { Control[UMFPACK_PRL] = print_lvl; }
220 
221  // This determines the action of MultTranspose (see below for details)
222  void SetTransposeSolve(bool transa_) { transa = transa_; }
223 
224  /** @brief This is solving the system A x = b */
225  virtual void Mult(const Vector &b, Vector &x) const;
226 
227  /** @brief
228  This is solving the system:
229  A^H x = b (when transa = false)
230  This is equivalent to solving the transpose block system for the
231  case of Convention = HERMITIAN
232  A^T x = b (when transa = true)
233  This is equivalent to solving the transpose block system for the
234  case of Convention = BLOCK_SYMMETRIC */
235  virtual void MultTranspose(const Vector &b, Vector &x) const;
236 
237  virtual ~ComplexUMFPackSolver();
238 };
239 
240 #endif
241 
242 #ifdef MFEM_USE_MPI
243 
244 /** @brief Specialization of the ComplexOperator built from a pair of
245  HypreParMatrices.
246 
247  The purpose of this specialization is to construct a single HypreParMatrix
248  object which is equivalent to the 2x2 block system that the ComplexOperator
249  mimics. The resulting HypreParMatrix can then be passed along to solvers
250  which require access to the CSR matrix data such as SuperLU, STRUMPACK, or
251  similar sparse linear solvers.
252 
253  See ComplexOperator documentation above for more information.
254  */
256 {
257 public:
259  bool ownReal, bool ownImag,
260  Convention convention = HERMITIAN);
261 
262  virtual HypreParMatrix & real();
263  virtual HypreParMatrix & imag();
264 
265  virtual const HypreParMatrix & real() const;
266  virtual const HypreParMatrix & imag() const;
267 
268  /** Combine the blocks making up this complex operator into a single
269  HypreParMatrix. The resulting matrix can be passed to solvers which
270  require access to the matrix entries themselves, such as sparse direct
271  solvers or Hypre preconditioners, rather than simply the action of the
272  operator. Note that this combined operator requires roughly twice the
273  memory of the block structured operator. */
275 
276  virtual Type GetType() const { return Complex_Hypre_ParCSR; }
277 
278 private:
279  void getColStartStop(const HypreParMatrix * A_r,
280  const HypreParMatrix * A_i,
281  int & num_recv_procs,
282  HYPRE_Int *& offd_col_start_stop) const;
283 
284  MPI_Comm comm_;
285  int myid_;
286  int nranks_;
287 };
288 
289 #endif // MFEM_USE_MPI
290 
291 }
292 
293 #endif // MFEM_COMPLEX_OPERATOR
double Control[UMFPACK_CONTROL]
virtual Operator & real()
Real or imaginary part accessor methods.
virtual Type GetType() const
void SetPrintLevel(int print_lvl)
ComplexSparseMatrix * mat
Interface with UMFPack solver specialized for ComplexSparseMatrix This approach avoids forming a mono...
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.
ID for class ComplexHypreParMatrix.
Definition: operator.hpp:250
ComplexSparseMatrix(SparseMatrix *A_Real, SparseMatrix *A_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
Alternate convention for damping operators.
ComplexOperator(Operator *Op_Real, Operator *Op_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
Constructs complex operator object.
bool hasRealPart() const
Check for existence of real or imaginary part of the operator.
Data type sparse matrix.
Definition: sparsemat.hpp:46
double b
Definition: lissajous.cpp:42
virtual Type GetType() const
virtual SparseMatrix & real()
Real or imaginary part accessor methods.
SparseMatrix * GetSystemMatrix() const
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:237
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void SetTransposeSolve(bool transa_)
Specialization of the ComplexOperator built from a pair of Sparse Matrices.
virtual Operator & imag()
Convention GetConvention() const
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a ComplexSparseMatrix.
virtual void Mult(const Vector &b, Vector &x) const
This is solving the system A x = b.
Specialization of the ComplexOperator built from a pair of HypreParMatrices.
virtual SparseMatrix & imag()
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...
Native convention for Hermitian 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...
virtual Type GetType() const
ID for class ComplexSparseMatrix.
Definition: operator.hpp:249
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 ...
Vector data type.
Definition: vector.hpp:51
ComplexUMFPackSolver(bool _use_long_ints=false, bool transa_=false)
For larger matrices, if the solver fails, set the parameter _use_long_ints = true.
HypreParMatrix * GetSystemMatrix() const
Base class for solvers.
Definition: operator.hpp:634
virtual HypreParMatrix & imag()
ID for class ComplexOperator.
Definition: operator.hpp:248
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:181
virtual HypreParMatrix & real()
Real or imaginary part accessor methods.