diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp index 73fe3452d2885ff538535788aff16fed54c8b596..4e54e12c0bf7b236e36dbeff7137d125a06525b0 100644 --- a/Geo/GFaceCompound.cpp +++ b/Geo/GFaceCompound.cpp @@ -29,7 +29,7 @@ #include "CreateFile.h" #include "Context.h" #include "discreteFace.h" -#include "eigenSolve.h" +#include "eigenSolver.h" static void fixEdgeToValue(GEdge *ed, double value, dofManager<double, double> &myAssembler) { diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt index 48363d84575b286d01d5d05e7b3b64a1a2c45850..2fe862d6ba39e46c3275425cb360e4c1cf60646e 100644 --- a/Solver/CMakeLists.txt +++ b/Solver/CMakeLists.txt @@ -9,7 +9,7 @@ set(SRC elasticityTerm.cpp elasticitySolver.cpp SElement.cpp - eigenSolve.cpp + eigenSolver.cpp ) file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h) diff --git a/Solver/eigenSolve.cpp b/Solver/eigenSolver.cpp similarity index 85% rename from Solver/eigenSolve.cpp rename to Solver/eigenSolver.cpp index d305f799b5a3d2c95b36604e157541223d9a520a..64176c5764e8c040755be1404b30e70a223a9d33 100644 --- a/Solver/eigenSolve.cpp +++ b/Solver/eigenSolver.cpp @@ -3,14 +3,14 @@ // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. -#include "eigenSolve.h" +#include "eigenSolver.h" #if defined(HAVE_SLEPC) #include <slepceps.h> -eigenSolve::eigenSolve(dofManager<double, double> *manager, std::string A, - std::string B) : _A(0), _B(0) +eigenSolver::eigenSolver(dofManager<double, double> *manager, std::string A, + std::string B) : _A(0), _B(0) { if(A.size()){ _A = dynamic_cast<linearSystemPETSc<double>*>(manager->getLinearSystem(A)); @@ -22,7 +22,7 @@ eigenSolve::eigenSolve(dofManager<double, double> *manager, std::string A, } } -void eigenSolve::solve() +void eigenSolver::solve() { if(!_A) return; Mat A = _A->getMatrix(); @@ -38,16 +38,18 @@ void eigenSolve::solve() _try(EPSCreate(PETSC_COMM_WORLD, &eps)); _try(EPSSetOperators(eps, A, B)); // FIXME: check if hermitian (EPS_HEP or EPS_GHEP) - _try(EPSSetProblemType(eps, _B ? EPS_GNHEP : EPS_NHEP)); + //_try(EPSSetProblemType(eps, _B ? EPS_GNHEP : EPS_NHEP)); + _try(EPSSetProblemType(eps, _B ? EPS_GHEP : EPS_HEP)); // set options - _try(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE)); - //_try(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE)); - //_try(EPSSetTarget(eps, 2.)); // find eigenvalues close to given target + //_try(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE)); + //_try(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL)); + _try(EPSSetTarget(eps, 0.)); // 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)); + //_try(EPSSetType(eps, EPSPOWER)); // override options at runtime, petsc-style _try(EPSSetFromOptions(eps)); @@ -98,7 +100,6 @@ void eigenSolve::solve() " 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; @@ -119,9 +120,13 @@ void eigenSolve::solve() 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])); + 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); } diff --git a/Solver/eigenSolve.h b/Solver/eigenSolver.h similarity index 77% rename from Solver/eigenSolve.h rename to Solver/eigenSolver.h index 1ee5299f30fae60592654667f2efc9dad79fdb6d..a58deb308fb4f745893db64f133925b78e90efd5 100644 --- a/Solver/eigenSolve.h +++ b/Solver/eigenSolver.h @@ -17,29 +17,31 @@ #include "linearSystemPETSc.h" #include <slepc.h> -class eigenSolve{ +class eigenSolver{ private: 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, std::string A, - std::string B=""); + eigenSolver(dofManager<double, double> *manager, std::string A, + std::string B=""); void solve(); + 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]; } }; #else -class eigenSolve{ +class eigenSolver{ private: std::vector<std::complex<double> > _dummy; public: - eigenSolve(dofManager<double, double> *manager, std::string A, - std::string B=""){} + eigenSolver(dofManager<double, double> *manager, std::string A, + std::string B=""){} void solve(){ Msg::Error("Eigen solver requires SLEPc"); } + int getNumEigenValues(){ return 0.; } std::complex<double> getEigenValue(int num){ return 0.; } std::vector<std::complex<double> > &getEigenVector(int num){ return _dummy; } };