]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
New function set_target_eigenvalue and minor API tweaks.
authoryoung <young@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Apr 2013 11:08:30 +0000 (11:08 +0000)
committeryoung <young@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Apr 2013 11:08:30 +0000 (11:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@29370 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/lac/slepc_solver.h
deal.II/source/lac/slepc_solver.cc

index 79856231147c24149985f8b23c577f7e7febd031..d28d46c02801b5b2c1d7d2861c5fdcbb5e807d7b 100644 (file)
@@ -63,7 +63,7 @@ DEAL_II_NAMESPACE_OPEN
  * SLEPc solvers before calling <code>solve()</code>. For example, if
  * the matrices of the general eigenspectrum problem are not hermitian
  * and the lower eigenvalues are wanted only, the following code can
- * be implemented:
+ * be implemented before calling <code>solve()</code>:
  * @code
  *  system.set_problem_type (EPS_NHEP);
  *  system.set_which_eigenpairs (EPS_SMALLEST_REAL);
@@ -102,7 +102,6 @@ DEAL_II_NAMESPACE_OPEN
  * @note Various tweaks and enhancments contributed by Eloy Romero and
  * Jose E. Roman 2009, 2010.
  */
-
 namespace SLEPcWrappers
 {
 
@@ -142,9 +141,9 @@ namespace SLEPcWrappers
      * Set up matrices for solving; Actually solve the system; Gather
      * the solution(s); and reset.
      *
-     * Note that the number of converged eigenvectors can be larger
-     * than the number of eigenvectors requested; this is due to a
-     * round off error (success) of the eigenproblem solver
+     * @note Note that the number of converged eigenvectors can be
+     * larger than the number of eigenvectors requested; this is due
+     * to a round off error (success) of the eigenproblem 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 than specified by ignoring any extras. By default
@@ -169,21 +168,6 @@ namespace SLEPcWrappers
            std::vector<OutputVector>       &vr,
            const unsigned int               n_eigenvectors = 1);
 
-    /**
-     * Initialize solver for the linear system $Ax=\lambda x$. (Note:
-     * this is required before calling solve ())
-     */
-    void
-    set_matrices (const PETScWrappers::MatrixBase &A);
-
-    /**
-     * Same as above, but here initialize solver for the linear system
-     * $A x=\lambda B x$.
-     */
-    void
-    set_matrices (const PETScWrappers::MatrixBase &A,
-                  const PETScWrappers::MatrixBase &B);
-
     /**
      * Set the initial vector for the solver.
      */
@@ -197,6 +181,13 @@ namespace SLEPcWrappers
     void
     set_transformation (SLEPcWrappers::TransformationBase &this_transformation);
 
+    /**
+     * Set target eigenvalues in the spectrum to be computed. By
+     * default, no target is set.
+     */
+    void
+    set_target_eigenvalue (const double &this_target);
+
     /**
      * Indicate which part of the spectrum is to be computed. By
      * default largest magnitude eigenvalues are computed.
@@ -217,52 +208,13 @@ namespace SLEPcWrappers
     void
     set_problem_type (EPSProblemType set_problem);
 
-    /**
-     * Solve the linear system for n_eigenvectors
-     * eigenstates. Parameter <code>n_converged</code> 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, $\text{index}\,\in\,0\hdots
-     * \text{n\_converged}-1$.
-     */
-    void
-    get_eigenpair (const unsigned int index,
-#ifndef PETSC_USE_COMPLEX
-                   double                    &kr,
-#else
-                   std::complex<double>      &kr,
-#endif
-                   PETScWrappers::VectorBase &vr);
-
-    /**
-     * Reset the solver, and return memory for eigenvectors
-     */
-    void
-    reset ();
-
-    /**
-     * Retrieve the SLEPc solver object that is internally used.
-     */
-    EPS *get_eps ();
-
     /**
      * Take the information provided from SLEPc and checks it against
      * deal.II's own SolverControl objects to see if convergence has
      * been reached.
      */
-    void get_solver_state (const SolverControl::State state);
-
-    /**
-     * Access to the object that controls convergence.
-     */
-    SolverControl &control () const;
+    void 
+    get_solver_state (const SolverControl::State state);
 
     /**
      * Exception. Standard exception.
@@ -285,6 +237,11 @@ namespace SLEPcWrappers
                     << "    The number of converged eigenvectors is " << arg1
                     << " but " << arg2 << " were requested. ");
 
+    /**
+     * Access to the object that controls convergence.
+     */
+    SolverControl &control () const;
+
   protected:
 
     /**
@@ -305,6 +262,61 @@ namespace SLEPcWrappers
      */
     virtual void set_solver_type (EPS &eps) const = 0;
 
+    /**
+     * Reset the solver, and return memory for eigenvectors
+     */
+    void
+    reset ();
+
+    /**
+     * Retrieve the SLEPc solver object that is internally used.
+     */
+    EPS *get_eps ();
+
+    /**
+     * Solve the linear system for n_eigenvectors
+     * eigenstates. Parameter <code>n_converged</code> 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, $\text{index}\,\in\,0\hdots
+     * \text{n\_converged}-1$.
+     */
+    void
+    get_eigenpair (const unsigned int index,
+#ifndef PETSC_USE_COMPLEX
+                   double                    &kr,
+#else
+                   std::complex<double>      &kr,
+#endif
+                   PETScWrappers::VectorBase &vr);
+
+    /**
+     * Initialize solver for the linear system $Ax=\lambda x$. (Note:
+     * this is required before calling solve ())
+     */
+    void
+    set_matrices (const PETScWrappers::MatrixBase &A);
+
+    /**
+     * Same as above, but here initialize solver for the linear system
+     * $A x=\lambda B x$.
+     */
+    void
+    set_matrices (const PETScWrappers::MatrixBase &A,
+                  const PETScWrappers::MatrixBase &B);
+
+    /**
+     * Target eigenvalue to solve for.
+     */
+    PetscScalar target_eigenvalue;
+
     /**
      * Which portion of the spectrum to solve from.
      */
@@ -315,6 +327,8 @@ namespace SLEPcWrappers
      */
     EPSProblemType set_problem;
 
+  private:
+
     /**
      * The matrix $A$ of the generalized eigenvalue problem
      * $Ax=B\lambda x$, or the standard eigenvalue problem $Ax=\lambda
@@ -339,8 +353,6 @@ namespace SLEPcWrappers
      */
     SLEPcWrappers::TransformationBase *transformation;
 
-  private:
-
     /**
      * A function that can be used in SLEPc as a callback to check on
      * convergence.
@@ -718,7 +730,7 @@ namespace SLEPcWrappers
     unsigned int n_converged = 0;
 
     // Set the matrices of the problem
-    set_matrices (A);
+    set_matrices (A); 
 
     // and solve
     solve (n_eigenvectors, &n_converged);
@@ -738,11 +750,11 @@ namespace SLEPcWrappers
 
   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)
+    SolverBase::solve (const PETScWrappers::MatrixBase &A,
+                      const PETScWrappers::MatrixBase &B,
+                      std::vector<double>             &kr,
+                      std::vector<OutputVector>       &vr,
+                      const unsigned int               n_eigenvectors)
   {
     // Panic if no eigenpairs are wanted.
     AssertThrow (n_eigenvectors != 0, ExcSLEPcWrappersUsageError());
index c76b6c43319e58f549643245b934ceeb4ecb6694..ded535c6bcf72acd45bb76dc659356ec50f06a43 100644 (file)
@@ -49,6 +49,7 @@ namespace SLEPcWrappers
     :
     solver_control (cn),
     mpi_communicator (mpi_communicator),
+    target_eigenvalue (PETSC_NULL),
     set_which (EPS_LARGEST_MAGNITUDE),
     set_problem (EPS_NHEP),
     opA (NULL),
@@ -89,6 +90,12 @@ namespace SLEPcWrappers
     transformation = &this_transformation;
   }
 
+  void
+  SolverBase::set_target_eigenvalue (const double &this_target)
+  {
+    target_eigenvalue = this_target;
+  }
+
   void
   SolverBase::set_which_eigenpairs (const EPSWhich eps_which)
   {
@@ -147,6 +154,10 @@ namespace SLEPcWrappers
     if (transformation)
       transformation->set_context (solver_data->eps);
 
+    // set target eigenvalues to solve for
+    ierr = EPSSetTarget (solver_data->eps, target_eigenvalue);
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+    
     // set which portion of the eigenspectrum to solve for
     ierr = EPSSetWhichEigenpairs (solver_data->eps, set_which);
     AssertThrow (ierr == 0, ExcSLEPcError(ierr));

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.