]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added SLEPcWrappers which give a handle on SLEPc: the resultant two base classes...
authorToby D. Young <tyoung@ippt.pan.pl>
Wed, 24 Jun 2009 12:24:08 +0000 (12:24 +0000)
committerToby D. Young <tyoung@ippt.pan.pl>
Wed, 24 Jun 2009 12:24:08 +0000 (12:24 +0000)
git-svn-id: https://svn.dealii.org/trunk@18957 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/lac/include/lac/slepc_solver.h [new file with mode: 0755]
deal.II/lac/include/lac/slepc_spectral_transformation.h [new file with mode: 0755]
deal.II/lac/source/slepc_solver.cc [new file with mode: 0644]
deal.II/lac/source/slepc_spectral_transformation.cc [new file with mode: 0644]

index eed45d43f2e65320a143702545d3d30002249770..e7b452b24c33702ca044e48ba881fed9f5290e7b 100644 (file)
@@ -129,6 +129,15 @@ inconvenience this causes.
   </li>
 </ol>
 
+<ol>
+  <li>
+  <p>
+  New: Based on work with Rickard Armiento, Francisco Alvaro, and Jose E. Roman, SLEPcWrappers that give a handle on some of the features of SLEPc (Scalable Library for Eigenvalue Problem Computations): (1) The SLEPcWrappers::SolverBase class can be used for specifying an eigenvalue problem, either in standard or generalized form, on serial or parallel architectures with support for a few solver types; and (2) The SLEPcWrappers::TransformationBase class encapsulates a variety of spectral transformations providing some functionality required for acceleration techniques based on the transformation of the spectrum.
+  <br>
+  (Toby D. Young 2009/06/25)
+  </p>
+  </li>
+</ol>
 
 
 <a name="deal.II"></a>
