include/deal.II/lac/slepc_solver.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: slepc_solver.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //    Author: Toby D. Young, Polish Academy of Sciences, 2008, 2009
00004 //
00005 //    Copyright (C) 2009, 2010, 2012 by the deal.II authors
00006 //
00007 //    This file is subject to QPL and may not be  distributed
00008 //    without copyright and license information. Please refer
00009 //    to the file deal.II/doc/license.html for the  text  and
00010 //    further information on this license.
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   // --------------------------- inline and template functions -----------
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 /*----------------------------   slepc_solver.h  ---------------------------*/
00734 
00735 #endif
00736 
00737 /*----------------------------   slepc_solver.h  ---------------------------*/
00738 
00739 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:39 by doxygen 1.7.3