initial version of eigensolver using slepc

Common/GmshConfig.h.in

@@ -38,7 +38,7 @@
 #cmakedefine HAVE_OSMESA
 #cmakedefine HAVE_PETSC
 #cmakedefine HAVE_QT
-#cmakedefine HAVE_SLEPSC
+#cmakedefine HAVE_SLEPC
 #cmakedefine HAVE_SOLVER
 #cmakedefine HAVE_TAUCS
 #cmakedefine HAVE_TETGEN

Solver/CMakeLists.txt

@@ -9,6 +9,7 @@ set(SRC
   elasticityTerm.cpp
   elasticitySolver.cpp
   SElement.cpp
+  eigenSolve.cpp
 )
 
 file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h) 

Solver/dofManager.h

@@ -73,7 +73,7 @@ class dofManager{
   std::map<Dof, Dof> associatedWith;
     
   // linearSystems
-  std::map<const std::string, linearSystem<dataMat>* > _linearSystems;
+  std::map<const std::string, linearSystem<dataMat>*> _linearSystems;
   linearSystem<dataMat>* _current;
    
  public:
@@ -236,7 +236,12 @@ class dofManager{
   }  
   linearSystem<dataMat> *getLinearSystem(std::string &name)
   {
-    return _linearSystems[name];
+    typename std::map<const std::string, linearSystem<dataMat>*>::iterator it =
+      _linearSystems.find(name);
+    if(it != _linearSystems.end())
+      return it->second;
+    else
+      return 0;
   }
 };
 

Solver/eigenSolve.cpp

@@ -9,15 +9,16 @@
 
 #include <slepceps.h>
 
-eigenSolve::eigenSolve(dofManager *manager, const char *A, const char *B)
+eigenSolve::eigenSolve(dofManager<double, double> *manager, std::string A, 
+                       std::string B) : _A(0), _B(0)
 {
-  if(A){
-    _A = manager->getLinearSystem(A);
-    if(!_A) Msg::Error("Could not find system A");
+  if(A.size()){
+    _A = dynamic_cast<linearSystemPETSc<double>*>(manager->getLinearSystem(A));
+    if(!_A) Msg::Error("Could not find PETSc system '%s'", A.c_str());
   }
-  if(B){
-    _B = manager->getLinearSystem(B);
-    if(!_B) Msg::Error("Could not find system B");
+  if(B.size()){
+    _B = dynamic_cast<linearSystemPETSc<double>*>(manager->getLinearSystem(B));
+    if(!_B) Msg::Error("Could not find PETSc system '%s'", B.c_str());
   }
 }
 
@@ -25,41 +26,63 @@ void eigenSolve::solve()
 {
   if(!_A) return;
   Mat A = _A->getMatrix();
-  Mat B = _B ? _b->getMatrix() : PETSC_NULL;
+  Mat B = _B ? _B->getMatrix() : PETSC_NULL;
+
+  _try(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+  _try(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+  PetscInt N, M;
+  _try(MatGetSize(A, &N, &M));
 
   // treatment of a generalized eigenvalue problem A x - \lambda B x = 0
   EPS eps;
   _try(EPSCreate(PETSC_COMM_WORLD, &eps));
   _try(EPSSetOperators(eps, A, B));
-  _try(EPSSetProblemType(eps, EPS_GNHEP));
+  // FIXME: check if hermitian (EPS_HEP or EPS_GHEP)
+  _try(EPSSetProblemType(eps, _B ? EPS_GNHEP : EPS_NHEP));
+
+  // set options
+  _try(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE));
+  //_try(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE));
+  //_try(EPSSetTarget(eps, 2.)); // find eigenvalues close to given target
+  PetscInt nev = 2; // number of eigenvalues to compute
+  PetscInt mpd = nev; // max dim allowed for the projected problem
+  PetscInt ncv = nev + mpd; // max dim of the subspace to be used by the solver
+  _try(EPSSetDimensions(eps, nev, ncv, mpd));
 
-  // set solver parameters at runtime
+  // override options at runtime, petsc-style
   _try(EPSSetFromOptions(eps));
+
+  // print info
+  const EPSType type;
+  _try(EPSGetType(eps, &type));
+  Msg::Info("SLEPc solution method: %s", type);
+  _try(EPSGetDimensions(eps, &nev, &ncv, &mpd));
+  Msg::Info("SLEPc number of requested eigenvalues: %d", nev);
+  PetscReal tol;
+  PetscInt maxit;
+  _try(EPSGetTolerances(eps, &tol, &maxit));
+  Msg::Info("SLEPc stopping condition: tol=%g, maxit=%d", tol, maxit);
  	
   Msg::Info("SLEPc solving...");
   _try(EPSSolve(eps));
+
   int its;
   _try(EPSGetIterationNumber(eps, &its));
-  Msg::Info("Number of iterations of the method: %d", its);
+  EPSConvergedReason reason;
+  _try(EPSGetConvergedReason(eps, &reason));
+  if(reason == EPS_CONVERGED_TOL)
+    Msg::Info("SLEPc converged in %d iterations", its);
+  else if(reason == EPS_DIVERGED_ITS)
+    Msg::Error("SLEPc did not converge in %d iterations", its);
+  else if(reason == EPS_DIVERGED_BREAKDOWN)
+    Msg::Error("SLEPc generic breakdown in method");
+  else if(reason == EPS_DIVERGED_NONSYMMETRIC)
+    Msg::Error("The operator is nonsymmetric");
   
-  // get some information from the solver and display it
-  EPSType type;
-  _try(EPSGetType(eps, &type));
-  Msg::Info("Solution method: %s", type);
-  int nev;
-  _try(EPSGetDimensions(eps, &nev, PETSC_NULL));
-  Msg::Info("Number of requested eigenvalues: %d", nev);
-  PetscReal tol;
-  int maxit;
-  _try(EPSGetTolerances(eps, &tol, &maxit));
-  Msg::Info("Stopping condition: tol=%.4g, maxit=%d", tol, maxit);
-   
   // get number of converged approximate eigenpairs
-  int nconv;
+  PetscInt nconv;
   _try(EPSGetConverged(eps, &nconv));
-  Msg::Info("Number of converged eigenpairs: %d", nconv);
-
-  Msg::Info("Re[Eigenvalue] Im[Eigenvalue] Rel. error = ||Ax - eig*Bx||/||eig*x||");
+  Msg::Info("SLEPc number of converged eigenpairs: %d", nconv);
 
   // if number of converged solutions is greated than requested
   // discarded surplus solutions
@@ -71,7 +94,11 @@ void eigenSolve::solve()
   Vec xr, xi;
   _try(MatGetVecs(A, PETSC_NULL, &xr));
   _try(MatGetVecs(A, PETSC_NULL, &xi));
+  Msg::Info("         Re[EigenValue]          Im[EigenValue]"
+            "          Relative error");
   for (int i = nconv - 1; i > nconv - nev - 1; i--){
+    std::vector<std::complex<double> > ev(N);
+
     PetscScalar kr, ki;
     _try(EPSGetEigenpair(eps, i, &kr, &ki, xr, xi));
     PetscReal error;
@@ -83,14 +110,19 @@ void eigenSolve::solve()
     PetscReal re = kr;
     PetscReal im = ki;
 #endif
-    // Output eigenvalues
-    Msg::Info("EIG %03d  %s  %.16e  %s  %.16e  %3.6e", 
-              (nconv-i), (re<0)? "" : " ", re, (im<0)? "" : " ", im,  error);
+    // Output eigenvalues and relative error (= ||Ax - eig*Bx||/||eig*x||)
+    Msg::Info("EIG %03d %s%.16e %s%.16e  %3.6e", (nconv - i), 
+              (re < 0) ? "" : " ", re, (im < 0) ? "" : " ", im, error);
+    _eigenValues.push_back(std::complex<double>(re, im));
 
     // Output eigenvector
-    PetscScalar *real_tmp, *imag_tmp;
-    _try(VecGetArray(xr, &real_tmp));
-    _try(VecGetArray(xi, &imag_tmp));
+    PetscScalar *tmpr, *tmpi;
+    _try(VecGetArray(xr, &tmpr));
+    _try(VecGetArray(xi, &tmpi));
+    for(int i = 0; i < N; i++)
+      ev[i] = std::complex<double>(PetscRealPart(tmpr[i]), 
+                                   PetscRealPart(tmpi[i]));
+    _eigenVectors.push_back(ev);
   }
   
   // Free work space

Solver/eigenSolve.h

@@ -6,6 +6,7 @@
 #ifndef _EIGEN_SOLVE_H_
 #define _EIGEN_SOLVE_H_
 
+#include <string>
 #include <complex>
 #include "GmshConfig.h"
 #include "GmshMessage.h"
@@ -13,19 +14,20 @@
 
 #if defined(HAVE_SLEPC)
 
+#include "linearSystemPETSc.h"
 #include <slepc.h>
 
 class eigenSolve{
  private:
-  linearSolver *_A, *_B;
+  linearSystemPETSc<double> *_A, *_B;
   std::vector<std::complex<double> > _eigenValues;
   std::vector<std::vector<std::complex<double> > > _eigenVectors;
   void _try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); }
  public:
-  eigenSolve(dofManager<double, double> *manager, const char *A, 
-             const char *B=0);
-  solve();
-  std::complex<double> getEigenValue(int num){ return _eigenValuesp[num]; }
+  eigenSolve(dofManager<double, double> *manager, std::string A, 
+             std::string B="");
+  void solve();
+  std::complex<double> getEigenValue(int num){ return _eigenValues[num]; }
   std::vector<std::complex<double> > &getEigenVector(int num){ return _eigenVectors[num]; }
 };
 
@@ -35,11 +37,11 @@ class eigenSolve{
  private:
   std::vector<std::complex<double> > _dummy;
  public:
-  eigenSolve(dofManager<double, double> *manager, const char *A, 
-             const char *B=0){}
-  void solve(){}{ Msg::Error("Eigen solver requires SLEPc"); }
+  eigenSolve(dofManager<double, double> *manager, std::string A, 
+             std::string B=""){}
+  void solve(){ Msg::Error("Eigen solver requires SLEPc"); }
   std::complex<double> getEigenValue(int num){ return 0.; }
-  std::vector<complex<double> > &getEigenVector(int num){ return _dummy; }
+  std::vector<std::complex<double> > &getEigenVector(int num){ return _dummy; }
 };
 
 #endif

Solver/linearSystemCSR.cpp

@@ -6,10 +6,10 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
-#include <time.h>
 #include "GmshConfig.h"
 #include "GmshMessage.h"
 #include "linearSystemCSR.h"
+#include "OS.h"
 
 #define SWAP(a, b)  temp = (a); (a) = (b); (b) = temp;
 #define SWAPI(a, b) tempi = (a); (a) = (b); (b) = tempi;
@@ -343,7 +343,7 @@ int linearSystemCSRTaucs<double>::systemSolve()
   myVeryCuteTaucsMatrix.flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
 
   char* options[] = { "taucs.factor.LLT=true", NULL };
-  clock_t t1 = clock();
+  double t1 = Cpu();
   int result = taucs_linsolve(&myVeryCuteTaucsMatrix,
                               NULL,
                               1,
@@ -351,9 +351,8 @@ int linearSystemCSRTaucs<double>::systemSolve()
                               &(*_b)[0],
                               options,
                               NULL);
-  clock_t t2 = clock();
-  Msg::Info("TAUCS has solved %d unknowns in %8.3f seconds",
-            _b->size(), (double)(t2 - t1) / CLOCKS_PER_SEC);
+  double t2 = Cpu();
+  Msg::Info("TAUCS has solved %d unknowns in %8.3f seconds", _b->size(), t2 - t1);
   if (result != TAUCS_SUCCESS){
     Msg::Error("Taucs Was Not Successfull %d", result);
   }