]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Get imaginary parts of eigenpairs.
authoryoung <young@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 10 May 2013 11:51:21 +0000 (11:51 +0000)
committeryoung <young@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 10 May 2013 11:51:21 +0000 (11:51 +0000)
git-svn-id: https://svn.dealii.org/trunk@29491 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 402f4ebf69d4cfd42d516306f0867b98ed776658..fb74ab55f774b267a021dacc343511d1983c0c7a 100644 (file)
@@ -84,7 +84,7 @@ 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_eigenpairs)
   { ... }
  @endcode
  * as an example on how to do this.
@@ -134,11 +134,11 @@ namespace SLEPcWrappers
      * 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, the fairly
-     * standard calling sequence executed here is used: Initialise;
-     * Set up matrices for solving; Actually solve the system; Gather
-     * the solution(s); and reset.
+     * eigenvectors, at least twice the memory size of
+     * <tt>eigenvectors</tt> is being used (and can be more). To avoid
+     * doing this, the fairly standard calling sequence executed here
+     * is used: Initialise; Set up matrices for solving; Actually
+     * solve the system; Gather the solution(s); and reset.
      *
      * @note Note that the number of converged eigenvectors can be
      * larger than the number of eigenvectors requested; this is due
@@ -151,8 +151,8 @@ namespace SLEPcWrappers
     template <typename OutputVector>
     void
     solve (const PETScWrappers::MatrixBase &A,
-           std::vector<double>             &r_eigenvalues,
-           std::vector<OutputVector>       &r_eigenvectors = std::vector<OutputVector> (),
+           std::vector<double>             &eigenvalues,
+           std::vector<OutputVector>       &eigenvectors,
            const unsigned int               n_eigenvectors = 1);
 
     /**
@@ -163,8 +163,23 @@ namespace SLEPcWrappers
     void
     solve (const PETScWrappers::MatrixBase &A,
            const PETScWrappers::MatrixBase &B,
-           std::vector<double>             &r_eigenvalues,
-           std::vector<OutputVector>       &r_eigenvectors = std::vector<OutputVector> (),
+           std::vector<double>             &eigenvalues,
+           std::vector<OutputVector>       &eigenvectors,
+           const unsigned int               n_eigenvectors = 1);
+
+    /**
+     * Same as above, but here a composite method for solving the
+     * system $A x=\lambda B x$ with real and imaginary parts to the
+     * eigenvalues/eigenvectors.
+     */
+    template <typename OutputVector>
+    void
+    solve (const PETScWrappers::MatrixBase &A,
+           const PETScWrappers::MatrixBase &B,
+           std::vector<double>             &real_eigenvalues,
+           std::vector<double>             &imag_eigenvalues,
+           std::vector<OutputVector>       &real_eigenvectors,
+           std::vector<OutputVector>       &imag_eigenvectors,
            const unsigned int               n_eigenvectors = 1);
 
     /**
@@ -283,14 +298,26 @@ namespace SLEPcWrappers
     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
+     * Access the real parts of 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,
-                   double                    &kr,
-                   PETScWrappers::VectorBase &vr);
+                   double                    &eigenvalues,
+                   PETScWrappers::VectorBase &eigenvectors);
+
+    /**
+     * Access the real and imaginary parts of 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,
+                   double                    &real_eigenvalues,
+                   double                    &imag_eigenvalues,
+                   PETScWrappers::VectorBase &real_eigenvectors,
+                   PETScWrappers::VectorBase &imag_eigenvectors);
 
     /**
      * Initialize solver for the linear system $Ax=\lambda x$. (Note:
@@ -711,74 +738,126 @@ namespace SLEPcWrappers
    * This is declared here to make it possible to take a std::vector
    * of different PETScWrappers vector types
    */
+  // todo: The logic of these functions can be simplified without breaking backward compatibility...
 
   template <typename OutputVector>
