Skip to content

Commit

Permalink
Merge pull request #6794 from KratosMultiphysics/eigen-solvers/update…
Browse files Browse the repository at this point in the history 
…-feast-settings

[EigenSolversApplication] Extend FEAST interface
  • Loading branch information
@qaumann
qaumann authored May 4, 2020
2 parents 5beae57 + 30b50c0 commit 5f7736b
Show file tree
Hide file tree
Showing 3 changed files with 122 additions and 34 deletions.
130 changes: 108 additions & 22 deletions applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,15 +46,18 @@ namespace { // helpers namespace
Parameters SettingsHelper<double>::GetDefaultParameters()
{
return Parameters(R"({
"e_min" : 0.0,
"e_max" : 0.0
})");
}

template<>
Parameters SettingsHelper<std::complex<double>>::GetDefaultParameters()
{
return Parameters(R"({
"e_mid_re" : 0.0,
"e_mid_im" : 0.0,
"e_r" : 0.0
})");
}

Expand All @@ -66,10 +69,6 @@ namespace { // helpers namespace

KRATOS_ERROR_IF( mParam["e_max"].GetDouble() <= mParam["e_min"].GetDouble() ) <<
"Invalid eigenvalue limits provided" << std::endl;

KRATOS_INFO_IF( "FEASTEigensystemSolver",
mParam["e_mid_re"].GetDouble() != 0.0 || mParam["e_mid_im"].GetDouble() != 0.0 || mParam["e_r"].GetDouble() != 0.0 ) <<
"Manually defined e_mid_re, e_mid_im, e_r are not used for real symmetric matrices" << std::endl;
}

template<>
Expand All @@ -78,20 +77,101 @@ namespace { // helpers namespace
KRATOS_ERROR_IF( mParam["e_r"].GetDouble() <= 0.0 ) <<
"Invalid search radius provided" << std::endl;

KRATOS_INFO_IF( "FEASTEigensystemSolver",
mParam["e_min"].GetDouble() != 0.0 || mParam["e_max"].GetDouble() != 0.0 ) <<
"Manually defined e_min, e_max are not used for complex symmetric matrices" << std::endl;

KRATOS_INFO_IF( "FEASTEigensystemSolver",
mParam["search_lowest_eigenvalues"].GetBool() || mParam["search_highest_eigenvalues"].GetBool() ) <<
"Search for extremal eigenvalues is only available for Hermitian problems" << std::endl;
KRATOS_ERROR_IF( mParam["search_lowest_eigenvalues"].GetBool() || mParam["search_highest_eigenvalues"].GetBool() ) <<
"Search for extremal eigenvalues is only available for real symmetric problems" << std::endl;
}

template<> double SettingsHelper<double>::GetE1() {return mParam["e_min"].GetDouble();}
template<> std::complex<double> SettingsHelper<std::complex<double>>::GetE1() {return std::complex<double>(mParam["e_mid_re"].GetDouble(), mParam["e_mid_im"].GetDouble());}

template<>double SettingsHelper<double>::GetE2() {return mParam["e_max"].GetDouble();}
template<> double SettingsHelper<std::complex<double>>::GetE2() {return mParam["e_r"].GetDouble();}

template<typename TScalar>
struct SortingHelper
{
SortingHelper(std::string Order) : mOrder(Order) {};
// void Check();
template<typename MatrixType, typename VectorType>
void SortEigenvalues(VectorType&, MatrixType&);
private:
std::string mOrder;
};

template<> template<typename MatrixType, typename VectorType>
void SortingHelper<double>::SortEigenvalues(VectorType &rEigenvalues, MatrixType &rEigenvectors)
{
KRATOS_WARNING_IF("FeastEigensystemSolver", mOrder == "si") << "Attempting to sort by imaginary value. Falling back on \"sr\"" << std::endl;
KRATOS_WARNING_IF("FeastEigensystemSolver", mOrder == "li") << "Attempting to sort by imaginary value. Falling back on \"lr\"" << std::endl;

std::vector<std::size_t> idx(rEigenvalues.size());
std::iota(idx.begin(), idx.end(), 0);

if( mOrder == "sr" || mOrder == "si" ) {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](std::size_t i1, std::size_t i2) {return rEigenvalues[i1] < rEigenvalues[i2];});
} else if( mOrder == "sm") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](std::size_t i1, std::size_t i2) {return std::abs(rEigenvalues[i1]) < std::abs(rEigenvalues[i2]);});
} else if( mOrder == "lr" || mOrder == "li" ) {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](std::size_t i1, std::size_t i2) {return rEigenvalues[i1] > rEigenvalues[i2];});
} else if( mOrder == "lm") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](std::size_t i1, std::size_t i2) {return std::abs(rEigenvalues[i1]) > std::abs(rEigenvalues[i2]);});
} else {
KRATOS_ERROR << "Invalid sort type. Allowed are sr, sm, si, lr, lm, li" << std::endl;
}

