From c2a7e0d7d86094976d48cbb20574cc7937eb705c Mon Sep 17 00:00:00 2001
From: Axel Modave <axel.modave@gmail.com>
Date: Tue, 24 Jun 2014 10:52:16 +0000
Subject: [PATCH] pp
---
Solver/eigenSolver.cpp | 89 ++++++++++++++++++++----------------------
Solver/eigenSolver.h | 55 ++++++++++++++------------
2 files changed, 72 insertions(+), 72 deletions(-)
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index 90f8089a3a..17ee43190e 100644
--- a/Solver/eigenSolver.cpp
+++ b/Solver/eigenSolver.cpp
@@ -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
diff --git a/Solver/eigenSolver.h b/Solver/eigenSolver.h
index 98244a69f2..0a7a94e549 100644
--- a/Solver/eigenSolver.h
+++ b/Solver/eigenSolver.h
@@ -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
--
GitLab