]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New to the SLEPcWrappers: A shortcut for the (unlikely?) case where the mass matrix...
authorToby D. Young <tyoung@ippt.pan.pl>
Mon, 13 Jul 2009 13:23:50 +0000 (13:23 +0000)
committerToby D. Young <tyoung@ippt.pan.pl>
Mon, 13 Jul 2009 13:23:50 +0000 (13:23 +0000)
git-svn-id: https://svn.dealii.org/trunk@19070 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/lac/include/lac/slepc_solver.h
deal.II/lac/source/slepc_solver.cc

index 142dee9b73433f70c8c34e2d5a23499a6831a0d9..3ff6913670d1a16dd92893692b64771fe22aa79a 100644 (file)
@@ -180,6 +180,18 @@ inconvenience this causes.
 <a name="lac"></a>
 <h3>lac</h3>
 
+<ol>
+  <li>
+  <p>
+  New: Based on work by Francisco Alvaro, the existing SLEPcWrappers now 
+  have a handle on the generalized eigenvalue problem where B=I.
+  <br>
+  (Toby D. Young 2009/06/25)
+  </p>
+  </li>
+</ol>
+
+
 <ol>
   <li>
   <p>
@@ -214,14 +226,16 @@ inconvenience this causes.
 <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. 
+  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>
index fef1bbeffdb7a36a2775d4d4798bcb883926d5e1..ac5cc49ea741fe4c0cad54f5cdf4eafdd62c22fa 100755 (executable)
@@ -19,8 +19,9 @@ 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.
+ * solver. On the other hand, note that: the AdditionalData structure
+ * is a dummy structure and is there for backward/forward
+ * compatibility.
  *
  * SLEPcWrappers can be implemented in application codes in the
  * following way:
@@ -30,13 +31,14 @@ DEAL_II_NAMESPACE_OPEN
                         mpi_communicator);
   system.solve (A, B, lambda, x, n_eigenvectors);
  @endverbatim
- * for the generalized eigenvalue problem $Ax=B\lambda x$.
+ * for the generalized eigenvalue problem $Ax=B\lambda x$. See also
+ * @ref step 36 "step-36" for a hands-on example.
  *
  * 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
+ * rather than just one. This freedom is intended for use of the
+ * SLEPcWrappers that require a greater handle on the eigenvalue
  * problem solver context. See also:
  @verbatim
   template <typename OutputVector>
@@ -45,8 +47,8 @@ DEAL_II_NAMESPACE_OPEN
                     const PETScWrappers::MatrixBase &B,
                     std::vector<double>             &kr,
                     std::vector<OutputVector>       &vr,
-                    const unsigned int               n_eigenvectors
- )
+                    const unsigned int               n_eigenvectors
+  {<code>code...</code>}
  @endverbatim
  * as an example on how to do this. 
  *
@@ -92,13 +94,13 @@ namespace SLEPcWrappers
 
                                    /**
                                     * 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
+                                    * eigensystem $Ax=\lambda x$. 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
@@ -113,9 +115,9 @@ namespace SLEPcWrappers
                                     * reset.
                                     * 
                                     * Note that the number of
-                                    * converged eigenstates can be
+                                    * converged eigenvectors can be
                                     * larger than the number of
-                                    * eigenstates requested; this is
+                                    * eigenvectors requested; this is
                                     * due to a round off error
                                     * (success) of the eigenvalue
                                     * solver context. If this is found
@@ -124,6 +126,18 @@ namespace SLEPcWrappers
                                     * requested but handle that it may
                                     * be more by ignoring any extras.
                                    */