VectorType tmp_eigenvalues(rEigenvalues.size());
MatrixType tmp_eigenvectors(rEigenvectors.size1(), rEigenvectors.size2());

for( std::size_t i=0; i<rEigenvalues.size(); ++i ) {
tmp_eigenvalues[i] = rEigenvalues[idx[i]];
column(tmp_eigenvectors, i).swap(column(rEigenvectors, idx[i]));
}
rEigenvalues.swap(tmp_eigenvalues);
rEigenvectors.swap(tmp_eigenvectors);
}

template<> template<typename MatrixType, typename VectorType>
void SortingHelper<std::complex<double>>::SortEigenvalues(VectorType &rEigenvalues, MatrixType &rEigenvectors)
{
std::vector<std::size_t> idx(rEigenvalues.size());
std::iota(idx.begin(), idx.end(), 0);

if( mOrder == "sr" ) {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](std::size_t i1, std::size_t i2) {return std::real(rEigenvalues[i1]) < std::real(rEigenvalues[i2]);});
} else if( mOrder == "sm") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](std::size_t i1, std::size_t i2) {return std::abs(rEigenvalues[i1]) < std::abs(rEigenvalues[i2]);});
} else if( mOrder == "si") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](std::size_t i1, std::size_t i2) {return std::imag(rEigenvalues[i1]) < std::imag(rEigenvalues[i2]);});
} else if( mOrder == "lr" ) {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](std::size_t i1, std::size_t i2) {return std::real(rEigenvalues[i1]) > std::real(rEigenvalues[i2]);});
} else if( mOrder == "lm") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](std::size_t i1, std::size_t i2) {return std::abs(rEigenvalues[i1]) > std::abs(rEigenvalues[i2]);});
} else if( mOrder == "li") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](std::size_t i1, std::size_t i2) {return std::imag(rEigenvalues[i1]) > std::imag(rEigenvalues[i2]);});
} else {
KRATOS_ERROR << "Invalid sort type. Allowed are sr, sm, si, lr, lm, li" << std::endl;
}

VectorType tmp_eigenvalues(rEigenvalues.size());
MatrixType tmp_eigenvectors(rEigenvectors.size1(), rEigenvectors.size2());

for( std::size_t i=0; i<rEigenvalues.size(); ++i ) {
tmp_eigenvalues[i] = rEigenvalues[idx[i]];
column(tmp_eigenvectors, i).swap(column(rEigenvectors, idx[i]));
}
rEigenvalues.swap(tmp_eigenvalues);
rEigenvectors.swap(tmp_eigenvectors);
}
}

template<
Expand Down Expand Up @@ -135,11 +215,8 @@ class FEASTEigensystemSolver
"number_of_eigenvalues" : 0,
"search_lowest_eigenvalues" : false,
"search_highest_eigenvalues" : false,
"e_min" : 0.0,
"e_max" : 0.0,
"e_mid_re" : 0.0,
"e_mid_im" : 0.0,
"e_r" : 0.0,
"sort_eigenvalues" : false,
"sort_order" : "sr",
"subspace_size" : 0,
"max_iteration" : 20,
"tolerance" : 1e-12,
Expand Down Expand Up @@ -169,6 +246,10 @@ class FEASTEigensystemSolver
mParam["number_of_eigenvalues"].GetInt() > 0 && mParam["subspace_size"].GetInt() > 0 ) <<
"Manually defined subspace size will be overwritten to match the defined number of eigenvalues" << std::endl;

const std::string s = mParam["sort_order"].GetString();
KRATOS_ERROR_IF( !(s=="sr" || s=="sm" || s=="si" || s=="lr" || s=="lm" || s=="li") ) <<
"Invalid sort type. Allowed are sr, sm, si, lr, lm, li" << std::endl;

SettingsHelper<TScalarOut>(mParam).CheckParameters();
}

