eigen_solver.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Header file for a class that defines interfaces to Eigensolvers
27 
28 // Include guard to prevent multiple inclusions of the header
29 #ifndef OOMPH_EIGEN_SOLVER_HEADER
30 #define OOMPH_EIGEN_SOLVER_HEADER
31 
32 // Include the header generated by autoconfig
33 #ifdef HAVE_CONFIG_H
34 #include <oomph-lib-config.h>
35 #endif
36 
37 #ifdef OOMPH_HAS_MPI
38 #include "mpi.h"
39 #endif
40 
41 #include <complex>
42 #include "Vector.h"
43 #include "complex_matrices.h"
44 
45 namespace oomph
46 {
47  // Forward definition of problem class
48  class Problem;
49 
50  // Forward definition of matrix class
51  class DoubleMatrixBase;
52 
53  // Forward definition of linear solver class
54  class LinearSolver;
55 
56  //=======================================================================
57  /// Base class for all EigenProblem solves. This simply defines standard
58  /// interfaces so that different solvers can be used easily.
59  //=======================================================================
61  {
62  protected:
63  /// Double value that represents the real part of the shift in
64  /// shifted eigensolvers
65  double Sigma_real;
66 
67  public:
68  /// Empty constructor
69  EigenSolver() : Sigma_real(0.0) {}
70 
71  /// Empty copy constructor
73 
74  /// Empty destructor
75  virtual ~EigenSolver() {}
76 
77  /// Solve the real eigenproblem that is assembled by elements in
78  /// a mesh in a Problem object. Note that the assembled matrices include the
79  /// shift and are real. The eigenvalues and eigenvectors are,
80  /// in general, complex, and the eigenvalues can be NaNs or Infs.
81  /// This function is therefore merely provided as a convenience, to be
82  /// used if the user is confident that NaNs don't arise (e.g. in Arnoldi
83  /// based solvers where typically only a small number of (finite)
84  /// eigenvalues are computed), or if the users is happy to deal with NaNs in
85  /// the subsequent post-processing.
86  /// Function is virtual so it can be overloaded for Arnoldi type solvers
87  /// that compute the (finite) eigenvalues directly
88  /// At least n_eval eigenvalues are computed.
89  virtual void solve_eigenproblem(Problem* const& problem_pt,
90  const int& n_eval,
91  Vector<std::complex<double>>& eigenvalue,
92  Vector<DoubleVector>& eigenvector_real,
93  Vector<DoubleVector>& eigenvector_imag,
94  const bool& do_adjoint_problem = false)
95  {
97  Vector<double> beta;
98 
99  // Call the "safe" version
100  solve_eigenproblem(problem_pt,
101  n_eval,
102  alpha,
103  beta,
104  eigenvector_real,
105  eigenvector_imag,
106  do_adjoint_problem);
107 
108  // Now do the brute force conversion, possibly creating NaNs and Infs...
109  unsigned n = alpha.size();
110  eigenvalue.resize(n);
111  for (unsigned i = 0; i < n; i++)
112  {
113  eigenvalue[i] = alpha[i] / beta[i];
114  }
115  }
116 
117  /// Solve the real eigenproblem that is assembled by elements in
118  /// a mesh in a Problem object. Note that the assembled matrices include the
119  /// shift and are real. The eigenvalues and eigenvectors are,
120  /// in general, complex. Eigenvalues may be infinite and are therefore
121  /// returned as
122  /// \f$ \lambda_i = \alpha_i / \beta_i \f$ where \f$ \alpha_i \f$ is complex
123  /// while \f$ \beta_i \f$ is real. The actual eigenvalues may then be
124  /// computed by doing the division, checking for zero betas to avoid NaNs.
125  /// There's a convenience wrapper to this function that simply computes
126  /// these eigenvalues regardless. That version may die in NaN checking is
127  /// enabled (via the fenv.h header and the associated feenable function).
128  /// At least n_eval eigenvalues are computed.
129  virtual void solve_eigenproblem(Problem* const& problem_pt,
130  const int& n_eval,
131  Vector<std::complex<double>>& alpha,
132  Vector<double>& beta,
133  Vector<DoubleVector>& eigenvector_real,
134  Vector<DoubleVector>& eigenvector_imag,
135  const bool& do_adjoint_problem = false) = 0;
136 
137  /// Set the value of the (real) shift
138  void set_shift(const double& shift_value)
139  {
140  Sigma_real = shift_value;
141  }
142 
143  /// Return the value of the (real) shift (const version)
144  const double& get_shift() const
145  {
146  return Sigma_real;
147  }
148  };
149 
150  ///////////////////////////////////////////////////////////////////////
151  ///////////////////////////////////////////////////////////////////////
152  ///////////////////////////////////////////////////////////////////////
153 
154  //=====================================================================
155  /// Class for the LAPACK QZ eigensolver
156  //=====================================================================
157  class LAPACK_QZ : public EigenSolver
158  {
159  public:
160  /// Empty constructor
162 
163  /// Broken copy constructor
164  LAPACK_QZ(const LAPACK_QZ&) = delete;
165 
166  /// Broken assignment operator
167  void operator=(const LAPACK_QZ&) = delete;
168 
169  /// Empty desctructor
170  virtual ~LAPACK_QZ() {}
171 
172  /// Solve the real eigenproblem that is assembled by elements in
173  /// a mesh in a Problem object. Note that the assembled matrices include the
174  /// shift and are real. The eigenvalues and eigenvectors are,
175  /// in general, complex. Eigenvalues may be infinite and are therefore
176  /// returned as
177  /// \f$ \lambda_i = \alpha_i / \beta_i \f$ where \f$ \alpha_i \f$ is complex
178  /// while \f$ \beta_i \f$ is real. The actual eigenvalues may then be
179  /// computed by doing the division, checking for zero betas to avoid NaNs.
180  /// There's a convenience wrapper to this function that simply computes
181  /// these eigenvalues regardless. That version may die in NaN checking is
182  /// enabled (via the fenv.h header and the associated feenable function).
183  /// At least n_eval eigenvalues are computed.
184  void solve_eigenproblem(Problem* const& problem_pt,
185  const int& n_eval,
186  Vector<std::complex<double>>& alpha,
187  Vector<double>& beta,
188  Vector<DoubleVector>& eigenvector_real,
189  Vector<DoubleVector>& eigenvector_imag,
190  const bool& do_adjoint_problem = false);
191 
192  /// Find the eigenvalues of a complex generalised eigenvalue problem
193  /// specified by \f$ Ax = \lambda Mx \f$. Note: the (real) shift
194  /// that's specifiable via the EigenSolver base class is ignored here.
195  /// A warning gets issued if it's set to a nonzero value.
196  void find_eigenvalues(const ComplexMatrixBase& A,
197  const ComplexMatrixBase& M,
198  Vector<std::complex<double>>& eigenvalue,
199  Vector<Vector<std::complex<double>>>& eigenvector);
200 
201 
202  /// Access to tolerance for checking complex conjugateness of eigenvalues
203  /// (const version)
205  {
207  }
208 
209 
210  /// Access to tolerance for checking complex conjugateness of eigenvalues
212  {
214  }
215 
216  private:
217  /// Helper function called from "raw" lapack
218  /// code
219  void solve_eigenproblem_helper(Problem* const& problem_pt,
220  const int& n_eval,
221  Vector<std::complex<double>>& alpha,
222  Vector<double>& beta,
223  Vector<DoubleVector>& eigenvector);
224 
225 
226  /// Provide diagonstic for DGGEV error return
227  void DGGEV_error(const int& info, const int& n)
228  {
229  // Throw an error
230  std::ostringstream error_stream;
231  error_stream << "Failure in LAPACK_DGGEV(...).\n"
232  << "info = " << info << std::endl;
233  error_stream
234  << "Diagnostics below are from \n\n"
235  << "http://www.netlib.org/lapack/explore-html/d9/d8e/"
236  "group__double_g_eeigen_ga4f59e87e670a755b41cbdd7e97f36bea.html"
237  << std::endl;
238  if (info < 0)
239  {
240  error_stream << -info << "-th input arg had an illegal value\n";
241  }
242  else if (info <= n)
243  {
244  error_stream << "The QZ iteration failed. No eigenvectors have been\n"
245  << "calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)\n"
246  << "should be correct for j=INFO+1,...,N, where \n"
247  << "info = " << info << " and N = " << n << std::endl;
248  }
249  else if (info == (n + 1))
250  {
251  error_stream << "QZ iteration failed in DHGEQZ.\n";
252  }
253  else if (info == (n + 2))
254  {
255  error_stream << "error return from DTGEVC.\n";
256  }
257  error_stream
258  << "Aborting here; if you know how to proceed then\n"
259  << "then implement ability to catch this error and continue\n";
260 
261  throw OomphLibError(
262  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
263  }
264 
265  /// Provide diagonstic for ZGGEV error return
266  void ZGGEV_error(const int& info, const int& n)
267  {
268  // Throw an error
269  std::ostringstream error_stream;
270  error_stream << "Failure in LAPACK_ZGGEV(...).\n"
271  << "info = " << info << std::endl;
272  error_stream << "Diagnostics below are from \n\n"
273  << "http://www.netlib.org/lapack/explore-html/"
274  << "db/d55/group__complex16_g_eeigen_ga79fcce20c"
275  << "617429ccf985e6f123a6171.html" << std::endl;
276  if (info < 0)
277  {
278  error_stream << -info << "-th input arg had an illegal value\n";
279  }
280  else if (info <= n)
281  {
282  error_stream << "The QZ iteration failed. No eigenvectors have been\n"
283  << "calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)\n"
284  << "should be correct for j=INFO+1,...,N, where \n"
285  << "info = " << info << " and N = " << n << std::endl;
286  }
287  else if (info == (n + 1))
288  {
289  error_stream << "QZ iteration failed in ZHGEQZ.\n";
290  }
291  else if (info == (n + 2))
292  {
293  error_stream << "error return from ZTGEVC.\n";
294  }
295  error_stream
296  << "Aborting here; if you know how to proceed then\n"
297  << "then implement ability to catch this error and continue\n";
298 
299  throw OomphLibError(
300  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
301  }
302 
303  /// Tolerance for checking complex conjugateness of eigenvalues
305  };
306 
307 } // namespace oomph
308 
309 #endif
cstr elem_len * i
Definition: cfortran.h:603
Abstract base class for matrices of complex doubles – adds abstract interfaces for solving,...
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
Base class for all EigenProblem solves. This simply defines standard interfaces so that different sol...
Definition: eigen_solver.h:61
virtual ~EigenSolver()
Empty destructor.
Definition: eigen_solver.h:75
EigenSolver(const EigenSolver &)
Empty copy constructor.
Definition: eigen_solver.h:72
EigenSolver()
Empty constructor.
Definition: eigen_solver.h:69
void set_shift(const double &shift_value)
Set the value of the (real) shift.
Definition: eigen_solver.h:138
virtual void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double >> &alpha, Vector< double > &beta, Vector< DoubleVector > &eigenvector_real, Vector< DoubleVector > &eigenvector_imag, const bool &do_adjoint_problem=false)=0
Solve the real eigenproblem that is assembled by elements in a mesh in a Problem object....
double Sigma_real
Double value that represents the real part of the shift in shifted eigensolvers.
Definition: eigen_solver.h:65
virtual void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double >> &eigenvalue, Vector< DoubleVector > &eigenvector_real, Vector< DoubleVector > &eigenvector_imag, const bool &do_adjoint_problem=false)
Solve the real eigenproblem that is assembled by elements in a mesh in a Problem object....
Definition: eigen_solver.h:89
const double & get_shift() const
Return the value of the (real) shift (const version)
Definition: eigen_solver.h:144
Class for the LAPACK QZ eigensolver.
Definition: eigen_solver.h:158
void ZGGEV_error(const int &info, const int &n)
Provide diagonstic for ZGGEV error return.
Definition: eigen_solver.h:266
void solve_eigenproblem_helper(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double >> &alpha, Vector< double > &beta, Vector< DoubleVector > &eigenvector)
Helper function called from "raw" lapack code.
Definition: eigen_solver.cc:65
void find_eigenvalues(const ComplexMatrixBase &A, const ComplexMatrixBase &M, Vector< std::complex< double >> &eigenvalue, Vector< Vector< std::complex< double >>> &eigenvector)
Find the eigenvalues of a complex generalised eigenvalue problem specified by . Note: the (real) shif...
double tolerance_for_ccness_check() const
Access to tolerance for checking complex conjugateness of eigenvalues (const version)
Definition: eigen_solver.h:204
void operator=(const LAPACK_QZ &)=delete
Broken assignment operator.
void solve_eigenproblem(Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double >> &alpha, Vector< double > &beta, Vector< DoubleVector > &eigenvector_real, Vector< DoubleVector > &eigenvector_imag, const bool &do_adjoint_problem=false)
Solve the real eigenproblem that is assembled by elements in a mesh in a Problem object....
void DGGEV_error(const int &info, const int &n)
Provide diagonstic for DGGEV error return.
Definition: eigen_solver.h:227
double Tolerance_for_ccness_check
Tolerance for checking complex conjugateness of eigenvalues.
Definition: eigen_solver.h:304
virtual ~LAPACK_QZ()
Empty desctructor.
Definition: eigen_solver.h:170
LAPACK_QZ()
Empty constructor.
Definition: eigen_solver.h:161
LAPACK_QZ(const LAPACK_QZ &)=delete
Broken copy constructor.
double & tolerance_for_ccness_check()
Access to tolerance for checking complex conjugateness of eigenvalues.
Definition: eigen_solver.h:211
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:154
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...