+      template <typename OutputVector>
+       void
+       solve (const PETScWrappers::MatrixBase &A,
+              std::vector<double>             &kr,
+              std::vector<OutputVector>       &vr,
+              const unsigned int              n_eigenvectors);
+
+                                   /** 
+                                   * Same as above, but here is a
+                                    * composite method for solving the
+                                    * system $A x=\lambda B x$.
+                                    */
       template <typename OutputVector>
        void
        solve (const PETScWrappers::MatrixBase &A,
@@ -134,23 +148,32 @@ namespace SLEPcWrappers
       
                                    /**
                                     * Initialize solver for the linear
-                                    * system $Ax=\lambda
-                                    * Bx$. (required before calling
-                                    * solve)
+                                    * system $Ax=\lambda x$. (Note:
+                                    * this is required before calling
+                                    * solve ())
+                                    */
+      void
+       set_matrices (const PETScWrappers::MatrixBase &A);   
+
+                                   /** 
+                                   * Same as above, but here is a
+                                    * composite method for solving the
+                                    * system $A x=\lambda B x$.
                                     */
       void
        set_matrices (const PETScWrappers::MatrixBase &A,
                      const PETScWrappers::MatrixBase &B);   
 
                                    /**
-                                    * Set the initial vector for the solver.
+                                    * 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
+                                   * to be used.
                                     */
       void
        set_transformation (SLEPcWrappers::TransformationBase &trans);
@@ -170,9 +193,10 @@ namespace SLEPcWrappers
                                     * Solve the linear system for
                                    * n_eigenvectors
                                    * eigenstates. Parameter
-                                   * n_converged contains the actual
-                                   * number of eigenstates that have
-                                   * .  converged; this can be both
+                                   * <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.
@@ -211,7 +235,7 @@ namespace SLEPcWrappers
                                     * Access to object that controls
                                     * convergence.
                                     */
-      SolverControl & control() const;
+      SolverControl &control() const;
       
                                    /**
                                     * Exceptions.
@@ -458,44 +482,70 @@ namespace SLEPcWrappers
 
 
   // --------------------------- inline and template functions -----------                                                                       
-
-  /**                                                                                                           
-   * This is declared here to make it                                                                           
-   * possible to take a std::vector                                                                             
-   * of different PETScWrappers vector                                                                          
-   * types                                                                                                      
+  /**
+   * 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,
+                       std::vector<double>             &kr,
+                       std::vector<OutputVector>       &vr,
+                       const unsigned int               n_eigenvectors = 1)
+    {
+      unsigned int n_converged;
+      
+      set_matrices(A);
+      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]);
+        }
+    }
+  
+  
   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)
+                       const unsigned int               n_eigenvectors = 1)
     {
       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
index f9042bece1b7198bc77fa361da338e4bc6515bfa..943cdc73b9dac97d7da982fa154cc6f9a2e893da 100644 (file)
@@ -44,7 +44,14 @@ namespace SLEPcWrappers
     if( solver_data != 0 )
       solver_data.reset ();
   }
-  
+
+  void
+  SolverBase::set_matrices (const PETScWrappers::MatrixBase &A)
+  {
+    opA = &A;
+    opB = NULL;
+  }
+
   void
   SolverBase::set_matrices (const PETScWrappers::MatrixBase &A,
                            const PETScWrappers::MatrixBase &B)
@@ -84,19 +91,22 @@ namespace SLEPcWrappers
     ierr = EPSCreate (mpi_communicator, &solver_data->eps);
     AssertThrow (ierr == 0, ExcSLEPcError(ierr));
 
-    AssertThrow (opA && opB, ExcSLEPcWrappersUsageError());
-    ierr = EPSSetOperators (solver_data->eps, *opA, *opB);
+    AssertThrow (opA, ExcSLEPcWrappersUsageError());
+    if (opB)
+      ierr = EPSSetOperators (solver_data->eps, *opA, *opB);
+    else
+      ierr = EPSSetOperators (solver_data->eps, *opA, PETSC_NULL);
     AssertThrow (ierr == 0, ExcSLEPcError(ierr));
-
-    if( ini_vec && ini_vec->size() != 0 
+    
+    if (ini_vec && ini_vec->size() != 0
       {
        ierr = EPSSetInitialVector(solver_data->eps, *ini_vec);
        AssertThrow (ierr == 0, ExcSLEPcError(ierr));
       }
     
-    if( transform )
+    if (transform)
       transform->set_context(solver_data->eps);
-
+    
                                     // set runtime options.
     set_solver_type (solver_data->eps);
 
@@ -118,7 +128,8 @@ namespace SLEPcWrappers
 
                                     // get number of converged
                                     // eigenstates
-    ierr = EPSGetConverged (solver_data->eps, reinterpret_cast<int *>(n_converged));
+    ierr = EPSGetConverged (solver_data->eps, 
+                           reinterpret_cast<int *>(n_converged));
     AssertThrow (ierr == 0, ExcSLEPcError(ierr));
   }
 
@@ -277,8 +288,6 @@ namespace SLEPcWrappers
 
 }
 
-
-
 DEAL_II_NAMESPACE_CLOSE
 
 #else

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.