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 : #include "libmesh/libmesh_common.h"
21 :
22 : #if defined(LIBMESH_HAVE_SLEPC) && defined(LIBMESH_HAVE_PETSC)
23 :
24 : // Local Includes
25 : #include "libmesh/libmesh_logging.h"
26 : #include "libmesh/petsc_matrix_base.h"
27 : #include "libmesh/petsc_vector.h"
28 : #include "libmesh/slepc_eigen_solver.h"
29 : #include "libmesh/shell_matrix.h"
30 : #include "libmesh/enum_to_string.h"
31 : #include "libmesh/solver_configuration.h"
32 : #include "libmesh/enum_eigen_solver_type.h"
33 : #include "libmesh/petsc_shell_matrix.h"
34 :
35 : // PETSc 3.15 includes a non-release SLEPc 3.14.2 that has already
36 : // deprecated STPrecondSetMatForPC but that hasn't upgraded its
37 : // version number to let us switch to the non-deprecated version. If
38 : // we're in that situation we need to ignore the deprecated warning.
39 : #if SLEPC_VERSION_LESS_THAN(3,15,0) && !SLEPC_VERSION_LESS_THAN(3,14,2)
40 : #include "libmesh/ignore_warnings.h"
41 : #endif
42 :
43 : namespace libMesh
44 : {
45 :
46 :
47 :
48 : template <typename T>
49 348 : SlepcEigenSolver<T>::SlepcEigenSolver (const Parallel::Communicator & comm_in) :
50 : EigenSolver<T>(comm_in),
51 348 : _initial_space(nullptr)
52 : {
53 348 : this->_eigen_solver_type = ARNOLDI;
54 348 : this->_eigen_problem_type = NHEP;
55 348 : }
56 :
57 :
58 :
59 : template <typename T>
60 696 : SlepcEigenSolver<T>::~SlepcEigenSolver ()
61 : {
62 348 : this->SlepcEigenSolver::clear ();
63 696 : }
64 :
65 :
66 :
67 : template <typename T>
68 707 : void SlepcEigenSolver<T>::clear () noexcept
69 : {
70 707 : if (this->initialized())
71 : {
72 359 : this->_is_initialized = false;
73 :
74 359 : PetscErrorCode ierr = LibMeshEPSDestroy(&_eps);
75 : if (ierr)
76 : libmesh_warning("Warning: EPSDestroy returned a non-zero error code which we ignored.");
77 :
78 : // SLEPc default eigenproblem solver
79 359 : this->_eigen_solver_type = KRYLOVSCHUR;
80 : }
81 707 : }
82 :
83 :
84 :
85 : template <typename T>
86 359 : void SlepcEigenSolver<T>::init ()
87 : {
88 : // Initialize the data structures if not done so already.
89 359 : if (!this->initialized())
90 : {
91 359 : this->_is_initialized = true;
92 :
93 : // Create the eigenproblem solver context
94 359 : LibmeshPetscCall(EPSCreate (this->comm().get(), &_eps));
95 :
96 : // Set user-specified solver
97 359 : set_slepc_solver_type();
98 : }
99 359 : }
100 :
101 :
102 :
103 : template <typename T>
104 : std::pair<unsigned int, unsigned int>
105 70 : SlepcEigenSolver<T>::solve_standard (SparseMatrix<T> & matrix_A_in,
106 : int nev, // number of requested eigenpairs
107 : int ncv, // number of basis vectors
108 : const double tol, // solver tolerance
109 : const unsigned int m_its) // maximum number of iterations
110 : {
111 4 : LOG_SCOPE("solve_standard()", "SlepcEigenSolver");
112 :
113 70 : this->clear ();
114 :
115 70 : this->init ();
116 :
117 : // Make sure the SparseMatrix passed in is really a PetscMatrix
118 2 : auto * const matrix_A = cast_ptr<PetscMatrixBase<T> *>(&matrix_A_in);
119 :
120 2 : libmesh_error_msg_if(!matrix_A, "Error: input matrix to solve_standard() must be a PetscMatrixBase.");
121 :
122 : // Close the matrix and vectors in case this wasn't already done.
123 70 : if (this->_close_matrix_before_solve)
124 70 : matrix_A->close ();
125 :
126 74 : return _solve_standard_helper(matrix_A->mat(), nullptr, nev, ncv, tol, m_its);
127 : }
128 :
129 :
130 : template <typename T>
131 : std::pair<unsigned int, unsigned int>
132 0 : SlepcEigenSolver<T>::solve_standard (ShellMatrix<T> & shell_matrix,
133 : int nev, // number of requested eigenpairs
134 : int ncv, // number of basis vectors
135 : const double tol, // solver tolerance
136 : const unsigned int m_its) // maximum number of iterations
137 : {
138 0 : this->clear ();
139 :
140 0 : this->init ();
141 :
142 : // Prepare the matrix. Note that the const_cast is only necessary
143 : // because PETSc does not accept a const void *. Inside the member
144 : // function _petsc_shell_matrix() below, the pointer is casted back
145 : // to a const ShellMatrix<T> *.
146 : Mat mat;
147 0 : LibmeshPetscCall(MatCreateShell(this->comm().get(),
148 : shell_matrix.m(), // Specify the number of local rows
149 : shell_matrix.n(), // Specify the number of local columns
150 : PETSC_DETERMINE,
151 : PETSC_DETERMINE,
152 : const_cast<void *>(static_cast<const void *>(&shell_matrix)),
153 : &mat));
154 :
155 0 : LibmeshPetscCall(MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult)));
156 0 : LibmeshPetscCall(MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal)));
157 :
158 0 : return _solve_standard_helper(mat, nullptr, nev, ncv, tol, m_its);
159 : }
160 :
161 : template <typename T>
162 : std::pair<unsigned int, unsigned int>
163 0 : SlepcEigenSolver<T>::solve_standard (ShellMatrix<T> & shell_matrix,
164 : SparseMatrix<T> & precond_in,
165 : int nev, // number of requested eigenpairs
166 : int ncv, // number of basis vectors
167 : const double tol, // solver tolerance
168 : const unsigned int m_its) // maximum number of iterations
169 : {
170 0 : this->clear ();
171 :
172 0 : this->init ();
173 :
174 : // Make sure the SparseMatrix passed in is really a PetscMatrix
175 0 : auto * const precond = cast_ptr<PetscMatrixBase<T> *>(&precond_in);
176 :
177 0 : libmesh_error_msg_if(!precond, "Error: input preconditioning matrix to solve_standard() must be a PetscMatrixBase.");
178 :
179 0 : auto * const matrix = cast_ptr<PetscShellMatrix<T> *> (&shell_matrix);
180 :
181 0 : libmesh_error_msg_if(!matrix, "Error: input operator matrix to solve_standard() must be a PetscShellMatrix.");
182 :
183 0 : return _solve_standard_helper(matrix->mat(), precond->mat(), nev, ncv, tol, m_its);
184 : }
185 :
186 : template <typename T>
187 : std::pair<unsigned int, unsigned int>
188 0 : SlepcEigenSolver<T>::solve_standard (ShellMatrix<T> & shell_matrix,
189 : ShellMatrix<T> & precond_in,
190 : int nev, // number of requested eigenpairs
191 : int ncv, // number of basis vectors
192 : const double tol, // solver tolerance
193 : const unsigned int m_its) // maximum number of iterations
194 : {
195 0 : this->clear ();
196 :
197 0 : this->init ();
198 :
199 : // Make sure the SparseMatrix passed in is really a PetscMatrix
200 0 : auto * const precond = cast_ptr<PetscShellMatrix<T> *>(&precond_in);
201 :
202 0 : libmesh_error_msg_if(!precond, "Error: input preconditioning matrix to solve_standard() must be a PetscShellMatrix.");
203 :
204 0 : auto * const matrix = cast_ptr<PetscShellMatrix<T> *> (&shell_matrix);
205 :
206 0 : libmesh_error_msg_if(!matrix, "Error: input operator matrix to solve_standard() must be a PetscShellMatrix.");
207 :
208 0 : return _solve_standard_helper(matrix->mat(), precond->mat(), nev, ncv, tol, m_its);
209 : }
210 :
211 :
212 :
213 : template <typename T>
214 : std::pair<unsigned int, unsigned int>
215 70 : SlepcEigenSolver<T>::_solve_standard_helper(Mat mat,
216 : Mat precond,
217 : int nev, // number of requested eigenpairs
218 : int ncv, // number of basis vectors
219 : const double tol, // solver tolerance
220 : const unsigned int m_its) // maximum number of iterations
221 : {
222 4 : LOG_SCOPE("solve_standard()", "SlepcEigenSolver");
223 :
224 : // Set operators.
225 70 : LibmeshPetscCall(EPSSetOperators (_eps, mat, LIBMESH_PETSC_NULLPTR));
226 :
227 74 : return this->_solve_helper(precond, nev, ncv, tol, m_its);
228 : }
229 :
230 :
231 :
232 : template <typename T>
233 : std::pair<unsigned int, unsigned int>
234 223 : SlepcEigenSolver<T>::solve_generalized (SparseMatrix<T> & matrix_A_in,
235 : SparseMatrix<T> & matrix_B_in,
236 : int nev, // number of requested eigenpairs
237 : int ncv, // number of basis vectors
238 : const double tol, // solver tolerance
239 : const unsigned int m_its) // maximum number of iterations
240 : {
241 223 : this->clear ();
242 :
243 223 : this->init ();
244 :
245 : // Make sure the data passed in are really of Petsc types
246 6 : auto * const matrix_A = cast_ptr<PetscMatrixBase<T> *>(&matrix_A_in);
247 6 : auto * const matrix_B = cast_ptr<PetscMatrixBase<T> *>(&matrix_B_in);
248 :
249 6 : libmesh_error_msg_if(!matrix_A || !matrix_B,
250 : "Error: inputs to solve_generalized() must be of type PetscMatrixBase.");
251 :
252 : // Close the matrix and vectors in case this wasn't already done.
253 223 : if (this->_close_matrix_before_solve)
254 : {
255 223 : matrix_A->close ();
256 223 : matrix_B->close ();
257 : }
258 :
259 223 : return _solve_generalized_helper (matrix_A->mat(), matrix_B->mat(), nullptr, nev, ncv, tol, m_its);
260 : }
261 :
262 : template <typename T>
263 : std::pair<unsigned int, unsigned int>
264 0 : SlepcEigenSolver<T>::solve_generalized (ShellMatrix<T> & shell_matrix_A,
265 : SparseMatrix<T> & matrix_B_in,
266 : int nev, // number of requested eigenpairs
267 : int ncv, // number of basis vectors
268 : const double tol, // solver tolerance
269 : const unsigned int m_its) // maximum number of iterations
270 : {
271 0 : this->clear();
272 :
273 0 : this->init ();
274 :
275 : // Prepare the matrix. Note that the const_cast is only necessary
276 : // because PETSc does not accept a const void *. Inside the member
277 : // function _petsc_shell_matrix() below, the pointer is casted back
278 : // to a const ShellMatrix<T> *.
279 : Mat mat_A;
280 0 : LibmeshPetscCall(MatCreateShell(this->comm().get(),
281 : shell_matrix_A.m(), // Specify the number of local rows
282 : shell_matrix_A.n(), // Specify the number of local columns
283 : PETSC_DETERMINE,
284 : PETSC_DETERMINE,
285 : const_cast<void *>(static_cast<const void *>(&shell_matrix_A)),
286 : &mat_A));
287 :
288 0 : auto * const matrix_B = cast_ptr<PetscMatrixBase<T> *>(&matrix_B_in);
289 :
290 0 : libmesh_error_msg_if(!matrix_B, "Error: inputs to solve_generalized() must be of type PetscMatrixBase.");
291 :
292 : // Close the matrix and vectors in case this wasn't already done.
293 0 : if (this->_close_matrix_before_solve)
294 0 : matrix_B->close ();
295 :
296 0 : LibmeshPetscCall(MatShellSetOperation(mat_A,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult)));
297 0 : LibmeshPetscCall(MatShellSetOperation(mat_A,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal)));
298 :
299 0 : return _solve_generalized_helper (mat_A, matrix_B->mat(), nullptr, nev, ncv, tol, m_its);
300 : }
301 :
302 : template <typename T>
303 : std::pair<unsigned int, unsigned int>
304 0 : SlepcEigenSolver<T>::solve_generalized (SparseMatrix<T> & matrix_A_in,
305 : ShellMatrix<T> & shell_matrix_B,
306 : int nev, // number of requested eigenpairs
307 : int ncv, // number of basis vectors
308 : const double tol, // solver tolerance
309 : const unsigned int m_its) // maximum number of iterations
310 : {
311 0 : this->clear();
312 :
313 0 : this->init ();
314 :
315 0 : auto * const matrix_A = cast_ptr<PetscMatrixBase<T> *>(&matrix_A_in);
316 :
317 0 : libmesh_error_msg_if(!matrix_A, "Error: inputs to solve_generalized() must be of type PetscMatrixBase.");
318 :
319 : // Close the matrix and vectors in case this wasn't already done.
320 0 : if (this->_close_matrix_before_solve)
321 0 : matrix_A->close ();
322 :
323 : // Prepare the matrix. Note that the const_cast is only necessary
324 : // because PETSc does not accept a const void *. Inside the member
325 : // function _petsc_shell_matrix() below, the pointer is casted back
326 : // to a const ShellMatrix<T> *.
327 : Mat mat_B;
328 0 : LibmeshPetscCall(MatCreateShell(this->comm().get(),
329 : shell_matrix_B.m(), // Specify the number of local rows
330 : shell_matrix_B.n(), // Specify the number of local columns
331 : PETSC_DETERMINE,
332 : PETSC_DETERMINE,
333 : const_cast<void *>(static_cast<const void *>(&shell_matrix_B)),
334 : &mat_B));
335 :
336 0 : LibmeshPetscCall(MatShellSetOperation(mat_B,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult)));
337 0 : LibmeshPetscCall(MatShellSetOperation(mat_B,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal)));
338 :
339 0 : return _solve_generalized_helper (matrix_A->mat(), mat_B, nullptr, nev, ncv, tol, m_its);
340 : }
341 :
342 : template <typename T>
343 : std::pair<unsigned int, unsigned int>
344 0 : SlepcEigenSolver<T>::solve_generalized (ShellMatrix<T> & shell_matrix_A,
345 : ShellMatrix<T> & shell_matrix_B,
346 : int nev, // number of requested eigenpairs
347 : int ncv, // number of basis vectors
348 : const double tol, // solver tolerance
349 : const unsigned int m_its) // maximum number of iterations
350 : {
351 0 : this->clear();
352 :
353 0 : this->init ();
354 :
355 : // Prepare the matrices. Note that the const_casts are only
356 : // necessary because PETSc does not accept a const void *. Inside
357 : // the member function _petsc_shell_matrix() below, the pointer is
358 : // casted back to a const ShellMatrix<T> *.
359 : Mat mat_A;
360 0 : LibmeshPetscCall(MatCreateShell(this->comm().get(),
361 : shell_matrix_A.m(), // Specify the number of local rows
362 : shell_matrix_A.n(), // Specify the number of local columns
363 : PETSC_DETERMINE,
364 : PETSC_DETERMINE,
365 : const_cast<void *>(static_cast<const void *>(&shell_matrix_A)),
366 : &mat_A));
367 :
368 : Mat mat_B;
369 0 : LibmeshPetscCall(MatCreateShell(this->comm().get(),
370 : shell_matrix_B.m(), // Specify the number of local rows
371 : shell_matrix_B.n(), // Specify the number of local columns
372 : PETSC_DETERMINE,
373 : PETSC_DETERMINE,
374 : const_cast<void *>(static_cast<const void *>(&shell_matrix_B)),
375 : &mat_B));
376 :
377 0 : LibmeshPetscCall(MatShellSetOperation(mat_A,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult)));
378 0 : LibmeshPetscCall(MatShellSetOperation(mat_A,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal)));
379 :
380 0 : LibmeshPetscCall(MatShellSetOperation(mat_B,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult)));
381 0 : LibmeshPetscCall(MatShellSetOperation(mat_B,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal)));
382 :
383 0 : return _solve_generalized_helper (mat_A, mat_B, nullptr, nev, ncv, tol, m_its);
384 : }
385 :
386 : template <typename T>
387 : std::pair<unsigned int, unsigned int>
388 66 : SlepcEigenSolver<T>::solve_generalized (ShellMatrix<T> & shell_matrix_A,
389 : ShellMatrix<T> & shell_matrix_B,
390 : SparseMatrix<T> & precond_in,
391 : int nev, // number of requested eigenpairs
392 : int ncv, // number of basis vectors
393 : const double tol, // solver tolerance
394 : const unsigned int m_its) // maximum number of iterations
395 : {
396 66 : this->clear();
397 :
398 66 : this->init ();
399 :
400 : // Make sure the SparseMatrix passed in is really a PetscMatrix
401 0 : auto * const precond = cast_ptr<PetscMatrixBase<T> *>(&precond_in);
402 :
403 0 : libmesh_error_msg_if(!precond, "Error: input preconditioning matrix to solve_generalized() must be of type PetscMatrixBase.");
404 :
405 0 : auto * const matrix_A = cast_ptr<PetscShellMatrix<T> *> (&shell_matrix_A);
406 :
407 0 : libmesh_error_msg_if(!matrix_A, "Error: input operator A to solve_generalized() must be of type PetscShellMatrix.");
408 :
409 0 : auto * const matrix_B = cast_ptr<PetscShellMatrix<T> *> (&shell_matrix_B);
410 :
411 0 : libmesh_error_msg_if(!matrix_B, "Error: input operator B to solve_generalized() must be of type PetscShellMatrix.");
412 :
413 66 : return _solve_generalized_helper (matrix_A->mat(), matrix_B->mat(), precond->mat(), nev, ncv, tol, m_its);
414 : }
415 :
416 : template <typename T>
417 : std::pair<unsigned int, unsigned int>
418 0 : SlepcEigenSolver<T>::solve_generalized (ShellMatrix<T> & shell_matrix_A,
419 : ShellMatrix<T> & shell_matrix_B,
420 : ShellMatrix<T> & precond_in,
421 : int nev, // number of requested eigenpairs
422 : int ncv, // number of basis vectors
423 : const double tol, // solver tolerance
424 : const unsigned int m_its) // maximum number of iterations
425 : {
426 0 : this->clear();
427 :
428 0 : this->init ();
429 :
430 : // Make sure the ShellMatrix passed in is really a PetscShellMatrix
431 0 : auto * const precond = cast_ptr<PetscShellMatrix<T> *>(&precond_in);
432 :
433 0 : libmesh_error_msg_if(!precond, "Error: input preconditioning matrix to solve_generalized() must be of type PetscShellMatrix.");
434 :
435 0 : auto * const matrix_A = cast_ptr<PetscShellMatrix<T> *> (&shell_matrix_A);
436 :
437 0 : libmesh_error_msg_if(!matrix_A, "Error: input operator A to solve_generalized() must be of type PetscShellMatrix.");
438 :
439 0 : auto * const matrix_B = cast_ptr<PetscShellMatrix<T> *> (&shell_matrix_B);
440 :
441 0 : libmesh_error_msg_if(!matrix_B, "Error: input operator B to solve_generalized() must be of type PetscShellMatrix.");
442 :
443 0 : return _solve_generalized_helper (matrix_A->mat(), matrix_B->mat(), precond->mat(), nev, ncv, tol, m_its);
444 : }
445 :
446 : template <typename T>
447 : std::pair<unsigned int, unsigned int>
448 289 : SlepcEigenSolver<T>::_solve_generalized_helper (Mat mat_A,
449 : Mat mat_B,
450 : Mat precond,
451 : int nev, // number of requested eigenpairs
452 : int ncv, // number of basis vectors
453 : const double tol, // solver tolerance
454 : const unsigned int m_its) // maximum number of iterations
455 : {
456 12 : LOG_SCOPE("solve_generalized()", "SlepcEigenSolver");
457 :
458 : // Set operators.
459 289 : LibmeshPetscCall(EPSSetOperators (_eps, mat_A, mat_B));
460 :
461 301 : return this->_solve_helper(precond, nev, ncv, tol, m_its);
462 : }
463 :
464 :
465 :
466 : template <typename T>
467 : std::pair<unsigned int, unsigned int>
468 359 : SlepcEigenSolver<T>::_solve_helper(Mat precond,
469 : int nev, // number of requested eigenpairs
470 : int ncv, // number of basis vectors
471 : const double tol, // solver tolerance
472 : const unsigned int m_its) // maximum number of iterations
473 : {
474 : // converged eigen pairs and number of iterations
475 359 : PetscInt nconv=0;
476 359 : PetscInt its=0;
477 359 : ST st=nullptr;
478 :
479 : //set the problem type and the position of the spectrum
480 359 : set_slepc_problem_type();
481 359 : set_slepc_position_of_spectrum();
482 :
483 : // Set eigenvalues to be computed.
484 : #if SLEPC_VERSION_LESS_THAN(3,0,0)
485 : LibmeshPetscCall(EPSSetDimensions (_eps, nev, ncv));
486 : #else
487 359 : LibmeshPetscCall(EPSSetDimensions (_eps, nev, ncv, PETSC_DECIDE));
488 : #endif
489 : // Set the tolerance and maximum iterations.
490 359 : LibmeshPetscCall(EPSSetTolerances (_eps, tol, m_its));
491 :
492 : // Set runtime options, e.g.,
493 : // -eps_type <type>, -eps_nev <nev>, -eps_ncv <ncv>
494 : // Similar to PETSc, these options will override those specified
495 : // above as long as EPSSetFromOptions() is called _after_ any
496 : // other customization routines.
497 359 : LibmeshPetscCall(EPSSetFromOptions (_eps));
498 :
499 : // Set a preconditioning matrix to ST
500 359 : if (precond) {
501 66 : LibmeshPetscCall(EPSGetST(_eps,&st));
502 : #if SLEPC_VERSION_LESS_THAN(3,15,0)
503 0 : LibmeshPetscCall(STPrecondSetMatForPC(st, precond));
504 : #else
505 66 : LibmeshPetscCall(STSetPreconditionerMat(st, precond));
506 : #endif
507 : }
508 :
509 : // If the SolverConfiguration object is provided, use it to override
510 : // solver options.
511 359 : if (this->_solver_configuration)
512 : {
513 0 : this->_solver_configuration->configure_solver();
514 : }
515 :
516 : // If an initial space is provided, let us attach it to EPS
517 359 : if (_initial_space) {
518 : // Get a handle for the underlying Vec.
519 136 : Vec initial_vector = _initial_space->vec();
520 :
521 136 : LibmeshPetscCall(EPSSetInitialSpace(_eps, 1, &initial_vector));
522 : }
523 :
524 : // Solve the eigenproblem.
525 359 : LibmeshPetscCall(EPSSolve (_eps));
526 :
527 : // Get the number of iterations.
528 359 : LibmeshPetscCall(EPSGetIterationNumber (_eps, &its));
529 :
530 : // Get number of converged eigenpairs.
531 359 : LibmeshPetscCall(EPSGetConverged(_eps,&nconv));
532 :
533 : // return the number of converged eigenpairs
534 : // and the number of iterations
535 367 : return std::make_pair(nconv, its);
536 : }
537 :
538 :
539 : template <typename T>
540 0 : void SlepcEigenSolver<T>::print_eigenvalues() const
541 : {
542 : // converged eigen pairs and number of iterations
543 0 : PetscInt nconv=0;
544 :
545 : // The relative error.
546 : PetscReal error, re, im;
547 :
548 : // Pointer to vectors of the real parts, imaginary parts.
549 : PetscScalar kr, ki;
550 :
551 : // Get number of converged eigenpairs.
552 0 : LibmeshPetscCall(EPSGetConverged(_eps,&nconv));
553 :
554 : // Display eigenvalues and relative errors.
555 0 : LibmeshPetscCall(PetscPrintf(this->comm().get(),
556 : " k ||Ax-kx||/|kx|\n"
557 : " ----------------- -----------------\n" ));
558 :
559 0 : for (PetscInt i=0; i<nconv; i++ )
560 : {
561 0 : LibmeshPetscCall(EPSGetEigenpair(_eps, i, &kr, &ki, LIBMESH_PETSC_NULLPTR,
562 : LIBMESH_PETSC_NULLPTR));
563 :
564 : #if SLEPC_VERSION_LESS_THAN(3,6,0)
565 : LibmeshPetscCall(EPSComputeRelativeError(_eps, i, &error));
566 : #else
567 0 : LibmeshPetscCall(EPSComputeError(_eps, i, EPS_ERROR_RELATIVE, &error));
568 : #endif
569 :
570 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
571 : re = PetscRealPart(kr);
572 : im = PetscImaginaryPart(kr);
573 : #else
574 0 : re = kr;
575 0 : im = ki;
576 : #endif
577 :
578 0 : if (im != .0)
579 0 : LibmeshPetscCall(PetscPrintf(this->comm().get()," %9f%+9f i %12f\n", double(re), double(im), double(error)));
580 : else
581 0 : LibmeshPetscCall(PetscPrintf(this->comm().get()," %12f %12f\n", double(re), double(error)));
582 : }
583 :
584 0 : LibmeshPetscCall(PetscPrintf(this->comm().get(),"\n" ));
585 0 : }
586 :
587 :
588 : template <typename T>
589 359 : void SlepcEigenSolver<T>::set_slepc_solver_type()
590 : {
591 359 : switch (this->_eigen_solver_type)
592 : {
593 0 : case POWER:
594 0 : LibmeshPetscCall(EPSSetType (_eps, EPSPOWER)); return;
595 0 : case SUBSPACE:
596 0 : LibmeshPetscCall(EPSSetType (_eps, EPSSUBSPACE)); return;
597 0 : case LAPACK:
598 0 : LibmeshPetscCall(EPSSetType (_eps, EPSLAPACK)); return;
599 347 : case ARNOLDI:
600 347 : LibmeshPetscCall(EPSSetType (_eps, EPSARNOLDI)); return;
601 0 : case LANCZOS:
602 0 : LibmeshPetscCall(EPSSetType (_eps, EPSLANCZOS)); return;
603 12 : case KRYLOVSCHUR:
604 12 : LibmeshPetscCall(EPSSetType (_eps, EPSKRYLOVSCHUR)); return;
605 : // case ARPACK:
606 : // LibmeshPetscCall(EPSSetType (_eps, (char *) EPSARPACK)); return;
607 :
608 0 : default:
609 0 : libMesh::err << "ERROR: Unsupported SLEPc Eigen Solver: "
610 0 : << Utility::enum_to_string(this->_eigen_solver_type) << std::endl
611 0 : << "Continuing with SLEPc defaults" << std::endl;
612 : }
613 : }
614 :
615 :
616 :
617 :
618 : template <typename T>
619 359 : void SlepcEigenSolver<T>:: set_slepc_problem_type()
620 : {
621 359 : switch (this->_eigen_problem_type)
622 : {
623 70 : case NHEP:
624 70 : LibmeshPetscCall(EPSSetProblemType (_eps, EPS_NHEP)); return;
625 66 : case GNHEP:
626 66 : LibmeshPetscCall(EPSSetProblemType (_eps, EPS_GNHEP)); return;
627 0 : case HEP:
628 0 : LibmeshPetscCall(EPSSetProblemType (_eps, EPS_HEP)); return;
629 223 : case GHEP:
630 223 : LibmeshPetscCall(EPSSetProblemType (_eps, EPS_GHEP)); return;
631 : #if !SLEPC_VERSION_LESS_THAN(3,3,0)
632 : // EPS_GHIEP added in 3.3.0
633 0 : case GHIEP:
634 0 : LibmeshPetscCall(EPSSetProblemType (_eps, EPS_GHIEP)); return;
635 : #endif
636 :
637 0 : default:
638 0 : libMesh::err << "ERROR: Unsupported SLEPc Eigen Problem: "
639 0 : << this->_eigen_problem_type << std::endl
640 0 : << "Continuing with SLEPc defaults" << std::endl;
641 : }
642 : }
643 :
644 :
645 :
646 : template <typename T>
647 359 : void SlepcEigenSolver<T>:: set_slepc_position_of_spectrum()
648 : {
649 359 : switch (this->_position_of_spectrum)
650 : {
651 136 : case LARGEST_MAGNITUDE:
652 : {
653 136 : LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_LARGEST_MAGNITUDE));
654 2 : return;
655 : }
656 0 : case SMALLEST_MAGNITUDE:
657 : {
658 0 : LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_MAGNITUDE));
659 0 : return;
660 : }
661 3 : case LARGEST_REAL:
662 : {
663 3 : LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_LARGEST_REAL));
664 0 : return;
665 : }
666 10 : case SMALLEST_REAL:
667 : {
668 10 : LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_REAL));
669 0 : return;
670 : }
671 0 : case LARGEST_IMAGINARY:
672 : {
673 0 : LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_LARGEST_IMAGINARY));
674 0 : return;
675 : }
676 0 : case SMALLEST_IMAGINARY:
677 : {
678 0 : LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_IMAGINARY));
679 0 : return;
680 : }
681 :
682 : // The EPS_TARGET_XXX enums were added in SLEPc 3.1
683 : #if !SLEPC_VERSION_LESS_THAN(3,1,0)
684 70 : case TARGET_MAGNITUDE:
685 : {
686 70 : LibmeshPetscCall(EPSSetTarget(_eps, PS(this->_target_val)));
687 70 : LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_TARGET_MAGNITUDE));
688 2 : return;
689 : }
690 140 : case TARGET_REAL:
691 : {
692 140 : LibmeshPetscCall(EPSSetTarget(_eps, PS(this->_target_val)));
693 140 : LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_TARGET_REAL));
694 4 : return;
695 : }
696 0 : case TARGET_IMAGINARY:
697 : {
698 0 : LibmeshPetscCall(EPSSetTarget(_eps, PS(this->_target_val)));
699 0 : LibmeshPetscCall(EPSSetWhichEigenpairs (_eps, EPS_TARGET_IMAGINARY));
700 0 : return;
701 : }
702 : #endif
703 :
704 0 : default:
705 0 : libmesh_error_msg("ERROR: Unsupported SLEPc position of spectrum: " << this->_position_of_spectrum);
706 : }
707 : }
708 :
709 :
710 :
711 :
712 :
713 :
714 : template <typename T>
715 919 : std::pair<Real, Real> SlepcEigenSolver<T>::get_eigenpair(dof_id_type i,
716 : NumericVector<T> & solution_in)
717 : {
718 : PetscReal re, im;
719 :
720 : // Make sure the NumericVector passed in is really a PetscVector
721 919 : PetscVector<T> * solution = dynamic_cast<PetscVector<T> *>(&solution_in);
722 :
723 919 : libmesh_error_msg_if(!solution, "Error getting eigenvector: input vector must be a PetscVector.");
724 :
725 : // real and imaginary part of the ith eigenvalue.
726 : PetscScalar kr, ki;
727 :
728 919 : solution->close();
729 :
730 919 : LibmeshPetscCall(EPSGetEigenpair(_eps, i, &kr, &ki, solution->vec(),
731 : LIBMESH_PETSC_NULLPTR));
732 :
733 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
734 : re = PetscRealPart(kr);
735 : im = PetscImaginaryPart(kr);
736 : #else
737 919 : re = kr;
738 919 : im = ki;
739 : #endif
740 :
741 943 : return std::make_pair(re, im);
742 : }
743 :
744 :
745 : template <typename T>
746 0 : std::pair<Real, Real> SlepcEigenSolver<T>::get_eigenvalue(dof_id_type i)
747 : {
748 : PetscReal re, im;
749 :
750 : // real and imaginary part of the ith eigenvalue.
751 : PetscScalar kr, ki;
752 :
753 0 : LibmeshPetscCall(EPSGetEigenvalue(_eps, i, &kr, &ki));
754 :
755 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
756 : re = PetscRealPart(kr);
757 : im = PetscImaginaryPart(kr);
758 : #else
759 0 : re = kr;
760 0 : im = ki;
761 : #endif
762 :
763 0 : return std::make_pair(re, im);
764 : }
765 :
766 :
767 : template <typename T>
768 0 : Real SlepcEigenSolver<T>::get_relative_error(unsigned int i)
769 : {
770 : PetscReal error;
771 :
772 : #if SLEPC_VERSION_LESS_THAN(3,6,0)
773 : LibmeshPetscCall(EPSComputeRelativeError(_eps, i, &error));
774 : #else
775 0 : LibmeshPetscCall(EPSComputeError(_eps, i, EPS_ERROR_RELATIVE, &error));
776 : #endif
777 :
778 0 : return error;
779 : }
780 :
781 :
782 : template <typename T>
783 0 : void SlepcEigenSolver<T>::attach_deflation_space(NumericVector<T> & deflation_vector_in)
784 : {
785 0 : this->init();
786 :
787 : // Make sure the input vector is actually a PetscVector
788 0 : PetscVector<T> * deflation_vector_petsc_vec =
789 0 : dynamic_cast<PetscVector<T> *>(&deflation_vector_in);
790 :
791 0 : libmesh_error_msg_if(!deflation_vector_petsc_vec, "Error attaching deflation space: input vector must be a PetscVector.");
792 :
793 : // Get a handle for the underlying Vec.
794 0 : Vec deflation_vector = deflation_vector_petsc_vec->vec();
795 :
796 : #if SLEPC_VERSION_LESS_THAN(3,1,0)
797 : LibmeshPetscCall(EPSAttachDeflationSpace(_eps, 1, &deflation_vector, PETSC_FALSE));
798 : #else
799 0 : LibmeshPetscCall(EPSSetDeflationSpace(_eps, 1, &deflation_vector));
800 : #endif
801 0 : }
802 :
803 : template <typename T>
804 136 : void SlepcEigenSolver<T>::set_initial_space(NumericVector<T> & initial_space_in)
805 : {
806 : #if SLEPC_VERSION_LESS_THAN(3,1,0)
807 : libmesh_error_msg("SLEPc 3.1 is required to call EigenSolver::set_initial_space()");
808 : #else
809 : // Make sure the input vector (which is still owned by caller) is
810 : // actually a PetscVector
811 136 : _initial_space = cast_ptr<PetscVector<T> *>(&initial_space_in);
812 : #endif
813 136 : }
814 :
815 : template <typename T>
816 0 : PetscErrorCode SlepcEigenSolver<T>::_petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest)
817 : {
818 : PetscFunctionBegin;
819 :
820 : // Get the matrix context.
821 : void * ctx;
822 0 : LibmeshPetscCallQ(MatShellGetContext(mat,&ctx));
823 :
824 : // Get user shell matrix object.
825 0 : const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
826 :
827 : // Make \p NumericVector instances around the vectors.
828 0 : PetscVector<T> arg_global(arg, shell_matrix.comm());
829 0 : PetscVector<T> dest_global(dest, shell_matrix.comm());
830 :
831 : // Call the user function.
832 0 : shell_matrix.vector_mult(dest_global,arg_global);
833 :
834 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
835 0 : }
836 :
837 : template <typename T>
838 0 : PetscErrorCode SlepcEigenSolver<T>::_petsc_shell_matrix_get_diagonal(Mat mat, Vec dest)
839 : {
840 : PetscFunctionBegin;
841 :
842 : // Get the matrix context.
843 : void * ctx;
844 0 : LibmeshPetscCallQ(MatShellGetContext(mat,&ctx));
845 :
846 : // Get user shell matrix object.
847 0 : const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
848 :
849 : // Make \p NumericVector instances around the vector.
850 0 : PetscVector<T> dest_global(dest, shell_matrix.comm());
851 :
852 : // Call the user function.
853 0 : shell_matrix.get_diagonal(dest_global);
854 :
855 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
856 0 : }
857 :
858 : //------------------------------------------------------------------
859 : // Explicit instantiations
860 : template class LIBMESH_EXPORT SlepcEigenSolver<Number>;
861 :
862 : } // namespace libMesh
863 :
864 :
865 : // In case we do unity builds someday
866 : #if SLEPC_VERSION_LESS_THAN(3,15,0) && !SLEPC_VERSION_LESS_THAN(3,14,2)
867 : #include "libmesh/restore_warnings.h"
868 : #endif
869 :
870 :
871 : #endif // #ifdef LIBMESH_HAVE_SLEPC
|
