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