]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Improve SLEPcWrappers. Fix usage of SolverControl.
authoryoung <young@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Jan 2013 18:38:44 +0000 (18:38 +0000)
committeryoung <young@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Jan 2013 18:38:44 +0000 (18:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@28132 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 4e1b09fbc19298a54a92c6c51a632d0522fb5959..7e6593810c082e95fa014da5f646b91180cea296 100644 (file)
@@ -191,11 +191,11 @@ namespace SLEPcWrappers
     set_which_eigenpairs (EPSWhich set_which);
 
     /**
-     * Solve the linear system for <code>n_eigenvectors</code>
+     * 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 <code>n_eigenvectors</code>, depending
-     * on the SLEPc eigensolver used.
+     * 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);
@@ -218,13 +218,19 @@ namespace SLEPcWrappers
      * Reset the solver, and return memory for eigenvectors
      */
     void
-    reset();
+    reset ();
 
     /**
      * Retrieve the SLEPc solver object that is internally used.
      */
-    EPS *
-    get_EPS ();
+    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.
@@ -232,15 +238,21 @@ namespace SLEPcWrappers
     SolverControl &control () const;
 
     /**
-     * Exceptions.
+     * Exception. Standard exception.
      */
     DeclException0 (ExcSLEPcWrappersUsageError);
 
+    /**
+     * Exception. SLEPc error with error number.
+     */
     DeclException1 (ExcSLEPcError,
                     int,
                     << "    An error with error number " << arg1
                     << " occurred while calling a SLEPc function");
 
+    /**
+     * Exception. Convergence failure on the number of eigenvectors.
+     */
     DeclException2 (ExcSLEPcEigenvectorConvergenceMismatchError,
                     int, int,
                     << "    The number of converged eigenvectors is " << arg1
@@ -298,23 +310,20 @@ namespace SLEPcWrappers
   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.
+     * A function that can be used in SLEPc as a callback to check on
+     * convergence. 
+     *
+     * @note This function is redundant.
      */
-    static
+    static 
     int
-    convergence_test (EPS                 eps,
-#ifdef PETSC_USE_64BIT_INDICES
-                      const PetscInt      iteration,
-#else
-                      const int           iteration,
-#endif
-                      const PetscReal     residual_norm,
-                      EPSConvergedReason *reason,
-                      void               *solver_control);
-
+    convergence_test (EPS          eps,
+                     PetscScalar  kr,
+                     PetscScalar  ki,
+                     PetscReal    residual_norm,
+                     PetscReal   *estimated_error,
+                      void        *solver_control);
+    
     /**
      * Objects of this type are explicitly created, but are destroyed
      * when the surrounding solver object goes out of scope, or when
@@ -335,6 +344,11 @@ namespace SLEPcWrappers
        * Objects for Eigenvalue Problem Solver.
        */
       EPS eps;
+
+      /**
+       * Convergence.
+       */
+      EPSConvergedReason reason;
     };
 
     /**
@@ -349,6 +363,7 @@ namespace SLEPcWrappers
    * complex.
    *
    * @ingroup SLEPcWrappers
+   *
    * @author Toby D. Young 2008
    */
   class SolverKrylovSchur : public SolverBase
@@ -391,6 +406,7 @@ namespace SLEPcWrappers
    * solver. Usage: All spectrum, all problem types, complex.
    *
    * @ingroup SLEPcWrappers
+   *
    * @author Toby D. Young 2008, 2011
    */
   class SolverArnoldi : public SolverBase
@@ -436,7 +452,6 @@ namespace SLEPcWrappers
      * and sets the type of solver that is appropriate for this class.
      */
     virtual void set_solver_type (EPS &eps) const;
-
   };
 
   /**
@@ -444,6 +459,7 @@ namespace SLEPcWrappers
    * solver. Usage: All spectrum, all problem types, complex.
    *
    * @ingroup SLEPcWrappers
+   *
    * @author Toby D. Young 2009
    */
   class SolverLanczos : public SolverBase
@@ -478,7 +494,6 @@ namespace SLEPcWrappers
      * and sets the type of solver that is appropriate for this class.
      */
     virtual void set_solver_type (EPS &eps) const;
-
   };
 
   /**
@@ -486,7 +501,8 @@ namespace SLEPcWrappers
    * solver. Usage: Largest values of spectrum only, all problem
    * types, complex.
    *
-   * @ingroup SLEPcWrappers
+   * @ingroup SLEPcWrappers 
+   *
    * @author Toby D. Young 2010
    */
   class SolverPower : public SolverBase
@@ -521,14 +537,15 @@ namespace SLEPcWrappers
      * 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
