00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #ifndef __deal2__slepc_solver_h
00015 #define __deal2__slepc_solver_h
00016
00017
00018 #include <deal.II/base/config.h>
00019
00020 #ifdef DEAL_II_USE_SLEPC
00021
00022 # include <deal.II/base/std_cxx1x/shared_ptr.h>
00023 # include <deal.II/lac/exceptions.h>
00024 # include <deal.II/lac/solver_control.h>
00025 # include <deal.II/lac/petsc_matrix_base.h>
00026 # include <deal.II/lac/slepc_spectral_transformation.h>
00027
00028 # include <petscksp.h>
00029 # include <slepceps.h>
00030
00031 DEAL_II_NAMESPACE_OPEN
00032
00093 namespace SLEPcWrappers
00094 {
00095
00110 class SolverBase
00111 {
00112 public:
00113
00119 SolverBase (SolverControl &cn,
00120 const MPI_Comm &mpi_communicator);
00121
00125 virtual ~SolverBase ();
00126
00164 template <typename OutputVector>
00165 void
00166 solve (const PETScWrappers::MatrixBase &A,
00167 std::vector<double> &kr,
00168 std::vector<OutputVector> &vr,
00169 const unsigned int n_eigenvectors);
00170
00176 template <typename OutputVector>
00177 void
00178 solve (const PETScWrappers::MatrixBase &A,
00179 const PETScWrappers::MatrixBase &B,
00180 std::vector<double> &kr,
00181 std::vector<OutputVector> &vr,
00182 const unsigned int n_eigenvectors);
00183
00190 void
00191 set_matrices (const PETScWrappers::MatrixBase &A);
00192
00198 void
00199 set_matrices (const PETScWrappers::MatrixBase &A,
00200 const PETScWrappers::MatrixBase &B);
00201
00206 void
00207 set_initial_vector
00208 (const PETScWrappers::VectorBase &this_initial_vector);
00209
00214 void
00215 set_transformation (SLEPcWrappers::TransformationBase &this_transformation);
00216
00225 void
00226 set_which_eigenpairs (EPSWhich set_which);
00227
00240 void
00241 solve (const unsigned int n_eigenvectors, unsigned int *n_converged);
00242
00243
00251 void
00252 get_eigenpair (const unsigned int index,
00253 #ifndef DEAL_II_USE_PETSC_COMPLEX
00254 double &kr,
00255 #else
00256 std::complex<double> &kr,
00257 #endif
00258 PETScWrappers::VectorBase &vr);
00259
00264 void
00265 reset();
00266
00271 EPS *
00272 get_EPS ();
00273
00274
00279 SolverControl &control () const;
00280
00284 DeclException0 (ExcSLEPcWrappersUsageError);
00285
00286 DeclException1 (ExcSLEPcError,
00287 int,
00288 << " An error with error number " << arg1
00289 << " occured while calling a SLEPc function");
00290
00291 DeclException2 (ExcSLEPcEigenvectorConvergenceMismatchError,
00292 int, int,
00293 << " The number of converged eigenvectors is " << arg1
00294 << " but " << arg2 << " were requested. ");
00295
00296 protected:
00297
00303 SolverControl &solver_control;
00304
00310 const MPI_Comm mpi_communicator;
00311
00319 virtual void set_solver_type (EPS &eps) const = 0;
00320
00325 EPSWhich set_which;
00326
00334 const PETScWrappers::MatrixBase *opA;
00335
00341 const PETScWrappers::MatrixBase *opB;
00342
00343 const PETScWrappers::VectorBase *initial_vector;
00344
00351 SLEPcWrappers::TransformationBase *transformation;
00352
00353 private:
00354
00364 static
00365 int
00366 convergence_test (EPS eps,
00367 #ifdef PETSC_USE_64BIT_INDICES
00368 const PetscInt iteration,
00369 #else
00370 const int iteration,
00371 #endif
00372 const PetscReal residual_norm,
00373 EPSConvergedReason *reason,
00374 void *solver_control);
00375
00389 struct SolverData
00390 {
00391
00395 ~SolverData ();
00396
00401 EPS eps;
00402 };
00403
00404 std_cxx1x::shared_ptr<SolverData> solver_data;
00405 };
00406
00415 class SolverKrylovSchur : public SolverBase
00416 {
00417 public:
00418
00424 struct AdditionalData
00425 {};
00426
00436 SolverKrylovSchur (SolverControl &cn,
00437 const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
00438 const AdditionalData &data = AdditionalData());
00439
00440 protected:
00441
00446 const AdditionalData additional_data;
00447
00454 virtual void set_solver_type (EPS &eps) const;
00455 };
00456
00464 class SolverArnoldi : public SolverBase
00465 {
00466 public:
00472 struct AdditionalData
00473 {
00480 AdditionalData (const bool delayed_reorthogonalization = false);
00481
00486 bool delayed_reorthogonalization;
00487 };
00488
00498 SolverArnoldi (SolverControl &cn,
00499 const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
00500 const AdditionalData &data = AdditionalData());
00501
00502 protected:
00503
00508 const AdditionalData additional_data;
00509
00516 virtual void set_solver_type (EPS &eps) const;
00517
00518 };
00519
00527 class SolverLanczos : public SolverBase
00528 {
00529 public:
00535 struct AdditionalData
00536 {};
00537
00547 SolverLanczos (SolverControl &cn,
00548 const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
00549 const AdditionalData &data = AdditionalData());
00550
00551 protected:
00552
00557 const AdditionalData additional_data;
00558
00565 virtual void set_solver_type (EPS &eps) const;
00566
00567 };
00568
00577 class SolverPower : public SolverBase
00578 {
00579 public:
00585 struct AdditionalData
00586 {};
00587
00597 SolverPower (SolverControl &cn,
00598 const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
00599 const AdditionalData &data = AdditionalData());
00600
00601 protected:
00602
00607 const AdditionalData additional_data;
00608
00615 virtual void set_solver_type (EPS &eps) const;
00616
00617 };
00618
00619 #if DEAL_II_PETSC_VERSION_GTE(3,1,0)
00620
00627 class SolverDavidson : public SolverBase
00628 {
00629 public:
00635 struct AdditionalData
00636 {};
00637
00647 SolverDavidson (SolverControl &cn,
00648 const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
00649 const AdditionalData &data = AdditionalData());
00650
00651 protected:
00652
00657 const AdditionalData additional_data;
00658
00665 virtual void set_solver_type (EPS &eps) const;
00666
00667 };
00668 #endif
00669
00670
00676 template <typename OutputVector>
00677 void
00678 SolverBase::solve (const PETScWrappers::MatrixBase &A,
00679 std::vector<double> &kr,
00680 std::vector<OutputVector> &vr,
00681 const unsigned int n_eigenvectors = 1)
00682 {
00683 unsigned int n_converged;
00684
00685 set_matrices (A);
00686 solve (n_eigenvectors,&n_converged);
00687
00688 if (n_converged > n_eigenvectors)
00689 n_converged = n_eigenvectors;
00690 AssertThrow (n_converged == n_eigenvectors,
00691 ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenvectors));
00692
00693 AssertThrow (vr.size() >= 1, ExcSLEPcWrappersUsageError());
00694 vr.resize (n_converged, vr.front());
00695 kr.resize (n_converged);
00696
00697 for (unsigned int index=0; index < n_converged; ++index)
00698 get_eigenpair (index, kr[index], vr[index]);
00699 }
00700
00701 template <typename OutputVector>
00702 void
00703 SolverBase::solve (const PETScWrappers::MatrixBase &A,
00704 const PETScWrappers::MatrixBase &B,
00705 std::vector<double> &kr,
00706 std::vector<OutputVector> &vr,
00707 const unsigned int n_eigenvectors = 1)
00708 {
00709 unsigned int n_converged;
00710
00711 set_matrices (A, B);
00712 solve (n_eigenvectors, &n_converged);
00713
00714 if (n_converged >= n_eigenvectors)
00715 n_converged = n_eigenvectors;
00716
00717 AssertThrow (n_converged == n_eigenvectors,
00718 ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenvectors));
00719 AssertThrow (vr.size() >= 1, ExcSLEPcWrappersUsageError());
00720
00721 vr.resize (n_converged, vr.front());
00722 kr.resize (n_converged);
00723
00724 for (unsigned int index=0; index < n_converged; ++index)
00725 get_eigenpair (index, kr[index], vr[index]);
00726 }
00727 }
00728
00729 DEAL_II_NAMESPACE_CLOSE
00730
00731 #endif // DEAL_II_USE_SLEPC
00732
00733
00734
00735 #endif
00736
00737
00738
00739