diff --git a/deal.II/lac/include/lac/slepc_solver.h b/deal.II/lac/include/lac/slepc_solver.h
new file mode 100755 (executable)
index 0000000..fef1bbe
--- /dev/null
@@ -0,0 +1,510 @@
+
+#ifndef __deal2__slepc_solver_h
+#define __deal2__slepc_solver_h
+
+#include <base/config.h>
+#include <lac/exceptions.h>
+#include <lac/solver_control.h>
+#include <lac/slepc_spectral_transformation.h>
+#include <boost/shared_ptr.hpp>
+
+#ifdef DEAL_II_USE_SLEPC
+
+#include <petscksp.h>
+#include <slepceps.h>
+
+DEAL_II_NAMESPACE_OPEN 
+
+/**
+ * Base class for solver classes using the SLEPc solvers which are
+ * selected based on flags passed to the eigenvalue problem solver
+ * context. Derived classes set the right flags to set the right
+ * solver -- in particular, note that the AdditionalData structure is
+ * a dummy structure and is there for backward compatibility.
+ *
+ * SLEPcWrappers can be implemented in application codes in the
+ * following way:
+ @verbatim
+  SolverControl solver_control (1000, 1e-10);
+  SolverArnoldi system (solver_control, 
+                        mpi_communicator);
+  system.solve (A, B, lambda, x, n_eigenvectors);
+ @endverbatim
+ * for the generalized eigenvalue problem $Ax=B\lambda x$.
+ *
+ * An alternative implementation to the one above is to use the API
+ * internals directly within the application code. In this way the
+ * calling sequence requires calling several of SolverBase functions
+ * rather than just one. This freedom is intended for use by users of
+ * the SLEPcWrappers that require a greater handle on the eigenvalue
+ * problem solver context. See also:
+ @verbatim
+  template <typename OutputVector>
+  void
+  SolverBase::solve (const PETScWrappers::MatrixBase &A,
+                    const PETScWrappers::MatrixBase &B,
+                    std::vector<double>             &kr,
+                    std::vector<OutputVector>       &vr,
+                    const unsigned int               n_eigenvectors
+ )
+ @endverbatim
+ * as an example on how to do this. 
+ *
+ * See also the @ref PETScWrappers "PETScWrappers", on which the
+ * SLEPcWrappers depend, for additional explanations.
+ *
+ * @ingroup SLEPcWrappers
+ * @author Toby D. Young 2008, 2009
+ */
+namespace SLEPcWrappers
+{
+
+                                   /**
+                                    * Base class for solver classes
+                                    * using the SLEPc solvers. Since
+                                    * solvers in SLEPc are selected
+                                    * based on flags passed to a
+                                    * generic solver object, basically
+                                    * all the actual solver calls
+                                    * happen in this class, and
+                                    * derived classes simply set the
+                                    * right flags to select one solver
+                                    * or another, or to set certain
+                                    * parameters for individual
+                                    * solvers.
+                                    */
+  class SolverBase
+    {
+    public:
+
+                                   /**
+                                    * Constructor. Takes the MPI
+                                    * communicator over which parallel
+                                    * computations are to happen.
+                                    */
+      SolverBase (SolverControl &cn,
+                 const MPI_Comm &mpi_communicator);
+      
+                                   /**
+                                    * Destructor.
+                                    */
+      virtual ~SolverBase ();      
+
+                                   /**
+                                    * Composite method that solves the
+                                    * linear system $Ax=\lambda
+                                    * Bx$. The eigenvector sent in has
+                                    * to have at least one element
+                                    * that we can use as a template
+                                    * when resizing, since we do not
+                                    * know the parameters of the
+                                    * specific vector class used
+                                    * (i.e. local_dofs for MPI
+                                    * vectors). However, while copying
+                                    * eigenvectors, at least twice the
+                                    * memory size of <tt>vr</tt> is
+                                    * being used (and can be more). To
+                                    * avoid doing this, use instead
+                                    * the calling sequence used here
+                                    * is fairly standard: Initialise;
+                                    * set up matrices for solving;
+                                    * actually solve the system;
+                                    * gather the solution(s); and
+                                    * reset.
+                                    * 
+                                    * Note that the number of
+                                    * converged eigenstates can be
+                                    * larger than the number of
+                                    * eigenstates requested; this is
+                                    * due to a round off error
+                                    * (success) of the eigenvalue
+                                    * solver context. If this is found
+                                    * to be the case, we simply do not
+                                    * bother with more eigenpairs than
+                                    * requested but handle that it may
+                                    * be more by ignoring any extras.
+                                   */
+      template <typename OutputVector>
+       void
+       solve (const PETScWrappers::MatrixBase &A,
+              const PETScWrappers::MatrixBase &B,
+              std::vector<double>             &kr,
+              std::vector<OutputVector>       &vr,
+              const unsigned int               n_eigenvectors);
+      
+                                   /**
+                                    * Initialize solver for the linear
+                                    * system $Ax=\lambda
+                                    * Bx$. (required before calling
+                                    * solve)
+                                    */
+      void
+       set_matrices (const PETScWrappers::MatrixBase &A,
+                     const PETScWrappers::MatrixBase &B);   
+
+                                   /**
+                                    * Set the initial vector for the solver.
+                                    */
+      void
+       set_initial_vector (const PETScWrappers::VectorBase &initial_vec);
+
+                                   /**
+                                   * Set the spectral transformation
+                                   * to be used.  By default SLEPc
+                                    */
+      void
+       set_transformation (SLEPcWrappers::TransformationBase &trans);
+
+                                   /**
+                                   * Indicate which part of the
+                                   * spectrum is to be computed. By
+                                   * default largest magnitude
+                                   * eigenvalues are computed.  For
+                                   * other allowed values see the
+                                   * SLEPc documentation.
+                                    */
+      void
+       set_which_eigenpairs (EPSWhich set_which);
+
+                                   /**
+                                    * Solve the linear system for
+                                   * n_eigenvectors
+                                   * eigenstates. Parameter
+                                   * n_converged contains the actual
+                                   * number of eigenstates that have
+                                   * .  converged; this can be both
+                                   * fewer or more than
+                                   * n_eigenvectors, depending on the
+                                   * SLEPc eigensolver used.
+                                    */
+      void
+       solve (const unsigned int n_eigenvectors, unsigned int *n_converged);   
+
+      
+                                   /**
+                                    * Access the solutions for a
+                                    * solved eigenvector problem, pair
+                                    * index solutions, index = 0
+                                    * ... n_converged-1
+                                    */
+      void
+       get_eigenpair (const unsigned int index,
+                      double &kr,
+                      PETScWrappers::VectorBase &vr);
+      
+                                   /**
+                                    * Reset the solver, and return
+                                    * memory for eigenvectors
+                                    */
+      void
+       reset();
+
+                                   /**
+                                   * Retrieve the SLEPc solver object
+                                   * used internally.
+                                    */
+      EPS *
+       get_EPS ();
+
+
+                                   /**
+                                    * Access to object that controls
+                                    * convergence.
+                                    */
+      SolverControl & control() const;
+      
+                                   /**
+                                    * Exceptions.
+                                    */
+      DeclException0 (ExcSLEPcWrappersUsageError);
+      DeclException1 (ExcSLEPcError, 
+                     int,
+                     << "An error with error number " << arg1
+                     << " occured while calling a SLEPc function");
+      
+    protected:
+      
+                                   /**
+                                    * Reference to the object that
+                                    * controls convergence of the
+                                    * iterative solver.
+                                    */
+      SolverControl &solver_control;
+
+                                   /**
+                                    * Copy of the MPI communicator
+                                    * object to be used for the
+                                    * solver.
+                                    */
+      const MPI_Comm mpi_communicator;  
+
+                                   /**
+                                    * Function that takes an
+                                    * Eigenvalue Problem Solver
+                                    * context object, and sets the
+                                    * type of solver that is requested
+                                    * by the derived class.
+                                    */
+      virtual void set_solver_type (EPS &eps) const = 0;
+
+                                   /**
+                                    * Attributes that store the
+                                   * relevant information for the
+                                   * eigenproblem solver context.
+                                    */
+      EPSWhich                           set_which;
+      const PETScWrappers::MatrixBase   *opA;
+      const PETScWrappers::MatrixBase   *opB;
+      const PETScWrappers::VectorBase   *ini_vec;
+      SLEPcWrappers::TransformationBase *transform;
+      
+    private:
+
+                                   /**
+                                    * A function that is used in SLEPc
+                                    * as a callback to check on
+                                    * convergence. It takes the
+                                    * information provided from SLEPc
+                                    * and checks it against deal.II's
+                                    * own SolverControl objects to see
+                                    * if convergence has been reached.
+                                    */
+      static
+       int
+       convergence_test (EPS                 eps,
+                         const int           iteration,
+                         const PetscScalar   residual_norm,
+                         EPSConvergedReason *reason,
+                         void               *solver_control);
+
+                                   /**
+                                    * Objects of this type are
+                                    * explicitly created, but are
+                                    * destroyed when the surrounding
+                                    * solver object goes out of scope,
+                                    * or when we assign a new value to
+                                    * the pointer to this object. The
+                                    * respective Destroy functions are
+                                    * therefore written into the
+                                    * destructor of this object, even
+                                    * though the object does not have
+                                    * a constructor.
+                                    */
+      struct SolverData
+      {
+
+                                   /**
+                                   * Destructor.
+                                   */
+       ~SolverData ();
+
+                                   /**
+                                    * Objects for Eigenvalue Problem
+                                    * Solver.
+                                    */
+       EPS eps;  
+      };
+
+      boost::shared_ptr<SolverData> solver_data;
+    };
+
+/**
+ * An implementation of the solver interface using the SLEPc
+ * Krylov-Schur solver. Usage: All spectrum, all problem types,
+ * complex.
+ *
+ * @ingroup SLEPcWrappers
+ * @author Toby D. Young 2008
+ */
+  class SolverKrylovSchur : public SolverBase
+    {
+    public:
+
+                                   /**
+                                    * Standardized data struct to pipe
+                                    * additional data to the solver,
+                                    * should it be needed.
+                                    */
+      struct AdditionalData 
+      {};      
+
+                                   /**
+                                    * SLEPc solvers will want to have
+                                    * an MPI communicator context over
+                                    * which computations are
+                                    * parallelized. By default, this
+                                    * carries the same behaviour has
+                                    * the PETScWrappers, but you can
+                                    * change that.
+                                    */
+       SolverKrylovSchur (SolverControl        &cn,
+                         const MPI_Comm       &mpi_communicator = PETSC_COMM_WORLD,
+                         const AdditionalData &data = AdditionalData());
+      
+    protected:
+
+                                   /**
+                                    * Store a copy of the flags for
+                                    * this particular solver.
+                                    */
+      const AdditionalData additional_data;
+
+                                   /**
+                                    * Function that takes a Eigenvalue
+                                    * Problem Solver context object,
+                                    * and sets the type of solver that
+                                    * is appropriate for this class.
+                                    */
+      virtual void set_solver_type (EPS &eps) const; 
+    };
+
+/**
+ * An implementation of the solver interface using the SLEPc Arnoldi
+ * solver. Usage: All spectrum, all problem types, complex.
+ *
+ * @ingroup SLEPcWrappers
+ * @author Toby D. Young 2008
+ */
+  class SolverArnoldi : public SolverBase
+    {
+    public:
+                                   /**
+                                    * Standardized data struct to pipe
+                                    * additional data to the solver,
+                                    * should it be needed.
+                                    */
+      struct AdditionalData 
+      {};      
+
+                                   /**
+                                    * SLEPc solvers will want to have
+                                    * an MPI communicator context over
+                                    * which computations are
+                                    * parallelized. By default, this
+                                    * carries the same behaviour has
+                                    * the PETScWrappers, but you can
+                                    * change that.
+                                    */
+      SolverArnoldi (SolverControl        &cn,
+                    const MPI_Comm       &mpi_communicator = PETSC_COMM_WORLD,
+                    const AdditionalData &data = AdditionalData());
+      
+    protected:
+
+                                   /**
+                                    * Store a copy of the flags for
+                                    * this particular solver.
+                                    */
+      const AdditionalData additional_data;
+
+                                   /**
+                                    * Function that takes a Eigenvalue
+                                    * Problem Solver context object,
+                                    * and sets the type of solver that
+                                    * is appropriate for this class.
+                                    */
+      virtual void set_solver_type (EPS &eps) const; 
+   
+    };
+
+/**
+ * An implementation of the solver interface using the SLEPc Lanczos
+ * solver. Usage: All spectrum, all problem types, complex.
+ *
+ * @ingroup SLEPcWrappers
+ * @author Toby D. Young 2009
+ */
+  class SolverLanczos : public SolverBase
+    {
+    public:
+                                   /**
+                                    * Standardized data struct to pipe
+                                    * additional data to the solver,
+                                    * should it be needed.
+                                    */
+      struct AdditionalData 
+      {};      
+      
+                                   /**
+                                    * SLEPc solvers will want to have
+                                    * an MPI communicator context over
+                                    * which computations are
+                                    * parallelized. By default, this
+                                    * carries the same behaviour has
+                                    * the PETScWrappers, but you can
+                                    * change that.
+                                    */
+      SolverLanczos (SolverControl        &cn,
+                    const MPI_Comm       &mpi_communicator = PETSC_COMM_WORLD,
+                    const AdditionalData &data = AdditionalData());
+      
+    protected:
+
+                                   /**
+                                    * Store a copy of the flags for
+                                    * this particular solver.
+                                    */
+      const AdditionalData additional_data;
+
+                                   /**
+                                    * Function that takes a Eigenvalue
+                                    * Problem Solver context object,
+                                    * and sets the type of solver that
+                                    * is appropriate for this class.
+                                    */
+      virtual void set_solver_type (EPS &eps) const; 
+   
+    };
+
+
+  // --------------------------- inline and template functions -----------                                                                       
+
+  /**                                                                                                           
+   * This is declared here to make it                                                                           
+   * possible to take a std::vector                                                                             
+   * of different PETScWrappers vector                                                                          
+   * types                                                                                                      
+   */
+  template <typename OutputVector>
+    void
+    SolverBase::solve (const PETScWrappers::MatrixBase &A,
+                       const PETScWrappers::MatrixBase &B,
+                       std::vector<double>             &kr,
+                       std::vector<OutputVector>       &vr,
+                       const unsigned int               n_eigenvectors = 0)
+    {
+      unsigned int n_converged;
+
+      set_matrices(A,B);
+
+      solve(n_eigenvectors,&n_converged);
+
+      if (n_converged > n_eigenvectors)
+        {
+          n_converged = n_eigenvectors;
+        }
+
+      AssertThrow (vr.size() >= 1, ExcSLEPcWrappersUsageError());
+      vr.resize(n_converged, vr.front());
+      kr.resize(n_converged);
+
+      for (unsigned int index=0; index < n_converged; 
+          ++index)
+        {
+          get_eigenpair(index, kr[index], vr[index]);
+        }
+    }
+
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_SLEPC
+
+/*----------------------------   slepc_solver.h  ---------------------------*/
+
+#endif 
+
+/*----------------------------   slepc_solver.h  ---------------------------*/
+
diff --git a/deal.II/lac/include/lac/slepc_spectral_transformation.h b/deal.II/lac/include/lac/slepc_spectral_transformation.h
new file mode 100755 (executable)
index 0000000..0b229a3
--- /dev/null
@@ -0,0 +1,317 @@
+
+#ifndef __deal2__slepc_spectral_transformation_h
+#define __deal2__slepc_spectral_transformation_h
+
+#include <base/config.h>
+#include <lac/exceptions.h>
+#include <lac/solver_control.h>
+#include <lac/slepc_solver.h>
+#include <boost/shared_ptr.hpp>
+
+#ifdef DEAL_II_USE_SLEPC
+
+#include <petscksp.h>
+#include <slepceps.h>
+
+DEAL_II_NAMESPACE_OPEN 
+
+
+namespace SLEPcWrappers
+{
+
+/**
+ * Base class for spectral transformation classes using the SLEPc
+ * solvers which are selected based on flags passed to the spectral
+ * transformation.
+ *
+ * @ingroup SLEPcWrappers
+ * @author Toby D. Young 2009
+ **/
+  class TransformationBase
+    {
+    public:
+
+                                   /**
+                                    * Constructor. Takes the MPI
+                                    * communicator over which parallel
+                                    * computations are to happen.
+                                    */
+      TransformationBase ();
+
+                                   /**
+                                    * Destructor.
+                                    */
+      virtual ~TransformationBase ();      
+
+                                   /**
+                                   * Record the EPS object that is associated
+                                   * to the spectral transformation
+                                    */
+      void set_context (EPS &eps);
+
+    protected:
+
+      virtual void set_transformation_type (ST &st) const = 0;
+
+    private:
+
+                                   /**
+                                    * Objects of this type are
+                                    * explicitly created, but are
+                                    * destroyed when the surrounding
+                                    * solver object goes out of scope,
+                                    * or when we assign a new value to
+                                    * the pointer to this object. The
+                                    * respective Destroy functions are
+                                    * therefore written into the
+                                    * destructor of this object, even
+                                    * though the object does not have
+                                    * a constructor.
+                                    */
+      struct TransformationData
+      {
+       
+                                   /**
+                                   * Destructor.
+                                   */
+       ~TransformationData ();
+
+                                   /**
+                                    * Objects for Eigenvalue Problem
+                                    * Solver.
+                                    */
+       ST st;  
+      };
+
+      boost::shared_ptr<TransformationData> transformation_data;
+    };
+
+/**
+ * An implementation of the transformation interface using the SLEPc 
+ * Shift.
+ *
+ * @ingroup SLEPcWrappers
+ * @author Toby D. Young 2009
+ */
+  class TransformationShift : public TransformationBase
+    {
+    public:
+      
+                                   /**
+                                    * Standardized data struct to
+                                    * pipe additional data to the
+                                    * solver.
+                                    */
+      struct AdditionalData
+      {
+
+                                   /**
+                                    * Constructor. By default, set the
+                                    * shift parameter to zero.
+                                    */
+       AdditionalData (const double shift_parameter = 0);
+
+                                   /**
+                                    * Shift parameter.
+                                    */
+       const double shift_parameter;
+      };
+
+        
+                                   /**
+                                    * Constructor.
+                                    */
+      TransformationShift (const AdditionalData &data = AdditionalData());
+
+
+    protected:
+
+                                   /**
+                                    * Store a copy of the flags for this
+                                    * particular solver.
+                                    */
+      const AdditionalData additional_data;
+
+                                   /**
+                                    * Function that takes a Spectral
+                                    * Transformation context object,
+                                    * and sets the type of spectral
+                                    * transformationthat is
+                                    * appropriate for this class.
+                                    */
+      virtual void set_transformation_type (ST &st) const; 
+    };
+
+/**
+ * An implementation of the transformation interface using the SLEPc 
+ * Shift and Invert.
+ *
+ * @ingroup SLEPcWrappers
+ * @author Toby D. Young 2009
+ */
+  class TransformationShiftInvert : public TransformationBase
+    {
+    public:
+     
+                                   /**
+                                    * Standardized data struct to
+                                    * pipe additional data to the
+                                    * solver.
+                                    */
+      struct AdditionalData
+      {
+                                   /**
+                                    * Constructor. By default, set the
+                                    * shift parameter to zero.
+                                    */
+       AdditionalData (const double shift_parameter = 0);
+
+                                   /**
+                                    * Shift parameter.
+                                    */
+       const double shift_parameter;
+      };
+
+        
+                                   /**
+                                    * Constructor.
+                                    */
+      TransformationShiftInvert (const AdditionalData &data = AdditionalData());
+
+    protected:
+
+                                   /**
+                                    * Store a copy of the flags for this
+                                    * particular solver.
+                                    */
+      const AdditionalData additional_data;
+
+                                   /**
+                                    * Function that takes a Spectral
+                                    * Transformation context object,
+                                    * and sets the type of spectral
+                                    * transformationthat is
+                                    * appropriate for this class.
+                                    */
+      virtual void set_transformation_type (ST &st) const; 
+    };
+
+/**
+ * An implementation of the transformation interface using the SLEPc 
+ * Spectrum Folding.
+ *
+ * @ingroup SLEPcWrappers
+ * @author Toby D. Young 2009
+ */
+  class TransformationSpectrumFolding : public TransformationBase
+    {
+    public:
+      
+                                   /**
+                                    * Standardized data struct to
+                                    * pipe additional data to the
+                                    * solver.
+                                    */
+      struct AdditionalData
+      {
+                                   /**
+                                    * Constructor. By default, set the
+                                    * shift parameter to zero.
+                                    */
+       AdditionalData (const double shift_parameter = 0);
+
+                                   /**
+                                    * Shift parameter.
+                                    */
+       const double shift_parameter;
+      };
+
+        
+                                   /**
+                                    * Constructor.
+                                    */
+      TransformationSpectrumFolding (const AdditionalData &data = AdditionalData());
+
+    protected:
+
+                                   /**
+                                    * Store a copy of the flags for this
+                                    * particular solver.
+                                    */
+      const AdditionalData additional_data;
+
+                                   /**
+                                    * Function that takes a Spectral
+                                    * Transformation context object,
+                                    * and sets the type of spectral
+                                    * transformationthat is
+                                    * appropriate for this class.
+                                    */
+      virtual void set_transformation_type (ST &st) const; 
+    };
+
+/**
+ * An implementation of the transformation interface using the SLEPc 
+ * Cayley.
+ *
+ * @ingroup SLEPcWrappers
+ * @author Toby D. Young 2009
+ */
+  class TransformationCayley : public TransformationBase
+    {
+    public:
+      
+                                   /**
+                                    * Standardized data struct to
+                                   * pipe additional data to the
+                                    * solver.
+                                    */
+      struct AdditionalData
+      {
+                                   /**
+                                    * Constructor. Requires two shift parameters
+                                    */
+       AdditionalData (const double shift_parameter = 0, const double antishift_parameter = 0);
+
+                                   /**
+                                    * Shift and antishift parameter.
+                                    */
+       const double shift_parameter;
+       const double antishift_parameter;
+      };
+
+        
+                                   /**
+                                    * Constructor.
+                                    */
+      TransformationCayley (const double shift, const double antishift);
+
+    protected:
+
+                                   /**
+                                    * Store a copy of the flags for this
+                                    * particular solver.
+                                    */
+      const AdditionalData additional_data;
+
+                                   /**
+                                    * Function that takes a Spectral
+                                    * Transformation context object,
+                                    * and sets the type of spectral
+                                    * transformationthat is
+                                    * appropriate for this class.
+                                    */
+      virtual void set_transformation_type (ST &st) const; 
+    };
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_SLEPC
+
+/*--------------------   slepc_spectral_transformation.h   ------------------*/
+
+#endif 
+
+/*--------------------   slepc_spectral_transformation.h   ------------------*/
diff --git a/deal.II/lac/source/slepc_solver.cc b/deal.II/lac/source/slepc_solver.cc
new file mode 100644 (file)
index 0000000..f9042be
--- /dev/null
@@ -0,0 +1,289 @@
+
+#include <lac/petsc_matrix_base.h>
+#include <lac/petsc_vector_base.h>
+#include <lac/petsc_vector.h>
+#include <lac/slepc_solver.h>
+#include <lac/slepc_spectral_transformation.h>
+
+#include <cmath>
+#include <vector>
+
+#ifdef DEAL_II_USE_SLEPC
+
+#if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2)
+#include <petscsles.h>
+#endif
+#include <petscversion.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace SLEPcWrappers
+{
+
+  SolverBase::SolverData::~SolverData ()
+  {
+                                   // Destroy the solver object.
+    int ierr = EPSDestroy (eps);
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr)); 
+  }
+  
+  SolverBase::SolverBase (SolverControl  &cn,
+                          const MPI_Comm &mpi_communicator)
+    :
+    solver_control (cn),
+    mpi_communicator (mpi_communicator),
+    set_which (EPS_LARGEST_MAGNITUDE),
+    opA (NULL), opB (NULL),
+    ini_vec (NULL),
+    transform (NULL)
+  {
+  }
+  
+  SolverBase::~SolverBase ()
+  {
+    if( solver_data != 0 )
+      solver_data.reset ();
+  }
+  
+  void
+  SolverBase::set_matrices (const PETScWrappers::MatrixBase &A,
+                           const PETScWrappers::MatrixBase &B)
+  {
+    opA = &A;
+    opB = &B;
+  }
+
+  void
+  SolverBase::set_initial_vector (const PETScWrappers::VectorBase &initial_vec)
+  {
+    ini_vec = &initial_vec;
+  }
+
+  void
+  SolverBase::set_transformation (SLEPcWrappers::TransformationBase &trans)
+  {
+    transform = &trans;
+  }
+
+  void
+  SolverBase::set_which_eigenpairs (const EPSWhich eps_which)
+  {
+    set_which = eps_which;
+  }
+
+  void
+  SolverBase::solve (const unsigned int n_eigenvectors, unsigned int *n_converged)
+  {
+    int ierr;
+    
+    AssertThrow (solver_data == 0, ExcSLEPcWrappersUsageError());    
+    solver_data.reset (new SolverData());
+    
+                                    // create eigensolver context and
+                                    // set operators.
+    ierr = EPSCreate (mpi_communicator, &solver_data->eps);
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+    AssertThrow (opA && opB, ExcSLEPcWrappersUsageError());
+    ierr = EPSSetOperators (solver_data->eps, *opA, *opB);
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+    if( ini_vec && ini_vec->size() != 0 ) 
+      {
+       ierr = EPSSetInitialVector(solver_data->eps, *ini_vec);
+       AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+      }
+    
+    if( transform )
+      transform->set_context(solver_data->eps);
+
+                                    // set runtime options.
+    set_solver_type (solver_data->eps);
+
+    ierr = EPSSetWhichEigenpairs (solver_data->eps, set_which);
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+                                    // set number of eigenvectors to
+                                    // compute
+    ierr = EPSSetDimensions (solver_data->eps, n_eigenvectors, 
+                            PETSC_DECIDE, PETSC_DECIDE);
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+    ierr = EPSSetFromOptions (solver_data->eps);
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+                                    // solve the eigensystem
+    ierr = EPSSolve (solver_data->eps);
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+                                    // get number of converged
+                                    // eigenstates
+    ierr = EPSGetConverged (solver_data->eps, reinterpret_cast<int *>(n_converged));
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+  }
+
+  void
+  SolverBase::get_eigenpair (const unsigned int index, 
+                            double &kr,
+                            PETScWrappers::VectorBase &vr) 
+  {
+    AssertThrow (solver_data != 0, ExcSLEPcWrappersUsageError());    
+
+                                    // get converged eigenpair
+    int ierr = EPSGetEigenpair(solver_data->eps, index, 
+                              &kr, PETSC_NULL, vr, PETSC_NULL);
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+  }
+
+
+  void
+  SolverBase::reset ()
+  {
+    AssertThrow (solver_data != 0, ExcSLEPcWrappersUsageError());
+
+                                    // destroy solver object.
+    solver_data.reset ();
+  }
+
+  EPS *
+  SolverBase::get_EPS ()
+  {
+    if( solver_data == 0 )
+      return NULL;
+    return &solver_data->eps;
+  }
+    
+  /* ---------------------- SolverControls ----------------------- */
+  SolverControl &
+  SolverBase::control () const
+  {
+    return solver_control;
+  }
+
+  int
+  SolverBase::convergence_test (EPS                 /*eps*/,
+                                const int           iteration,
+                                const PetscScalar   residual_norm,
+                                EPSConvergedReason *reason,
+                                void               *solver_control_x)
+  {
+    SolverControl &solver_control 
+      = *reinterpret_cast<SolverControl*>(solver_control_x);
+    
+    const SolverControl::State state
+      = solver_control.check (iteration, residual_norm);
+    
+    switch (state)
+      {
+      case ::dealii::SolverControl::iterate:
+       *reason = EPS_CONVERGED_ITERATING;
+       break;
+       
+      case ::dealii::SolverControl::success:
+       *reason = static_cast<EPSConvergedReason>(1);
+       break;
+       
+      case ::dealii::SolverControl::failure:
+       if (solver_control.last_step() > solver_control.max_steps())
+         *reason = EPS_DIVERGED_ITS;
+       break;
+       
+      default:
+       Assert (false, ExcNotImplemented());
+      }
+    
+                                    // return without failure.
+    return 0;
+  }
+
+  /* ---------------------- SolverKrylovSchur ------------------------ */
+  SolverKrylovSchur::SolverKrylovSchur (SolverControl        &cn,
+                                       const MPI_Comm       &mpi_communicator,
+                                       const AdditionalData &data)
+    :
+    SolverBase (cn, mpi_communicator),
+    additional_data (data)
+  {}
+    
+  void
+  SolverKrylovSchur::set_solver_type (EPS &eps) const
+  {
+    int ierr;
+    ierr = EPSSetType (eps, const_cast<char *>(EPSKRYLOVSCHUR));
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+                                    // hand over the absolute
+                                    // tolerance and the maximum
+                                    // number of iteration steps to
+                                    // the SLEPc convergence
+                                    // criterion.
+    ierr = EPSSetTolerances(eps, this->solver_control.tolerance(),
+                           this->solver_control.max_steps());
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+  }
+  
+  /* ---------------------- SolverArnoldi ------------------------ */
+  SolverArnoldi::SolverArnoldi (SolverControl        &cn,
+                               const MPI_Comm       &mpi_communicator,
+                               const AdditionalData &data)
+    :
+    SolverBase (cn, mpi_communicator),
+    additional_data (data)
+  {}
+    
+  void
+  SolverArnoldi::set_solver_type (EPS &eps) const
+  {
+    int ierr;
+    ierr = EPSSetType (eps, const_cast<char *>(EPSARNOLDI));
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+                                    // hand over the absolute
+                                    // tolerance and the maximum
+                                    // number of iteration steps to
+                                    // the SLEPc convergence
+                                    // criterion.
+    ierr = EPSSetTolerances(eps, this->solver_control.tolerance(),
+                           this->solver_control.max_steps());
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+  }
+
+  /* ---------------------- Lanczos ------------------------ */
+
+  SolverLanczos::SolverLanczos (SolverControl        &cn,
+                               const MPI_Comm       &mpi_communicator,
+                               const AdditionalData &data)
+    :
+    SolverBase (cn, mpi_communicator),
+    additional_data (data)
+  {}
+    
+  void
+  SolverLanczos::set_solver_type (EPS &eps) const
+  {
+    int ierr;
+    ierr = EPSSetType (eps, const_cast<char *>(EPSLANCZOS));
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+                                    // hand over the absolute
+                                    // tolerance and the maximum
+                                    // number of iteration steps to
+                                    // the SLEPc convergence
+                                    // criterion.
+    ierr = EPSSetTolerances(eps, this->solver_control.tolerance(),
+                           this->solver_control.max_steps());
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+  }
+
+}
+
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#else
+// On gcc2.95 on Alpha OSF1, the native assembler does not like empty
+// files, so provide some dummy code
+namespace { void dummy () {} }
+#endif // DEAL_II_USE_SLEPC
+
diff --git a/deal.II/lac/source/slepc_spectral_transformation.cc b/deal.II/lac/source/slepc_spectral_transformation.cc
new file mode 100644 (file)
index 0000000..6b4839c
--- /dev/null
@@ -0,0 +1,154 @@
+
+#include <lac/petsc_matrix_base.h>
+#include <lac/petsc_vector_base.h>
+#include <lac/petsc_vector.h>
+#include <lac/slepc_solver.h>
+#include <lac/slepc_spectral_transformation.h>
+
+#include <cmath>
+#include <vector>
+
+#ifdef DEAL_II_USE_SLEPC
+
+#if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2)
+#  include <petscsles.h>
+#endif
+#include <petscversion.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace SLEPcWrappers
+{
+  TransformationBase::TransformationData::~TransformationData ()
+  {
+  }
+
+  TransformationBase::TransformationBase ()
+  {
+  }
+
+  TransformationBase::~TransformationBase ()
+  {
+  }
+
+  void TransformationBase::set_context (EPS &eps)
+  {
+    AssertThrow (transformation_data == 0, 
+                SolverBase::ExcSLEPcWrappersUsageError());
+    transformation_data.reset (new TransformationData());
+    
+    int ierr = EPSGetST(eps, &transformation_data->st);
+    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
+    
+    set_transformation_type(transformation_data->st);
+  }
+  
+  /* ------------------- TransformationShift --------------------- */
+  TransformationShift::AdditionalData::
+  AdditionalData (const double shift_parameter)
+                  :
+                  shift_parameter (shift_parameter)
+  {}
+  
+  TransformationShift::TransformationShift (const AdditionalData &data)
+    :
+    additional_data (data)
+  {}
+
+  void
+  TransformationShift::set_transformation_type (ST &st) const
+  {
+    int ierr;
+    ierr = STSetType (st, const_cast<char *>(STSHIFT));
+    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
+
+    ierr = STSetShift (st, additional_data.shift_parameter);
+    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
+  }
+
+  /* ---------------- TransformationShiftInvert ------------------ */
+  TransformationShiftInvert::AdditionalData::
+  AdditionalData (const double shift_parameter)
+                  :
+                  shift_parameter (shift_parameter)
+  {}
+  
+  TransformationShiftInvert::TransformationShiftInvert (const AdditionalData &data)
+    :
+    additional_data (data)
+  {}
+
+  void
+  TransformationShiftInvert::set_transformation_type (ST &st) const
+  {
+    int ierr;
+    ierr = STSetType (st, const_cast<char *>(STSINV));
+    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
+
+    ierr = STSetShift (st, additional_data.shift_parameter);
+    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
+  }
+
+  /* --------------- TransformationSpectrumFolding ----------------- */
+  TransformationSpectrumFolding::AdditionalData::
+  AdditionalData (const double shift_parameter)
+                  :
+                  shift_parameter (shift_parameter)
+  {}
+  
+  TransformationSpectrumFolding::TransformationSpectrumFolding (const AdditionalData &data)
+    :
+    additional_data (data)
+  {}
+
+
+  void
+  TransformationSpectrumFolding::set_transformation_type (ST &st) const
+  {
+    int ierr;
+    ierr = STSetType (st, const_cast<char *>(STFOLD));
+    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
+
+    ierr = STSetShift (st, additional_data.shift_parameter);
+    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
+  }
+  
+  /* ------------------- TransformationCayley --------------------- */
+  TransformationCayley::AdditionalData::
+  AdditionalData (const double shift_parameter, 
+                 const double antishift_parameter)
+                  :
+                  shift_parameter (shift_parameter),
+                 antishift_parameter (antishift_parameter)
+  {
+  }
+  
+  TransformationCayley::TransformationCayley (const double shift, 
+                                             const double antishift)
+    :
+    additional_data (shift, antishift)
+  {}
+
+  void
+  TransformationCayley::set_transformation_type (ST &st) const
+  {
+    int ierr = STSetType (st, const_cast<char *>(STCAYLEY));
+    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
+
+    ierr = STSetShift (st, additional_data.shift_parameter);
+    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
+
+    ierr = STCayleySetAntishift (st, additional_data.antishift_parameter);
+    AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
+  }
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#else
+// On gcc2.95 on Alpha OSF1, the native assembler does not like empty
+// files, so provide some dummy code
+namespace { void dummy () {} }
+#endif // DEAL_II_USE_SLEPC
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.