-  void
-  SolverBase::solve (const PETScWrappers::MatrixBase &A,
-                     std::vector<double>             &kr,
-                     std::vector<OutputVector>       &vr,
-                     const unsigned int               n_eigenvectors)
-  {
-    // Panic if the number of eigenpairs wanted is out of bounds.
-    AssertThrow ((n_eigenvectors > 0) && (n_eigenvectors <= A.m ()), 
-                ExcSLEPcWrappersUsageError());
-
-    // Set the matrices of the problem
-    set_matrices (A); 
-
-    // and solve
-    unsigned int n_converged = 0;
-    solve (n_eigenvectors, &n_converged);
-
-    if (n_converged > n_eigenvectors)
-      n_converged = n_eigenvectors;
-    AssertThrow (n_converged == n_eigenvectors,
-                 ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenvectors));
-
-    AssertThrow (vr.size() != 0, 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]);
-  }
-
+    void
+    SolverBase::solve (const PETScWrappers::MatrixBase &A,
+                      std::vector<double>             &eigenvalues,
+                      std::vector<OutputVector>       &eigenvectors,
+                      const unsigned int               n_eigenvectors)
+    {
+      // Panic if the number of eigenpairs wanted is out of bounds.
+      AssertThrow ((n_eigenvectors > 0) && (n_eigenvectors <= A.m ()), 
+                  ExcSLEPcWrappersUsageError());
+      
+      // Set the matrices of the problem
+      set_matrices (A);  
+      
+      // and solve
+      unsigned int n_converged = 0; 
+      solve (n_eigenvectors, &n_converged); 
+      
+      if (n_converged > n_eigenvectors) 
+       n_converged = n_eigenvectors; 
+      AssertThrow (n_converged == n_eigenvectors, 
+                  ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenvectors)); 
+      
+      AssertThrow (eigenvectors.size() != 0, ExcSLEPcWrappersUsageError()); 
+      eigenvectors.resize (n_converged, eigenvectors.front()); 
+      eigenvalues.resize (n_converged); 
+      
+      for (unsigned int index=0; index<n_converged; ++index) 
+       get_eigenpair (index, eigenvalues[index], eigenvectors[index]); 
+    }
+  
   template <typename OutputVector>
-  void
+    void
     SolverBase::solve (const PETScWrappers::MatrixBase &A,
                       const PETScWrappers::MatrixBase &B,
-                      std::vector<double>             &kr,
-                      std::vector<OutputVector>       &vr,
+                      std::vector<double>             &eigenvalues,
+                      std::vector<OutputVector>       &eigenvectors,
                       const unsigned int               n_eigenvectors)
-  {
-    // Guard against incompatible matrix sizes:
-    AssertThrow (A.m() == B.m (), ExcDimensionMismatch(A.m(), B.m()));
-    AssertThrow (A.n() == B.n (), ExcDimensionMismatch(A.n(), B.n()));
+    { 
+      // Guard against incompatible matrix sizes:
+      AssertThrow (A.m() == B.m (), ExcDimensionMismatch(A.m(), B.m()));
+      AssertThrow (A.n() == B.n (), ExcDimensionMismatch(A.n(), B.n()));
+      
+      // Panic if the number of eigenpairs wanted is out of bounds.
+      AssertThrow ((n_eigenvectors > 0) && (n_eigenvectors <= A.m ()), 
+                  ExcSLEPcWrappersUsageError());
+      
+      // Set the matrices of the problem
+      set_matrices (A, B);
+      
+      // and solve
+      unsigned int n_converged = 0;
+      solve (n_eigenvectors, &n_converged);
+      
+      if (n_converged >= n_eigenvectors)
+       n_converged = n_eigenvectors;
+      
+      AssertThrow (n_converged == n_eigenvectors,
+                  ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenvectors));
+      AssertThrow (eigenvectors.size() != 0, ExcSLEPcWrappersUsageError());
+      
+      eigenvectors.resize (n_converged, eigenvectors.front());
+      eigenvalues.resize (n_converged);
+      
+      for (unsigned int index=0; index<n_converged; ++index)
+       get_eigenpair (index, eigenvalues[index], eigenvectors[index]);
+    }
 
-    // Panic if the number of eigenpairs wanted is out of bounds.
-    AssertThrow ((n_eigenvectors > 0) && (n_eigenvectors <= A.m ()), 
-                ExcSLEPcWrappersUsageError());
-
-    // Set the matrices of the problem
-    set_matrices (A, B);
-
-    // and solve
-    unsigned int n_converged = 0;
-    solve (n_eigenvectors, &n_converged);
-
-    if (n_converged >= n_eigenvectors)
-      n_converged = n_eigenvectors;
-
-    AssertThrow (n_converged == n_eigenvectors,
-                 ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenvectors));
-    AssertThrow (vr.size() != 0, ExcSLEPcWrappersUsageError());
-
-    vr.resize (n_converged, vr.front());
-    kr.resize (n_converged);
+  template <typename OutputVector>
+    void
+    SolverBase::solve (const PETScWrappers::MatrixBase &A,
+                      const PETScWrappers::MatrixBase &B,
+                      std::vector<double>             &real_eigenvalues,
+                      std::vector<double>             &imag_eigenvalues,
+                      std::vector<OutputVector>       &real_eigenvectors,
+                      std::vector<OutputVector>       &imag_eigenvectors,
+                      const unsigned int               n_eigenvectors)
+    { 
+      // Guard against incompatible matrix sizes:
+      AssertThrow (A.m() == B.m (), ExcDimensionMismatch(A.m(), B.m()));
+      AssertThrow (A.n() == B.n (), ExcDimensionMismatch(A.n(), B.n()));
+
+      // and incompatible eigenvalue/eigenvector sizes
+      AssertThrow (real_eigenvalues.size() == imag_eigenvalues.size(),  
+                          ExcDimensionMismatch(real_eigenvalues.size(), imag_eigenvalues.size())); 
+      AssertThrow (real_eigenvectors.size() == imag_eigenvectors.n (),  
+                          ExcDimensionMismatch(real_eigenvectors.size(), imag_eigenvectors.size())); 
+
+      // Panic if the number of eigenpairs wanted is out of bounds.
+      AssertThrow ((n_eigenvectors > 0) && (n_eigenvectors <= A.m ()), 
+                  ExcSLEPcWrappersUsageError());
+      
+      // Set the matrices of the problem
+      set_matrices (A, B);
+      
+      // and solve
+      unsigned int n_converged = 0;
+      solve (n_eigenvectors, &n_converged);
+      
+      if (n_converged >= n_eigenvectors)
+       n_converged = n_eigenvectors;
+      
+      AssertThrow (n_converged == n_eigenvectors,
+                  ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenvectors));
+      AssertThrow ((real_eigenvectors.size()!=0) && (imag_eigenvectors.size()!=0), 
+                  ExcSLEPcWrappersUsageError());
+      
+      real_eigenvectors.resize (n_converged, real_eigenvectors.front());
+      imag_eigenvectors.resize (n_converged, imag_eigenvectors.front());
+      real_eigenvalues.resize (n_converged);
+      imag_eigenvalues.resize (n_converged);
+      
+      for (unsigned int index=0; index<n_converged; ++index)
+       get_eigenpair (index, 
+                      real_eigenvalues[index], imag_eigenvalues[index], 
+                      real_eigenvectors[index], imag_eigenvectors[index]);
+    }
 
-    for (unsigned int index=0; index<n_converged; ++index)
-      get_eigenpair (index, kr[index], vr[index]);
-  }
 }
 
 DEAL_II_NAMESPACE_CLOSE
index 9d5877f2aaa7eb8b1eb79f0482c912162952bb37..906cfde2e49b57303f9c74d14e35f5be055bbd10 100644 (file)
@@ -231,6 +231,22 @@ namespace SLEPcWrappers
     AssertThrow (ierr == 0, ExcSLEPcError(ierr));
   }
 
+  void
+  SolverBase::get_eigenpair (const unsigned int         index,
+                             double                    &real_eigenvalues,
+                             double                    &imag_eigenvalues,
+                             PETScWrappers::VectorBase &real_eigenvectors,
+                             PETScWrappers::VectorBase &imag_eigenvectors)
+  {
+    AssertThrow (solver_data.get() != 0, ExcSLEPcWrappersUsageError());
+    
+    // get converged eigenpair
+    int ierr = EPSGetEigenpair (solver_data->eps, index,
+                               &real_eigenvalues, &imag_eigenvalues, 
+                               real_eigenvectors, imag_eigenvectors);
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+  }
+
   void
   SolverBase::reset ()
   {

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.