-   * Generalized Davidson solver. Usage: All problem types.
+   * Davidson solver. Usage (incomplete/untested): All problem types.
    *
-   * @ingroup SLEPcWrappers @author Toby D. Young 2013
+   * @ingroup SLEPcWrappers 
+   *
+   * @author Toby D. Young 2010
    */
   class SolverGeneralizedDavidson : public SolverBase
   {
@@ -562,14 +579,15 @@ namespace SLEPcWrappers
      * 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
    * Jacobi-Davidson solver. Usage: All problem types.
    *
-   * @ingroup SLEPcWrappers @author Toby D. Young 2013
+   * @ingroup SLEPcWrappers 
+   *
+   * @author Toby D. Young 2013
    */
   class SolverJacobiDavidson : public SolverBase
   {
@@ -603,7 +621,6 @@ namespace SLEPcWrappers
      * and sets the type of solver that is appropriate for this class.
      */
     virtual void set_solver_type (EPS &eps) const;
-
   };
 
   // --------------------------- inline and template functions -----------
@@ -619,21 +636,24 @@ namespace SLEPcWrappers
                      std::vector<OutputVector>       &vr,
                      const unsigned int               n_eigenvectors = 1)
   {
-    unsigned int n_converged;
+    unsigned int n_converged = 0;
 
+    // Set the matrices of the problem 
     set_matrices (A);
-    solve (n_eigenvectors,&n_converged);
 
+    // and solve    
+    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() >= 1, ExcSLEPcWrappersUsageError());
+    AssertThrow (vr.size() != 0, ExcSLEPcWrappersUsageError());
     vr.resize (n_converged, vr.front());
     kr.resize (n_converged);
 
-    for (unsigned int index=0; index < n_converged; ++index)
+    for (unsigned int index=0; index<n_converged; ++index)
       get_eigenpair (index, kr[index], vr[index]);
   }
 
@@ -645,9 +665,13 @@ namespace SLEPcWrappers
                      std::vector<OutputVector>       &vr,
                      const unsigned int               n_eigenvectors = 1)
   {
-    unsigned int n_converged;
+    int ierr;
+    unsigned int n_converged = 0;
 
+    // Set the matrices of the problem 
     set_matrices (A, B);
+
+    // and solve
     solve (n_eigenvectors, &n_converged);
 
     if (n_converged >= n_eigenvectors)
@@ -655,12 +679,12 @@ namespace SLEPcWrappers
 
     AssertThrow (n_converged == n_eigenvectors,
                  ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenvectors));
-    AssertThrow (vr.size() >= 1, ExcSLEPcWrappersUsageError());
+    AssertThrow (vr.size() != 0, ExcSLEPcWrappersUsageError());
 
     vr.resize (n_converged, vr.front());
     kr.resize (n_converged);
 
-    for (unsigned int index=0; index < n_converged; ++index)
+    for (unsigned int index=0; index<n_converged; ++index)
       get_eigenpair (index, kr[index], vr[index]);
   }
 }
index 7f73e250e87c3c0c24697265e1ebfc05745a68ed..926a76f1ade9e44791d705e418f0e61b9ab8c03d 100644 (file)
@@ -50,7 +50,8 @@ namespace SLEPcWrappers
     solver_control (cn),
     mpi_communicator (mpi_communicator),
     set_which (EPS_LARGEST_MAGNITUDE),
-    opA (NULL), opB (NULL),
+    opA (NULL), 
+    opB (NULL),
     initial_vector (NULL),
     transformation (NULL)
   {}
@@ -113,6 +114,9 @@ namespace SLEPcWrappers
       ierr = EPSSetOperators (solver_data->eps, *opA, PETSC_NULL);
     AssertThrow (ierr == 0, ExcSLEPcError(ierr));
 
+    // set runtime options
+    set_solver_type (solver_data->eps);
+
     // set the initial vector(s) if any
     if (initial_vector && initial_vector->size() != 0)
       {
@@ -129,10 +133,7 @@ namespace SLEPcWrappers
 
     // set transformation type if any
     if (transformation)
-      transformation->set_context(solver_data->eps);
-
-    // set runtime options
-    set_solver_type (solver_data->eps);
+      transformation->set_context (solver_data->eps);
 
     // set which portion of the eigenspectrum to solve for
     ierr = EPSSetWhichEigenpairs (solver_data->eps, set_which);
@@ -147,6 +148,15 @@ namespace SLEPcWrappers
     ierr = EPSSetFromOptions (solver_data->eps);
     AssertThrow (ierr == 0, ExcSLEPcError(ierr));
 
+    // Set convergence test to be absolute
+    ierr = EPSSetConvergenceTest (solver_data->eps, EPS_CONV_ABS);  
+    AssertThrow (ierr == 0, ExcSLEPcError(ierr));    
+
+    // Set the convergence test function
+    // ierr = EPSSetConvergenceTestFunction (solver_data->eps, &convergence_test, 
+    //                                           reinterpret_cast<void *>(&solver_control)); 
+    // AssertThrow (ierr == 0, ExcSLEPcError(ierr));   
+    
     // solve the eigensystem
     ierr = EPSSolve (solver_data->eps);
     AssertThrow (ierr == 0, ExcSLEPcError(ierr));
@@ -154,11 +164,39 @@ namespace SLEPcWrappers
     // get number of converged eigenstates
     ierr = EPSGetConverged (solver_data->eps,
 #ifdef PETSC_USE_64BIT_INDICES
-                            reinterpret_cast<PetscInt *>(n_converged));
+                            reinterpret_cast<PetscInt *>(n_converged)
 #else
-                            reinterpret_cast<int *>(n_converged));
+                           reinterpret_cast<int *>(n_converged)
 #endif
+                           );
     AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+    int n_iterations     = 0;
+    double residual_norm = 1e300;
+
+    // @todo Investigate elaborating on some of this to act on the
+    // complete eigenspectrum
+    {
+      // get the number of solver iterations
+      ierr = EPSGetIterationNumber (solver_data->eps, &n_iterations);  
+      AssertThrow (ierr == 0, ExcSLEPcError(ierr));  
+      
+      // get the residual norm of the most extreme eigenvalue
+      ierr = EPSComputeResidualNorm (solver_data->eps, 0, &residual_norm);
+      AssertThrow (ierr == 0, ExcSLEPcError(ierr));  
+
+      // check the solver state
+      const SolverControl::State state
+               = solver_control.check (n_iterations, residual_norm);
+
+      // get the solver state according to SLEPc
+      get_solver_state (state);
+      
+      // and in case of failure: throw exception
+      if (solver_control.last_check () != SolverControl::success) 
+       throw SolverControl::NoConvergence (solver_control.last_step (), 
+                                                   solver_control.last_value ()); 
+    }
   }
 
   void
@@ -178,7 +216,6 @@ namespace SLEPcWrappers
     AssertThrow (ierr == 0, ExcSLEPcError(ierr));
   }
 
-
   void
   SolverBase::reset ()
   {
@@ -189,14 +226,39 @@ namespace SLEPcWrappers
   }
 
   EPS *
-  SolverBase::get_EPS ()
+  SolverBase::get_eps ()
   {
-    if (solver_data.get() == 0)
+    if (solver_data.get () == 0)
       return NULL;
-
+    
     return &solver_data->eps;
   }
 
+  void
+  SolverBase::get_solver_state (const SolverControl::State state)
+  {
+    switch (state)
+      {
+      case ::dealii::SolverControl::iterate:
+       solver_data->reason = EPS_CONVERGED_ITERATING;
+       break;
+       
+      case ::dealii::SolverControl::success:
+       solver_data->reason = static_cast<EPSConvergedReason>(1);
+       break;
+       
+      case ::dealii::SolverControl::failure:
+       if (solver_control.last_step() > solver_control.max_steps())
+         solver_data->reason = EPS_DIVERGED_ITS;
+       else
+         solver_data->reason = EPS_DIVERGED_BREAKDOWN;
+       break;
+       
+      default:
+       Assert (false, ExcNotImplemented());
+      }
+  }
+
   /* ---------------------- SolverControls ----------------------- */
   SolverControl &
   SolverBase::control () const
@@ -205,42 +267,14 @@ namespace SLEPcWrappers
   }
 
   int
-  SolverBase::convergence_test (EPS                 /*eps*/,
-#ifdef PETSC_USE_64BIT_INDICES
-                                const PetscInt      iteration,
-#else
-                                const int           iteration,
-#endif
-                                const PetscReal     residual_norm,
-                                EPSConvergedReason *reason,
-                                void               *solver_control_x)
+  SolverBase::convergence_test (EPS          /*eps             */,
+                               PetscScalar  /*kr              */,
+                               PetscScalar  /*ki              */,
+                               PetscReal    /*residual_norm   */,
+                               PetscReal   */*estimated_error */,
+                                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;
-        else
-          *reason = EPS_DIVERGED_BREAKDOWN;
-        break;
-
-      default:
-        Assert (false, ExcNotImplemented());
-      }
+    // This function is undefined (future reference only).
 
     // return without failure.
     return 0;

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.