diff --git a/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h b/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h index 4df28e5afdd8..bdc5b4f0fd85 100644 --- a/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h +++ b/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h @@ -46,7 +46,8 @@ namespace { // helpers namespace Parameters SettingsHelper::GetDefaultParameters() { return Parameters(R"({ - + "e_min" : 0.0, + "e_max" : 0.0 })"); } @@ -54,7 +55,9 @@ namespace { // helpers namespace Parameters SettingsHelper>::GetDefaultParameters() { return Parameters(R"({ - + "e_mid_re" : 0.0, + "e_mid_im" : 0.0, + "e_r" : 0.0 })"); } @@ -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<> @@ -78,13 +77,8 @@ 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::GetE1() {return mParam["e_min"].GetDouble();} @@ -92,6 +86,92 @@ namespace { // helpers namespace template<>double SettingsHelper::GetE2() {return mParam["e_max"].GetDouble();} template<> double SettingsHelper>::GetE2() {return mParam["e_r"].GetDouble();} + + template + struct SortingHelper + { + SortingHelper(std::string Order) : mOrder(Order) {}; + // void Check(); + template + void SortEigenvalues(VectorType&, MatrixType&); + private: + std::string mOrder; + }; + + template<> template + void SortingHelper::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 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 template + void SortingHelper>::SortEigenvalues(VectorType &rEigenvalues, MatrixType &rEigenvectors) + { + std::vector 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 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(mParam).CheckParameters(); } @@ -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(mParam["number_of_eigenvalues"].GetInt()); + subspace_size = 2 * static_cast(mParam["number_of_eigenvalues"].GetInt()); } else if( mParam["subspace_size"].GetInt() == 0 ) { - subspace_size = 1.5 * static_cast(mParam["number_of_eigenvalues"].GetInt()); + subspace_size = 1.5 * static_cast(mParam["number_of_eigenvalues"].GetInt()); } else { - subspace_size = static_cast(mParam["subspace_size"].GetInt()); + subspace_size = static_cast(mParam["subspace_size"].GetInt()); } // create column based matrix for the fortran routine @@ -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(mParam["sort_order"].GetString()).SortEigenvalues(tmp_eigenvalues, tmp_eigenvectors); + } + // copy eigenvalues to result vector rEigenvalues.swap(tmp_eigenvalues); diff --git a/applications/EigenSolversApplication/tests/test_feast_eigensystem_solver.py b/applications/EigenSolversApplication/tests/test_feast_eigensystem_solver.py index 1cb5f857a8f4..4b5032f78b36 100644 --- a/applications/EigenSolversApplication/tests/test_feast_eigensystem_solver.py +++ b/applications/EigenSolversApplication/tests/test_feast_eigensystem_solver.py @@ -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 }''') @@ -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) @@ -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 }''') @@ -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 }''') @@ -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) diff --git a/kratos/python_scripts/eigen_solver_factory.py b/kratos/python_scripts/eigen_solver_factory.py index e6aef38864e5..d6d17df9e977 100644 --- a/kratos/python_scripts/eigen_solver_factory.py +++ b/kratos/python_scripts/eigen_solver_factory.py @@ -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") @@ -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")