Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #ifndef LIBMESH_SLEPC_EIGEN_SOLVER_H
21 : #define LIBMESH_SLEPC_EIGEN_SOLVER_H
22 :
23 : #include "libmesh/libmesh_config.h"
24 :
25 : #ifdef LIBMESH_HAVE_SLEPC
26 :
27 : // Local includes
28 : #include "libmesh/eigen_solver.h"
29 : #include "libmesh/slepc_macro.h"
30 :
31 : // SLEPc include files.
32 : EXTERN_C_FOR_SLEPC_BEGIN
33 : # include "libmesh/ignore_warnings.h"
34 : # include <slepceps.h>
35 : # include "libmesh/restore_warnings.h"
36 : EXTERN_C_FOR_SLEPC_END
37 :
38 : namespace libMesh
39 : {
40 : template <typename T> class PetscVector;
41 :
42 : /**
43 : * This class provides an interface to the SLEPc
44 : * eigenvalue solver library from http://slepc.upv.es/.
45 : *
46 : * \author Steffen Peterson
47 : * \date 2005
48 : * \brief EigenSolver implementation based on SLEPc.
49 : */
50 : template <typename T>
51 : class SlepcEigenSolver : public EigenSolver<T>
52 : {
53 :
54 : public:
55 :
56 : /**
57 : * Constructor. Initializes Petsc data structures
58 : */
59 : SlepcEigenSolver(const Parallel::Communicator & comm_in);
60 :
61 :
62 : /**
63 : * Destructor.
64 : */
65 : ~SlepcEigenSolver();
66 :
67 :
68 : /**
69 : * Release all memory and clear data structures.
70 : * clear() is called from the destructor, so it should not throw.
71 : */
72 : virtual void clear() noexcept override;
73 :
74 :
75 : /**
76 : * Initialize data structures if not done so already.
77 : */
78 : virtual void init() override;
79 :
80 :
81 : /**
82 : * This function calls the SLEPc solver to compute
83 : * the eigenpairs of the SparseMatrix matrix_A. \p nev is
84 : * the number of eigenpairs to be computed and
85 : * \p ncv is the number of basis vectors to be
86 : * used in the solution procedure. Return values
87 : * are the number of converged eigen values and the
88 : * number of the iterations carried out by the eigen
89 : * solver.
90 : */
91 : virtual std::pair<unsigned int, unsigned int>
92 : solve_standard (SparseMatrix<T> & matrix_A,
93 : int nev,
94 : int ncv,
95 : const double tol,
96 : const unsigned int m_its) override;
97 :
98 : /**
99 : * Same as above except that matrix_A is a ShellMatrix
100 : * in this case.
101 : */
102 : virtual std::pair<unsigned int, unsigned int>
103 : solve_standard (ShellMatrix<T> & shell_matrix,
104 : int nev,
105 : int ncv,
106 : const double tol,
107 : const unsigned int m_its) override;
108 :
109 : /**
110 : * Same as above except that matrix_A is a ShellMatrix
111 : * in this case.
112 : */
113 : virtual std::pair<unsigned int, unsigned int>
114 : solve_standard (ShellMatrix<T> & shell_matrix,
115 : SparseMatrix<T> & precond,
116 : int nev,
117 : int ncv,
118 : const double tol,
119 : const unsigned int m_its) override;
120 :
121 : /**
122 : * Same as above except that precond is a ShellMatrix
123 : * in this case.
124 : */
125 : virtual std::pair<unsigned int, unsigned int>
126 : solve_standard (ShellMatrix<T> & shell_matrix,
127 : ShellMatrix<T> & precond,
128 : int nev,
129 : int ncv,
130 : const double tol,
131 : const unsigned int m_its) override;
132 :
133 :
134 : /**
135 : * This function calls the SLEPc solver to compute
136 : * the eigenpairs for the generalized eigenproblem
137 : * defined by the matrix_A and matrix_B,
138 : * which are of type SparseMatrix. The argument
139 : * \p nev is the number of eigenpairs to be computed
140 : * and \p ncv is the number of basis vectors to be
141 : * used in the solution procedure. Return values
142 : * are the number of converged eigen values and the
143 : * number of the iterations carried out by the eigen
144 : * solver.
145 : */
146 : virtual std::pair<unsigned int, unsigned int>
147 : solve_generalized(SparseMatrix<T> & matrix_A,
148 : SparseMatrix<T> & matrix_B,
149 : int nev,
150 : int ncv,
151 : const double tol,
152 : const unsigned int m_its) override;
153 :
154 : /**
155 : * Solve generalized eigenproblem when matrix_A is of
156 : * type ShellMatrix, matrix_B is of type SparseMatrix.
157 : */
158 : virtual std::pair<unsigned int, unsigned int>
159 : solve_generalized(ShellMatrix<T> & matrix_A,
160 : SparseMatrix<T> & matrix_B,
161 : int nev,
162 : int ncv,
163 : const double tol,
164 : const unsigned int m_its) override;
165 :
166 : /**
167 : * Solve generalized eigenproblem when matrix_A is of
168 : * type SparseMatrix, matrix_B is of type ShellMatrix.
169 : * When using this function, one should use the
170 : * command line options:
171 : * -st_ksp_type gmres -st_pc_type none
172 : * or
173 : * -st_ksp_type gmres -st_pc_type jacobi
174 : * or similar.
175 : */
176 : virtual std::pair<unsigned int, unsigned int>
177 : solve_generalized(SparseMatrix<T> & matrix_A,
178 : ShellMatrix<T> & matrix_B,
179 : int nev,
180 : int ncv,
181 : const double tol,
182 : const unsigned int m_its) override;
183 :
184 : /**
185 : * Solve generalized eigenproblem when both matrix_A and
186 : * matrix_B are of type ShellMatrix.
187 : * When using this function, one should use the
188 : * command line options:
189 : * -st_ksp_type gmres -st_pc_type none
190 : * or
191 : * -st_ksp_type gmres -st_pc_type jacobi
192 : * or similar.
193 : */
194 : virtual std::pair<unsigned int, unsigned int>
195 : solve_generalized(ShellMatrix<T> & matrix_A,
196 : ShellMatrix<T> & matrix_B,
197 : int nev,
198 : int ncv,
199 : const double tol,
200 : const unsigned int m_its) override;
201 :
202 :
203 : virtual std::pair<unsigned int, unsigned int>
204 : solve_generalized(ShellMatrix<T> & matrix_A,
205 : ShellMatrix<T> & matrix_B,
206 : SparseMatrix<T> & precond,
207 : int nev,
208 : int ncv,
209 : const double tol,
210 : const unsigned int m_its) override;
211 :
212 : virtual std::pair<unsigned int, unsigned int>
213 : solve_generalized(ShellMatrix<T> & matrix_A,
214 : ShellMatrix<T> & matrix_B,
215 : ShellMatrix<T> & precond,
216 : int nev,
217 : int ncv,
218 : const double tol,
219 : const unsigned int m_its) override;
220 :
221 : /**
222 : * \returns The real and imaginary part of the ith eigenvalue and
223 : * copies the respective eigenvector to the solution vector.
224 : *
225 : * \note The eigenpair may be complex even for real-valued matrices.
226 : */
227 : virtual std::pair<Real, Real>
228 : get_eigenpair (dof_id_type i,
229 : NumericVector<T> & solution_in) override;
230 :
231 : /**
232 : * Same as above, but does not copy the eigenvector.
233 : */
234 : virtual std::pair<Real, Real>
235 : get_eigenvalue (dof_id_type i) override;
236 :
237 : /**
238 : * \returns The relative error \f$ ||A x - \lambda x|| / |\lambda x| \f$
239 : * of the ith eigenpair (or the equivalent for a general eigenvalue problem).
240 : */
241 : Real get_relative_error (unsigned int i);
242 :
243 : /**
244 : * Attach a deflation space defined by a single vector.
245 : */
246 : virtual void attach_deflation_space(NumericVector<T> & deflation_vector) override;
247 :
248 : /**
249 : * Use \p initial_space_in as the initial guess.
250 : */
251 : virtual void
252 : set_initial_space(NumericVector<T> & initial_space_in) override;
253 :
254 : /**
255 : * \returns The raw SLEPc \p EPS pointer.
256 : */
257 0 : EPS eps() { this->init(); return _eps; }
258 :
259 : /**
260 : * Print the eigenvalues and associated error
261 : */
262 : void print_eigenvalues() const;
263 :
264 : private:
265 :
266 : /**
267 : * Helper function that actually performs the standard eigensolve.
268 : */
269 : std::pair<unsigned int, unsigned int> _solve_standard_helper (Mat mat,
270 : Mat precond,
271 : int nev,
272 : int ncv,
273 : const double tol,
274 : const unsigned int m_its);
275 :
276 : /**
277 : * Helper function that actually performs the generalized eigensolve.
278 : */
279 : std::pair<unsigned int, unsigned int> _solve_generalized_helper (Mat mat_A,
280 : Mat mat_B,
281 : Mat precond,
282 : int nev,
283 : int ncv,
284 : const double tol,
285 : const unsigned int m_its);
286 :
287 : /**
288 : * Helper function that actually performs either eigensolve.
289 : */
290 : std::pair<unsigned int, unsigned int> _solve_helper (Mat precond,
291 : int nev,
292 : int ncv,
293 : const double tol,
294 : const unsigned int m_its);
295 :
296 : /**
297 : * Tells Slepc to use the user-specified solver stored in
298 : * \p _eigen_solver_type
299 : */
300 : void set_slepc_solver_type ();
301 :
302 : /**
303 : * Tells Slepc to deal with the type of problem stored in
304 : * \p _eigen_problem_type
305 : */
306 : void set_slepc_problem_type ();
307 :
308 : /**
309 : * Tells Slepc to compute the spectrum at the position
310 : * stored in \p _position_of_spectrum
311 : */
312 : void set_slepc_position_of_spectrum();
313 :
314 : /**
315 : * Internal function if shell matrix mode is used, this just
316 : * calls the shell matrix's matrix multiplication function.
317 : * See PetscLinearSolver for a similar implementation.
318 : */
319 : static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest);
320 :
321 : /**
322 : * Internal function if shell matrix mode is used, this just
323 : * calls the shell matrix's get_diagonal function.
324 : * Required in order to use Jacobi preconditioning.
325 : */
326 : static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest);
327 :
328 : /**
329 : * Eigenproblem solver context
330 : */
331 : EPS _eps;
332 :
333 : /**
334 : * A vector used for initial space. The vector will be used as the basis for EPS.
335 : */
336 : PetscVector<T>* _initial_space;
337 : };
338 :
339 : } // namespace libMesh
340 :
341 :
342 : #endif // #ifdef LIBMESH_HAVE_SLEPC
343 : #endif // LIBMESH_SLEPC_EIGEN_SOLVER_H
|
