Skip to content

Commit

Permalink
add function to sort eigensolver result
Browse files Browse the repository at this point in the history
  • Loading branch information
@qaumann
qaumann committed Apr 16, 2020
1 parent bcc3f19 commit cf1cae7
Show file tree
Hide file tree
Showing 2 changed files with 108 additions and 4 deletions.
98 changes: 98 additions & 0 deletions applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,93 @@ namespace { // helpers namespace

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<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](size_t i1, size_t i2) {return rEigenvalues[i1] < rEigenvalues[i2];});
} else if( mOrder == "sm") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](size_t i1, 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](size_t i1, size_t i2) {return rEigenvalues[i1] > rEigenvalues[i2];});
} else if( mOrder == "lm") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](size_t i1, 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( 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<size_t> idx(rEigenvalues.size());
std::iota(idx.begin(), idx.end(), 0);

// const std::string t = mParam["sort_order"].GetString();
if( mOrder == "sr" ) {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](size_t i1, size_t i2) {return std::real(rEigenvalues[i1]) < std::real(rEigenvalues[i2]);});
} else if( mOrder == "sm") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](size_t i1, size_t i2) {return std::abs(rEigenvalues[i1]) < std::abs(rEigenvalues[i2]);});
} else if( mOrder == "si") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](size_t i1, size_t i2) {return std::imag(rEigenvalues[i1]) < std::imag(rEigenvalues[i2]);});
} else if( mOrder == "lr" ) {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](size_t i1, size_t i2) {return std::real(rEigenvalues[i1]) > std::real(rEigenvalues[i2]);});
} else if( mOrder == "lm") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](size_t i1, size_t i2) {return std::abs(rEigenvalues[i1]) > std::abs(rEigenvalues[i2]);});
} else if( mOrder == "li") {
std::stable_sort(idx.begin(), idx.end(),
[&rEigenvalues](size_t i1, 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( 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 @@ -129,6 +216,8 @@ class FEASTEigensystemSolver
"number_of_eigenvalues" : 0,
"search_lowest_eigenvalues" : false,
"search_highest_eigenvalues" : false,
"sort_eigenvalues" : false,
"sort_order" : "sr",
"subspace_size" : 0,
"max_iteration" : 20,
"tolerance" : 1e-12,
Expand Down Expand Up @@ -158,6 +247,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 Down Expand Up @@ -260,6 +353,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

0 comments on commit cf1cae7

Please sign in to comment.