MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
complex_operator.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_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
25namespace 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*/
69{
70public:
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
121protected:
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
129protected:
132
135
137
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
153 */
155{
156public:
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. */
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{
187protected:
189 bool transa;
191
192 void *Numeric;
193 SuiteSparse_long *AI, *AJ;
194
195 void Init();
196
197public:
198 real_t Control[UMFPACK_CONTROL];
199 mutable real_t 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
257 */
259{
260public:
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
281private:
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
Specialization of the ComplexOperator built from a pair of HypreParMatrices.
ComplexHypreParMatrix(HypreParMatrix *A_Real, HypreParMatrix *A_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
virtual HypreParMatrix & imag()
HypreParMatrix * GetSystemMatrix() const
virtual HypreParMatrix & real()
Real or imaginary part accessor methods.
Mimic the action of a complex operator using two real operators.
virtual Operator & imag()
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.
virtual Type GetType() const
virtual Operator & real()
Real or imaginary part accessor methods.
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 ...
Convention GetConvention() const
@ HERMITIAN
Native convention for Hermitian operators.
@ BLOCK_SYMMETRIC
Alternate convention for damping operators.
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Specialization of the ComplexOperator built from a pair of Sparse Matrices.
virtual Type GetType() const
ComplexSparseMatrix(SparseMatrix *A_Real, SparseMatrix *A_Imag, bool ownReal, bool ownImag, Convention convention=HERMITIAN)
SparseMatrix * GetSystemMatrix() const
virtual SparseMatrix & real()
Real or imaginary part accessor methods.
virtual SparseMatrix & imag()
Interface with UMFPack solver specialized for ComplexSparseMatrix This approach avoids forming a mono...
virtual void Mult(const Vector &b, Vector &x) const
This is solving the system A x = b.
void SetTransposeSolve(bool transa_)
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a ComplexSparseMatrix.
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,...
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...
real_t Control[UMFPACK_CONTROL]
ComplexUMFPackSolver(bool use_long_ints_=false, bool transa_=false)
For larger matrices, if the solver fails, set the parameter use_long_ints_ = true.
void SetPrintLevel(int print_lvl)
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Abstract operator.
Definition operator.hpp:25
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Type
Enumeration defining IDs for some classes derived from Operator.
Definition operator.hpp:284
@ MFEM_ComplexSparseMat
ID for class ComplexSparseMatrix.
Definition operator.hpp:295
@ Complex_Operator
ID for class ComplexOperator.
Definition operator.hpp:294
@ Complex_Hypre_ParCSR
ID for class ComplexHypreParMatrix.
Definition operator.hpp:296
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
Base class for solvers.
Definition operator.hpp:683
Data type sparse matrix.
Definition sparsemat.hpp:51
Vector data type.
Definition vector.hpp:80
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
float real_t
Definition config.hpp:43