Expand All @@ -189,14 +270,14 @@ class FEASTEigensystemSolver
{
// settings
const std::size_t system_size = rK.size1();
size_t subspace_size;
std::size_t subspace_size;

if( mParam["search_lowest_eigenvalues"].GetBool() || mParam["search_highest_eigenvalues"].GetBool() ) {
subspace_size = 2 * static_cast<size_t>(mParam["number_of_eigenvalues"].GetInt());
subspace_size = 2 * static_cast<std::size_t>(mParam["number_of_eigenvalues"].GetInt());
} else if( mParam["subspace_size"].GetInt() == 0 ) {
subspace_size = 1.5 * static_cast<size_t>(mParam["number_of_eigenvalues"].GetInt());
subspace_size = 1.5 * static_cast<std::size_t>(mParam["number_of_eigenvalues"].GetInt());
} else {
subspace_size = static_cast<size_t>(mParam["subspace_size"].GetInt());
subspace_size = static_cast<std::size_t>(mParam["subspace_size"].GetInt());
}

// create column based matrix for the fortran routine
Expand Down Expand Up @@ -271,6 +352,11 @@ class FEASTEigensystemSolver
tmp_eigenvalues.resize(M, true);
tmp_eigenvectors.resize(system_size, M, true);

// sort if required
if( mParam["sort_eigenvalues"].GetBool() ) {
SortingHelper<TScalarOut>(mParam["sort_order"].GetString()).SortEigenvalues(tmp_eigenvalues, tmp_eigenvectors);
}

// copy eigenvalues to result vector
rEigenvalues.swap(tmp_eigenvalues);

Expand Down
14 changes: 10 additions & 4 deletions applications/EigenSolversApplication/tests/test_feast_eigensystem_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,8 @@ def test_real_general_gev_complex_result(self):
"e_mid_re": 0.0,
"e_mid_im": 0.0,
"e_r": 30.0,
"sort_eigenvalues": true,
"sort_order": "li",
"echo_level": 0
}''')

Expand All @@ -168,8 +170,8 @@ def test_real_general_gev_complex_result(self):
# Solve
eigen_solver.Solve(K, M, eigenvalues, eigenvectors)

self.assertAlmostEqual(eigenvalues[0], 0.0-1.0j, 7)
self.assertAlmostEqual(eigenvalues[1], 0.0+1.0j, 7)
self.assertAlmostEqual(eigenvalues[0], 0.0+1.0j, 7)
self.assertAlmostEqual(eigenvalues[1], 0.0-1.0j, 7)

Kc = KratosMultiphysics.ComplexCompressedMatrix(K)
Mc = KratosMultiphysics.ComplexCompressedMatrix(M)
Expand Down Expand Up @@ -199,6 +201,8 @@ def test_complex_symmetric_gev(self):
"e_mid_re": 0.0,
"e_mid_im": 0.0,
"e_r": 0.16,
"sort_eigenvalues": true,
"sort_order": "lm",
"echo_level": 0
}''')

Expand Down Expand Up @@ -258,6 +262,8 @@ def test_complex_general_gev(self):
"e_mid_re": 10.0,
"e_mid_im": 0.0,
"e_r": 3.0,
"sort_eigenvalues": true,
"sort_order": "sm",
"echo_level": 0
}''')

Expand Down Expand Up @@ -295,8 +301,8 @@ def test_complex_general_gev(self):
self.assertEqual(eigenvectors.Size1(), 2)
self.assertEqual(eigenvectors.Size2(), 5)

self.assertAlmostEqual(eigenvalues[0], 12.0+1.2j, 7)
self.assertAlmostEqual(eigenvalues[1], 10.5+1.05j, 7)
self.assertAlmostEqual(eigenvalues[0], 10.5+1.05j, 7)
self.assertAlmostEqual(eigenvalues[1], 12.0+1.2j, 7)

for i in range(eigenvalues.Size()):
eigenvector = KratosMultiphysics.ComplexVector(n)
Expand Down
12 changes: 4 additions & 8 deletions kratos/python_scripts/eigen_solver_factory.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,8 @@ def ConstructSolver(settings):
if kratos_utils.CheckIfApplicationsAvailable("EigenSolversApplication"):
import KratosMultiphysics.EigenSolversApplication as EiSA
if EiSA.HasFEAST():
if settings.Has("symmetric") and settings["symmetric"].GetBool():
eigen_solver = EiSA.FEASTSymmetricEigensystemSolver(settings)
else:
eigen_solver = EiSA.FEASTGeneralEigensystemSolver(settings)
is_symmetric = settings["symmetric"].GetBool() if settings.Has("symmetric") else True
eigen_solver = EiSA.FEASTSymmetricEigensystemSolver(settings) if is_symmetric else EiSA.FEASTGeneralEigensystemSolver(settings)
return eigen_solver
else:
raise Exception("FEAST not available in EigenSolversApplication")
Expand All @@ -36,10 +34,8 @@ def ConstructSolver(settings):
if kratos_utils.CheckIfApplicationsAvailable("EigenSolversApplication"):
import KratosMultiphysics.EigenSolversApplication as EiSA
if EiSA.HasFEAST():
if settings.Has("symmetric") and settings["symmetric"].GetBool():
eigen_solver = EiSA.ComplexFEASTSymmetricEigensystemSolver(settings)
else:
eigen_solver = EiSA.ComplexFEASTGeneralEigensystemSolver(settings)
is_symmetric = settings["symmetric"].GetBool() if settings.Has("symmetric") else True
eigen_solver = EiSA.ComplexFEASTSymmetricEigensystemSolver(settings) if is_symmetric else EiSA.ComplexFEASTGeneralEigensystemSolver(settings)
return eigen_solver
else:
raise Exception("FEAST not available in EigenSolversApplication")
Expand Down

0 comments on commit 5f7736b

Please sign in to comment.