slepc_eigen_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
34 
35 namespace libMesh
36 {
37 
38 extern "C"
39 {
40  // Older versions of PETSc do not have the different int typedefs.
41  // On 64-bit machines, PetscInt may actually be a long long int.
42  // This change occurred in Petsc-2.2.1.
43 #if PETSC_VERSION_LESS_THAN(2,2,1)
44  typedef int PetscErrorCode;
45  typedef int PetscInt;
46 #endif
47 }
48 
49 
50 /*----------------------- functions ----------------------------------*/
51 template <typename T>
53 {
54  if (this->initialized())
55  {
56  this->_is_initialized = false;
57 
59 
60  ierr = LibMeshEPSDestroy(&_eps);
61  LIBMESH_CHKERRABORT(ierr);
62 
63  // SLEPc default eigenproblem solver
64 #if SLEPC_VERSION_LESS_THAN(2,3,2)
65  this->_eigen_solver_type = ARNOLDI;
66 #else
67  // Krylov-Schur showed up as of Slepc 2.3.2
68  this->_eigen_solver_type = KRYLOVSCHUR;
69 #endif
70  }
71 }
72 
73 
74 
75 template <typename T>
77 {
78 
80 
81  // Initialize the data structures if not done so already.
82  if (!this->initialized())
83  {
84  this->_is_initialized = true;
85 
86  // Create the eigenproblem solver context
87  ierr = EPSCreate (this->comm().get(), &_eps);
88  LIBMESH_CHKERRABORT(ierr);
89 
90  // Set user-specified solver
91  set_slepc_solver_type();
92  }
93 }
94 
95 
96 
97 template <typename T>
98 std::pair<unsigned int, unsigned int>
100  int nev, // number of requested eigenpairs
101  int ncv, // number of basis vectors
102  const double tol, // solver tolerance
103  const unsigned int m_its) // maximum number of iterations
104 {
105 // START_LOG("solve_standard()", "SlepcEigenSolver");
106 
107  this->init ();
108 
109  // Make sure the SparseMatrix passed in is really a PetscMatrix
110  PetscMatrix<T>* matrix_A = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_A_in);
111 
112  // Close the matrix and vectors in case this wasn't already done.
113  matrix_A->close ();
114 
115  // just for debugging, remove this
116 // char mat_file[] = "matA.petsc";
117 // PetscViewer petsc_viewer;
118 // ierr = PetscViewerBinaryOpen(this->comm().get(), mat_file, PETSC_FILE_CREATE, &petsc_viewer);
119 // LIBMESH_CHKERRABORT(ierr);
120 // ierr = MatView(matrix_A->mat(),petsc_viewer);
121 // LIBMESH_CHKERRABORT(ierr);
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 
138 
139  // Prepare the matrix.
140  Mat mat;
141  ierr = MatCreateShell(this->comm().get(),
142  shell_matrix.m(), // Specify the number of local rows
143  shell_matrix.n(), // Specify the number of local columns
144  PETSC_DETERMINE,
145  PETSC_DETERMINE,
146  const_cast<void*>(static_cast<const void*>(&shell_matrix)),
147  &mat);
148 
149  /* Note that the const_cast above is only necessary because PETSc
150  does not accept a const void*. Inside the member function
151  _petsc_shell_matrix() below, the pointer is casted back to a
152  const ShellMatrix<T>*. */
153 
154  LIBMESH_CHKERRABORT(ierr);
155  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
156  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
157  LIBMESH_CHKERRABORT(ierr);
158 
159 
160  return _solve_standard_helper(mat, nev, ncv, tol, m_its);
161 }
162 
163 template <typename T>
164 std::pair<unsigned int, unsigned int>
166  (Mat mat,
167  int nev, // number of requested eigenpairs
168  int ncv, // number of basis vectors
169  const double tol, // solver tolerance
170  const unsigned int m_its) // maximum number of iterations
171 {
172  START_LOG("solve_standard()", "SlepcEigenSolver");
173 
175 
176  // converged eigen pairs and number of iterations
177  PetscInt nconv=0;
178  PetscInt its=0;
179 
180 #ifdef DEBUG
181  // The relative error.
182  PetscReal error, re, im;
183 
184  // Pointer to vectors of the real parts, imaginary parts.
185  PetscScalar kr, ki;
186 #endif
187 
188  // Set operators.
189  ierr = EPSSetOperators (_eps, mat, PETSC_NULL);
190  LIBMESH_CHKERRABORT(ierr);
191 
192  //set the problem type and the position of the spectrum
193  set_slepc_problem_type();
194  set_slepc_position_of_spectrum();
195 
196  // Set eigenvalues to be computed.
197 #if SLEPC_VERSION_LESS_THAN(3,0,0)
198  ierr = EPSSetDimensions (_eps, nev, ncv);
199 #else
200  ierr = EPSSetDimensions (_eps, nev, ncv, PETSC_DECIDE);
201 #endif
202  LIBMESH_CHKERRABORT(ierr);
203  // Set the tolerance and maximum iterations.
204  ierr = EPSSetTolerances (_eps, tol, m_its);
205  LIBMESH_CHKERRABORT(ierr);
206 
207  // Set runtime options, e.g.,
208  // -eps_type <type>, -eps_nev <nev>, -eps_ncv <ncv>
209  // Similar to PETSc, these options will override those specified
210  // above as long as EPSSetFromOptions() is called _after_ any
211  // other customization routines.
212  ierr = EPSSetFromOptions (_eps);
213  LIBMESH_CHKERRABORT(ierr);
214 
215  // Solve the eigenproblem.
216  ierr = EPSSolve (_eps);
217  LIBMESH_CHKERRABORT(ierr);
218 
219  // Get the number of iterations.
220  ierr = EPSGetIterationNumber (_eps, &its);
221  LIBMESH_CHKERRABORT(ierr);
222 
223  // Get number of converged eigenpairs.
224  ierr = EPSGetConverged(_eps,&nconv);
225  LIBMESH_CHKERRABORT(ierr);
226 
227 
228 #ifdef DEBUG
229  // ierr = PetscPrintf(this->comm().get(),
230  // "\n Number of iterations: %d\n"
231  // " Number of converged eigenpairs: %d\n\n", its, nconv);
232 
233  // Display eigenvalues and relative errors.
234  ierr = PetscPrintf(this->comm().get(),
235  " k ||Ax-kx||/|kx|\n"
236  " ----------------- -----------------\n" );
237  LIBMESH_CHKERRABORT(ierr);
238 
239  for(PetscInt i=0; i<nconv; i++ )
240  {
241  ierr = EPSGetEigenpair(_eps, i, &kr, &ki, PETSC_NULL, PETSC_NULL);
242  LIBMESH_CHKERRABORT(ierr);
243 
244  ierr = EPSComputeRelativeError(_eps, i, &error);
245  LIBMESH_CHKERRABORT(ierr);
246 
247 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
248  re = PetscRealPart(kr);
249  im = PetscImaginaryPart(kr);
250 #else
251  re = kr;
252  im = ki;
253 #endif
254 
255  if (im != .0)
256  {
257  ierr = PetscPrintf(this->comm().get()," %9f%+9f i %12f\n", re, im, error);
258  LIBMESH_CHKERRABORT(ierr);
259  }
260  else
261  {
262  ierr = PetscPrintf(this->comm().get()," %12f %12f\n", re, error);
263  LIBMESH_CHKERRABORT(ierr);
264  }
265  }
266 
267  ierr = PetscPrintf(this->comm().get(),"\n" );
268  LIBMESH_CHKERRABORT(ierr);
269 #endif // DEBUG
270 
271 
272  STOP_LOG("solve_standard()", "SlepcEigenSolver");
273 
274  // return the number of converged eigenpairs
275  // and the number of iterations
276  return std::make_pair(nconv, its);
277 
278 }
279 
280 
281 
282 
283 
284 template <typename T>
285 std::pair<unsigned int, unsigned int>
287  SparseMatrix<T> &matrix_B_in,
288  int nev, // number of requested eigenpairs
289  int ncv, // number of basis vectors
290  const double tol, // solver tolerance
291  const unsigned int m_its) // maximum number of iterations
292 {
293  this->init ();
294 
295  // Make sure the data passed in are really of Petsc types
296  PetscMatrix<T>* matrix_A = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_A_in);
297  PetscMatrix<T>* matrix_B = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_B_in);
298 
299  // Close the matrix and vectors in case this wasn't already done.
300  matrix_A->close ();
301  matrix_B->close ();
302 
303 
304  return _solve_generalized_helper (matrix_A->mat(), matrix_B->mat(), nev, ncv, tol, m_its);
305 }
306 
307 template <typename T>
308 std::pair<unsigned int, unsigned int>
310  SparseMatrix<T> &matrix_B_in,
311  int nev, // number of requested eigenpairs
312  int ncv, // number of basis vectors
313  const double tol, // solver tolerance
314  const unsigned int m_its) // maximum number of iterations
315 {
316  this->init ();
317 
319 
320  // Prepare the matrix.
321  Mat mat_A;
322  ierr = MatCreateShell(this->comm().get(),
323  shell_matrix_A.m(), // Specify the number of local rows
324  shell_matrix_A.n(), // Specify the number of local columns
325  PETSC_DETERMINE,
326  PETSC_DETERMINE,
327  const_cast<void*>(static_cast<const void*>(&shell_matrix_A)),
328  &mat_A);
329 
330  PetscMatrix<T>* matrix_B = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_B_in);
331 
332  // Close the matrix and vectors in case this wasn't already done.
333  matrix_B->close ();
334 
335  /* Note that the const_cast above is only necessary because PETSc
336  does not accept a const void*. Inside the member function
337  _petsc_shell_matrix() below, the pointer is casted back to a
338  const ShellMatrix<T>*. */
339 
340  LIBMESH_CHKERRABORT(ierr);
341  ierr = MatShellSetOperation(mat_A,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
342  ierr = MatShellSetOperation(mat_A,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
343  LIBMESH_CHKERRABORT(ierr);
344 
345  return _solve_generalized_helper (mat_A, matrix_B->mat(), nev, ncv, tol, m_its);
346 }
347 
348 template <typename T>
349 std::pair<unsigned int, unsigned int>
351  ShellMatrix<T> &shell_matrix_B,
352  int nev, // number of requested eigenpairs
353  int ncv, // number of basis vectors
354  const double tol, // solver tolerance
355  const unsigned int m_its) // maximum number of iterations
356 {
357  this->init ();
358 
360 
361  PetscMatrix<T>* matrix_A = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_A_in);
362 
363  // Close the matrix and vectors in case this wasn't already done.
364  matrix_A->close ();
365 
366  // Prepare the matrix.
367  Mat mat_B;
368  ierr = MatCreateShell(this->comm().get(),
369  shell_matrix_B.m(), // Specify the number of local rows
370  shell_matrix_B.n(), // Specify the number of local columns
371  PETSC_DETERMINE,
372  PETSC_DETERMINE,
373  const_cast<void*>(static_cast<const void*>(&shell_matrix_B)),
374  &mat_B);
375 
376 
377  /* Note that the const_cast above is only necessary because PETSc
378  does not accept a const void*. Inside the member function
379  _petsc_shell_matrix() below, the pointer is casted back to a
380  const ShellMatrix<T>*. */
381 
382  LIBMESH_CHKERRABORT(ierr);
383  ierr = MatShellSetOperation(mat_B,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
384  ierr = MatShellSetOperation(mat_B,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
385  LIBMESH_CHKERRABORT(ierr);
386 
387  return _solve_generalized_helper (matrix_A->mat(), mat_B, nev, ncv, tol, m_its);
388 }
389 
390 template <typename T>
391 std::pair<unsigned int, unsigned int>
393  ShellMatrix<T> &shell_matrix_B,
394  int nev, // number of requested eigenpairs
395  int ncv, // number of basis vectors
396  const double tol, // solver tolerance
397  const unsigned int m_its) // maximum number of iterations
398 {
399  this->init ();
400 
402 
403  // Prepare the matrix.
404  Mat mat_A;
405  ierr = MatCreateShell(this->comm().get(),
406  shell_matrix_A.m(), // Specify the number of local rows
407  shell_matrix_A.n(), // Specify the number of local columns
408  PETSC_DETERMINE,
409  PETSC_DETERMINE,
410  const_cast<void*>(static_cast<const void*>(&shell_matrix_A)),
411  &mat_A);
412 
413  Mat mat_B;
414  ierr = MatCreateShell(this->comm().get(),
415  shell_matrix_B.m(), // Specify the number of local rows
416  shell_matrix_B.n(), // Specify the number of local columns
417  PETSC_DETERMINE,
418  PETSC_DETERMINE,
419  const_cast<void*>(static_cast<const void*>(&shell_matrix_B)),
420  &mat_B);
421 
422  /* Note that the const_cast above is only necessary because PETSc
423  does not accept a const void*. Inside the member function
424  _petsc_shell_matrix() below, the pointer is casted back to a
425  const ShellMatrix<T>*. */
426 
427  LIBMESH_CHKERRABORT(ierr);
428  ierr = MatShellSetOperation(mat_A,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
429  ierr = MatShellSetOperation(mat_A,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
430  LIBMESH_CHKERRABORT(ierr);
431 
432  LIBMESH_CHKERRABORT(ierr);
433  ierr = MatShellSetOperation(mat_B,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
434  ierr = MatShellSetOperation(mat_B,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
435  LIBMESH_CHKERRABORT(ierr);
436 
437  return _solve_generalized_helper (mat_A, mat_B, nev, ncv, tol, m_its);
438 }
439 
440 
441 
442 template <typename T>
443 std::pair<unsigned int, unsigned int>
445  Mat mat_B,
446  int nev, // number of requested eigenpairs
447  int ncv, // number of basis vectors
448  const double tol, // solver tolerance
449  const unsigned int m_its) // maximum number of iterations
450 {
451  START_LOG("solve_generalized()", "SlepcEigenSolver");
452 
454 
455  // converged eigen pairs and number of iterations
456  PetscInt nconv=0;
457  PetscInt its=0;
458 
459 #ifdef DEBUG
460  // The relative error.
461  PetscReal error, re, im;
462 
463  // Pointer to vectors of the real parts, imaginary parts.
464  PetscScalar kr, ki;
465 #endif
466 
467  // Set operators.
468  ierr = EPSSetOperators (_eps, mat_A, mat_B);
469  LIBMESH_CHKERRABORT(ierr);
470 
471  //set the problem type and the position of the spectrum
472  set_slepc_problem_type();
473  set_slepc_position_of_spectrum();
474 
475  // Set eigenvalues to be computed.
476 #if SLEPC_VERSION_LESS_THAN(3,0,0)
477  ierr = EPSSetDimensions (_eps, nev, ncv);
478 #else
479  ierr = EPSSetDimensions (_eps, nev, ncv, PETSC_DECIDE);
480 #endif
481  LIBMESH_CHKERRABORT(ierr);
482 
483 
484  // Set the tolerance and maximum iterations.
485  ierr = EPSSetTolerances (_eps, tol, m_its);
486  LIBMESH_CHKERRABORT(ierr);
487 
488  // Set runtime options, e.g.,
489  // -eps_type <type>, -eps_nev <nev>, -eps_ncv <ncv>
490  // Similar to PETSc, these options will override those specified
491  // above as long as EPSSetFromOptions() is called _after_ any
492  // other customization routines.
493  ierr = EPSSetFromOptions (_eps);
494  LIBMESH_CHKERRABORT(ierr);
495 
496  // Solve the eigenproblem.
497  ierr = EPSSolve (_eps);
498  LIBMESH_CHKERRABORT(ierr);
499 
500  // Get the number of iterations.
501  ierr = EPSGetIterationNumber (_eps, &its);
502  LIBMESH_CHKERRABORT(ierr);
503 
504  // Get number of converged eigenpairs.
505  ierr = EPSGetConverged(_eps,&nconv);
506  LIBMESH_CHKERRABORT(ierr);
507 
508 
509 #ifdef DEBUG
510  // ierr = PetscPrintf(this->comm().get(),
511  // "\n Number of iterations: %d\n"
512  // " Number of converged eigenpairs: %d\n\n", its, nconv);
513 
514  // Display eigenvalues and relative errors.
515  ierr = PetscPrintf(this->comm().get(),
516  " k ||Ax-kx||/|kx|\n"
517  " ----------------- -----------------\n" );
518  LIBMESH_CHKERRABORT(ierr);
519 
520  for(PetscInt i=0; i<nconv; i++ )
521  {
522  ierr = EPSGetEigenpair(_eps, i, &kr, &ki, PETSC_NULL, PETSC_NULL);
523  LIBMESH_CHKERRABORT(ierr);
524 
525  ierr = EPSComputeRelativeError(_eps, i, &error);
526  LIBMESH_CHKERRABORT(ierr);
527 
528 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
529  re = PetscRealPart(kr);
530  im = PetscImaginaryPart(kr);
531 #else
532  re = kr;
533  im = ki;
534 #endif
535 
536  if (im != .0)
537  {
538  ierr = PetscPrintf(this->comm().get()," %9f%+9f i %12f\n", re, im, error);
539  LIBMESH_CHKERRABORT(ierr);
540  }
541  else
542  {
543  ierr = PetscPrintf(this->comm().get()," %12f %12f\n", re, error);
544  LIBMESH_CHKERRABORT(ierr);
545  }
546  }
547 
548  ierr = PetscPrintf(this->comm().get(),"\n" );
549  LIBMESH_CHKERRABORT(ierr);
550 #endif // DEBUG
551 
552  STOP_LOG("solve_generalized()", "SlepcEigenSolver");
553 
554  // return the number of converged eigenpairs
555  // and the number of iterations
556  return std::make_pair(nconv, its);
557 
558 }
559 
560 
561 
562 
563 
564 
565 
566 
567 
568 
569 
570 template <typename T>
572 {
573  PetscErrorCode ierr = 0;
574 
575  switch (this->_eigen_solver_type)
576  {
577  case POWER:
578  ierr = EPSSetType (_eps, (char*) EPSPOWER); LIBMESH_CHKERRABORT(ierr); return;
579  case SUBSPACE:
580  ierr = EPSSetType (_eps, (char*) EPSSUBSPACE); LIBMESH_CHKERRABORT(ierr); return;
581  case LAPACK:
582  ierr = EPSSetType (_eps, (char*) EPSLAPACK); LIBMESH_CHKERRABORT(ierr); return;
583  case ARNOLDI:
584  ierr = EPSSetType (_eps, (char*) EPSARNOLDI); LIBMESH_CHKERRABORT(ierr); return;
585  case LANCZOS:
586  ierr = EPSSetType (_eps, (char*) EPSLANCZOS); LIBMESH_CHKERRABORT(ierr); return;
587 #if !SLEPC_VERSION_LESS_THAN(2,3,2)
588  // EPSKRYLOVSCHUR added in 2.3.2
589  case KRYLOVSCHUR:
590  ierr = EPSSetType (_eps, (char*) EPSKRYLOVSCHUR); LIBMESH_CHKERRABORT(ierr); return;
591 #endif
592  // case ARPACK:
593  // ierr = EPSSetType (_eps, (char*) EPSARPACK); LIBMESH_CHKERRABORT(ierr); return;
594 
595  default:
596  libMesh::err << "ERROR: Unsupported SLEPc Eigen Solver: "
597  << Utility::enum_to_string(this->_eigen_solver_type) << std::endl
598  << "Continuing with SLEPc defaults" << std::endl;
599  }
600 }
601 
602 
603 
604 
605 template <typename T>
607 {
608  PetscErrorCode ierr = 0;
609 
610  switch (this->_eigen_problem_type)
611  {
612  case NHEP:
613  ierr = EPSSetProblemType (_eps, EPS_NHEP); LIBMESH_CHKERRABORT(ierr); return;
614  case GNHEP:
615  ierr = EPSSetProblemType (_eps, EPS_GNHEP); LIBMESH_CHKERRABORT(ierr); return;
616  case HEP:
617  ierr = EPSSetProblemType (_eps, EPS_HEP); LIBMESH_CHKERRABORT(ierr); return;
618  case GHEP:
619  ierr = EPSSetProblemType (_eps, EPS_GHEP); LIBMESH_CHKERRABORT(ierr); return;
620 
621  default:
622  libMesh::err << "ERROR: Unsupported SLEPc Eigen Problem: "
623  << this->_eigen_problem_type << std::endl
624  << "Continuing with SLEPc defaults" << std::endl;
625  }
626 }
627 
628 
629 
630 template <typename T>
632 {
633  PetscErrorCode ierr = 0;
634 
635  switch (this->_position_of_spectrum)
636  {
637  case LARGEST_MAGNITUDE:
638  ierr = EPSSetWhichEigenpairs (_eps, EPS_LARGEST_MAGNITUDE); LIBMESH_CHKERRABORT(ierr); return;
639  case SMALLEST_MAGNITUDE:
640  ierr = EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_MAGNITUDE); LIBMESH_CHKERRABORT(ierr); return;
641  case LARGEST_REAL:
642  ierr = EPSSetWhichEigenpairs (_eps, EPS_LARGEST_REAL); LIBMESH_CHKERRABORT(ierr); return;
643  case SMALLEST_REAL:
644  ierr = EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_REAL); LIBMESH_CHKERRABORT(ierr); return;
645  case LARGEST_IMAGINARY:
646  ierr = EPSSetWhichEigenpairs (_eps, EPS_LARGEST_IMAGINARY); LIBMESH_CHKERRABORT(ierr); return;
647  case SMALLEST_IMAGINARY:
648  ierr = EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_IMAGINARY); LIBMESH_CHKERRABORT(ierr); return;
649 
650 
651  default:
652  libMesh::err << "ERROR: Unsupported SLEPc position of spectrum: "
653  << this->_position_of_spectrum << std::endl;
654  libmesh_error();
655  }
656 }
657 
658 
659 
660 
661 
662 
663 template <typename T>
664 std::pair<Real, Real> SlepcEigenSolver<T>::get_eigenpair(unsigned int i,
665  NumericVector<T> &solution_in)
666 {
668 
669  PetscReal re, im;
670 
671  // Make sure the NumericVector passed in is really a PetscVector
672  PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
673 
674  // real and imaginary part of the ith eigenvalue.
675  PetscScalar kr, ki;
676 
677  solution->close();
678 
679  ierr = EPSGetEigenpair(_eps, i, &kr, &ki, solution->vec(), PETSC_NULL);
680  LIBMESH_CHKERRABORT(ierr);
681 
682 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
683  re = PetscRealPart(kr);
684  im = PetscImaginaryPart(kr);
685 #else
686  re = kr;
687  im = ki;
688 #endif
689 
690  return std::make_pair(re, im);
691 }
692 
693 
694 template <typename T>
695 std::pair<Real, Real> SlepcEigenSolver<T>::get_eigenvalue(unsigned int i)
696 {
698 
699  PetscReal re, im;
700 
701  // real and imaginary part of the ith eigenvalue.
702  PetscScalar kr, ki;
703 
704  ierr = EPSGetEigenvalue(_eps, i, &kr, &ki);
705  LIBMESH_CHKERRABORT(ierr);
706 
707 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
708  re = PetscRealPart(kr);
709  im = PetscImaginaryPart(kr);
710 #else
711  re = kr;
712  im = ki;
713 #endif
714 
715  return std::make_pair(re, im);
716 }
717 
718 
719 template <typename T>
721 {
723  PetscReal error;
724 
725  ierr = EPSComputeRelativeError(_eps, i, &error);
726  LIBMESH_CHKERRABORT(ierr);
727 
728  return error;
729 }
730 
731 
732 template <typename T>
734 {
735  this->init();
736 
737  PetscErrorCode ierr = 0;
738  Vec deflation_vector = (libmesh_cast_ptr<PetscVector<T>*>(&deflation_vector_in))->vec();
739  Vec* deflation_space = &deflation_vector;
740 #if SLEPC_VERSION_LESS_THAN(3,1,0)
741  ierr = EPSAttachDeflationSpace(_eps, 1, deflation_space, PETSC_FALSE);
742 #else
743  ierr = EPSSetDeflationSpace(_eps, 1, deflation_space);
744 #endif
745  LIBMESH_CHKERRABORT(ierr);
746 }
747 
748 template <typename T>
750 {
751  /* Get the matrix context. */
753  void* ctx;
754  ierr = MatShellGetContext(mat,&ctx);
755 
757  PetscObjectGetComm((PetscObject)mat,&comm);
758  CHKERRABORT(comm,ierr);
759 
760  /* Get user shell matrix object. */
761  const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
762 
763  /* Make \p NumericVector instances around the vectors. */
764  PetscVector<T> arg_global(arg, shell_matrix.comm());
765  PetscVector<T> dest_global(dest, shell_matrix.comm());
766 
767  /* Call the user function. */
768  shell_matrix.vector_mult(dest_global,arg_global);
769 
770  return ierr;
771 }
772 
773 template <typename T>
775 {
776  /* Get the matrix context. */
778  void* ctx;
779  ierr = MatShellGetContext(mat,&ctx);
780 
782  PetscObjectGetComm((PetscObject)mat,&comm);
783  CHKERRABORT(comm,ierr);
784 
785  /* Get user shell matrix object. */
786  const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
787 
788  /* Make \p NumericVector instances around the vector. */
789  PetscVector<T> dest_global(dest, shell_matrix.comm());
790 
791  /* Call the user function. */
792  shell_matrix.get_diagonal(dest_global);
793 
794  return ierr;
795 }
796 
797 //------------------------------------------------------------------
798 // Explicit instantiations
799 template class SlepcEigenSolver<Number>;
800 
801 } // namespace libMesh
802 
803 
804 
805 #endif // #ifdef LIBMESH_HAVE_SLEPC

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo