diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index 2fadf56577858393c8583ca1b60c799dc2319d95..3b484fee7b53c20e474eb4bc9e541c891cc8c735 100644
--- a/Common/GmshConfig.h.in
+++ b/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
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index 699c6edc02f67cf3445df34b25fbee8c770bcb08..48363d84575b286d01d5d05e7b3b64a1a2c45850 100644
--- a/Solver/CMakeLists.txt
+++ b/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)
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 525281004bd4cc699dbbf4ed6ba27da37b66fbc5..20f7aa5a5ef59b38fc5146135632cb35a8f864d4 100644
--- a/Solver/dofManager.h
+++ b/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;
}
};
diff --git a/Solver/eigenSolve.cpp b/Solver/eigenSolve.cpp
index 2a61b171621496799c48f903066394146a641a3b..d305f799b5a3d2c95b36604e157541223d9a520a 100644
--- a/Solver/eigenSolve.cpp
+++ b/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
diff --git a/Solver/eigenSolve.h b/Solver/eigenSolve.h
index 1e02ce35c7ba179e54e939ba28487a78a1498599..1ee5299f30fae60592654667f2efc9dad79fdb6d 100644
--- a/Solver/eigenSolve.h
+++ b/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
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 3f83b68c98f6acee7cb8b41290d0fb5646772573..4bceac89cab8bf2a9ca6605962e13a1763621c84 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/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);
}