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