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