Skip to content
Snippets Groups Projects
Commit c2a7e0d7 authored by Axel Modave's avatar Axel Modave
Browse files

pp

parent 6a2bcd55
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Solver/eigenSolver.cpp
+ 42
− 47
View file @ c2a7e0d7
......@@ -18,21 +18,21 @@ void eigenSolver::_try(int ierr) const
CHKERRABORT(PETSC_COMM_WORLD, ierr);
}
eigenSolver::eigenSolver(dofManager<double> *manager, std::string A,
std::string B, bool hermitian)
eigenSolver::eigenSolver(dofManager<double> *manager, std::string A, std::string B, bool hermitian)
: _A(0), _B(0), _hermitian(hermitian)
{
if(A.size()){
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.size()){
if (B.size()) {
_B = dynamic_cast<linearSystemPETSc<double>*>(manager->getLinearSystem(B));
if(!_B) Msg::Error("Could not find PETSc system '%s'", B.c_str());
}
}
eigenSolver::eigenSolver(linearSystemPETSc<double> *A,linearSystemPETSc<double> *B, bool hermitian) : _A(A), _B(B), _hermitian(hermitian){}
eigenSolver::eigenSolver(linearSystemPETSc<double> *A, linearSystemPETSc<double> *B,
bool hermitian) : _A(A), _B(B), _hermitian(hermitian) {}
bool eigenSolver::solve(int numEigenValues, std::string which, std::string method, double tolVal, int iterMax)
{
......@@ -79,24 +79,24 @@ bool eigenSolver::solve(int numEigenValues, std::string which, std::string metho
_try(EPSSetFromOptions(eps));
// force options specified directly as arguments
if(numEigenValues)
if (numEigenValues)
_try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
if(which == "smallest")
if (which=="smallest")
_try(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE));
else if(which == "smallestReal")
else if (which=="smallestReal")
_try(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL));
else if(which == "largest")
else if (which=="largest")
_try(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE));
// print info
#if (SLEPC_VERSION_RELEASE == 0 || (SLEPC_VERSION_MAJOR > 3 || (SLEPC_VERSION_MAJOR == 3 && SLEPC_VERSION_MINOR >= 4)))
#if (SLEPC_VERSION_RELEASE == 0 || (SLEPC_VERSION_MAJOR > 3 || (SLEPC_VERSION_MAJOR == 3 && SLEPC_VERSION_MINOR >= 4)))
EPSType type;
#else
#else
const EPSType type;
#endif
#endif
_try(EPSGetType(eps, &type));
Msg::Debug("SLEPc solution method: %s", type);
PetscInt nev;
_try(EPSGetDimensions(eps, &nev, PETSC_NULL, PETSC_NULL));
Msg::Debug("SLEPc number of requested eigenvalues: %d", nev);
......@@ -104,45 +104,45 @@ bool eigenSolver::solve(int numEigenValues, std::string which, std::string metho
PetscInt maxit;
_try(EPSGetTolerances(eps, &tol, &maxit));
Msg::Debug("SLEPc stopping condition: tol=%g, maxit=%d", tol, maxit);
// solve
Msg::Info("SLEPc solving...");
double t1 = Cpu();
_try(EPSSolve(eps));
// check convergence
int its;
_try(EPSGetIterationNumber(eps, &its));
EPSConvergedReason reason;
_try(EPSGetConvergedReason(eps, &reason));
if(reason == EPS_CONVERGED_TOL){
double t2 = Cpu();
if (reason==EPS_CONVERGED_TOL) {
double t2=Cpu();
Msg::Debug("SLEPc converged in %d iterations (%g s)", its, t2-t1);
}
else if(reason == EPS_DIVERGED_ITS)
else if (reason==EPS_DIVERGED_ITS)
Msg::Error("SLEPc diverged after %d iterations", its);
else if(reason == EPS_DIVERGED_BREAKDOWN)
else if (reason==EPS_DIVERGED_BREAKDOWN)
Msg::Error("SLEPc generic breakdown in method");
#if (SLEPC_VERSION_MAJOR < 3 || (SLEPC_VERSION_MAJOR == 3 && SLEPC_VERSION_MINOR < 2))
else if(reason == EPS_DIVERGED_NONSYMMETRIC)
else if (reason==EPS_DIVERGED_NONSYMMETRIC)
Msg::Error("The operator is nonsymmetric");
#endif
// get number of converged approximate eigenpairs
PetscInt nconv;
_try(EPSGetConverged(eps, &nconv));
Msg::Debug("SLEPc number of converged eigenpairs: %d", nconv);
// ignore additional eigenvalues if we get more than what we asked
if(nconv > nev) nconv = nev;
if (nconv>nev) nconv = nev;
if (nconv > 0) {
if (nconv>0) {
Vec xr, xi;
_try(MatGetVecs(A, PETSC_NULL, &xr));
_try(MatGetVecs(A, PETSC_NULL, &xi));
Msg::Debug(" Re[EigenValue] Im[EigenValue]"
" Relative error");
for (int i = 0; i < nconv; i++){
for (int i=0; i<nconv; i++) {
PetscScalar kr, ki;
_try(EPSGetEigenpair(eps, i, &kr, &ki, xr, xi));
PetscReal error;
......@@ -155,7 +155,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which, std::string metho
PetscReal im = ki;
#endif
Msg::Debug("EIG %03d %s%.16e %s%.16e %3.6e",
i, (re < 0) ? "" : " ", re, (im < 0) ? "" : " ", im, error);
i, (re<0) ? "" : " ", re, (im<0) ? "" : " ", im, error);
// store eigenvalues and eigenvectors
_eigenValues.push_back(std::complex<double>(re, im));
......@@ -163,14 +163,14 @@ bool eigenSolver::solve(int numEigenValues, std::string which, std::string metho
_try(VecGetArray(xr, &tmpr));
_try(VecGetArray(xi, &tmpi));
std::vector<std::complex<double> > ev(N);
for(int i = 0; i < N; i++){
for (int i=0; i<N; i++) {
#if defined(PETSC_USE_COMPLEX)
ev[i] = tmpr[i];
#else
ev[i] = std::complex<double>(tmpr[i], tmpi[i]);
#endif
}
_eigenVectors.push_back(ev);
_eigenVectors.push_back(ev);
}
_try(VecDestroy(&xr));
_try(VecDestroy(&xi));
......@@ -178,40 +178,35 @@ bool eigenSolver::solve(int numEigenValues, std::string which, std::string metho
_try(EPSDestroy(&eps));
if(reason == EPS_CONVERGED_TOL){
if (reason==EPS_CONVERGED_TOL) {
Msg::Debug("SLEPc done");
return true;
}
else{
else {
Msg::Warning("SLEPc failed");
return false;
}
}
void eigenSolver::normalize_mode(double scale){
void eigenSolver::normalize_mode(double scale) {
Msg::Info("Normalize all eigenvectors");
for (unsigned int i=0; i<_eigenVectors.size(); i++){
for (unsigned int i=0; i<_eigenVectors.size(); i++) {
double norm = 0.;
for (unsigned int j=0; j<_eigenVectors[i].size(); j++){
for (unsigned int j=0; j<_eigenVectors[i].size(); j++) {
std::complex<double> val = _eigenVectors[i][j];
double normval = std::abs(val);
if (normval >norm)
if (normval>norm)
norm = normval;
};
if (norm == 0) {
}
if (norm==0) {
Msg::Error("zero eigenvector");
return;
};
for (unsigned int j=0; j<_eigenVectors[i].size(); j++){
}
for (unsigned int j=0; j<_eigenVectors[i].size(); j++) {
_eigenVectors[i][j] *= (scale/norm);
};
};
};
}
}
}
#endif
Solver/eigenSolver.h
+ 30
− 25
View file @ c2a7e0d7
......@@ -16,7 +16,7 @@
#include "linearSystemPETSc.h"
class eigenSolver{
class eigenSolver {
private:
linearSystemPETSc<double> *_A, *_B;
bool _hermitian;
......@@ -30,46 +30,51 @@ class eigenSolver{
bool hermitian=true);
bool solve(int numEigenValues=0, std::string which="", std::string method="krylovschur",
double tolVal=1.e-7, int iterMax=20);
int getNumEigenValues(){ return _eigenValues.size(); }
std::complex<double> getEigenValue(int num){ return _eigenValues[num]; }
std::vector<std::complex<double> > &getEigenVector(int num){ return _eigenVectors[num]; }
void clear()
{
_eigenValues.clear();
_eigenVectors.clear();
};
std::complex<double> getEigenVectorComp(int num, int com)
{
int getNumEigenValues() {return _eigenValues.size();}
int getNumberEigenvectors() {return _eigenVectors.size();}
std::complex<double> getEigenValue(int num) {
return _eigenValues[num];
}
std::complex<double> getEigenVectorComp(int num, int com) {
return _eigenVectors[num][com];
};
int getNumberEigenvectors() {return _eigenVectors.size();};
std::vector<std::complex<double> > &getEigenVector(int num) {
return _eigenVectors[num];
}
void normalize_mode(double scale=1.);
void clear() {
_eigenValues.clear();
_eigenVectors.clear();
};
};
#else
#include "linearSystemPETSc.h"
class eigenSolver{
class eigenSolver {
private:
std::vector<std::complex<double> > _dummy;
public:
eigenSolver(dofManager<double> *manager, std::string A,
std::string B="", bool hermitian=false){}
eigenSolver(linearSystemPETSc<double> *A,linearSystemPETSc<double>* B = NULL,
bool hermitian=false){}
bool solve(int numEigenValues=0, std::string which="")
{
std::string B="", bool hermitian=false) {}
eigenSolver(linearSystemPETSc<double> *A, linearSystemPETSc<double>* B=NULL,
bool hermitian=false) {}
bool solve(int=0, std::string="", std::string="", double=0, int=0) {
Msg::Error("Eigen solver requires SLEPc");
return false;
}
int getNumEigenValues(){ return 0; }
std::complex<double> getEigenValue(int num){ return 0.; }
std::vector<std::complex<double> > &getEigenVector(int num){ return _dummy; }
void clear(){}
std::complex<double> getEigenVectorComp(int num, int com) { return 0.; }
int getNumberEigenvectors() {return 0;};
void normalize_mode(double scale=1.) {};
int getNumEigenValues() {return 0;}
int getNumberEigenvectors() {return 0;}
std::complex<double> getEigenValue(int num) {return 0.;}
std::complex<double> getEigenVectorComp(int num, int com) {return 0.;}
std::vector<std::complex<double> > &getEigenVector(int num) {return _dummy;}
void normalize_mode(double scale=1.) {}
void clear() {}
};
#endif
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment