From 45ef66954b05fd2de1b0dbce55e8559260bf26d6 Mon Sep 17 00:00:00 2001 From: qaumann Date: Tue, 7 Apr 2020 13:57:38 +0200 Subject: [PATCH 1/4] add default parameters depending on problem type --- .../custom_solvers/feast_eigensystem_solver.h | 25 ++++++------------- 1 file changed, 7 insertions(+), 18 deletions(-) diff --git a/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h b/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h index f3807971da18..d9dd74589b98 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();} @@ -135,11 +129,6 @@ 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, "subspace_size" : 0, "max_iteration" : 20, "tolerance" : 1e-12, From d6239e921cde7a3d586211e07ab37f7c30876a74 Mon Sep 17 00:00:00 2001 From: qaumann Date: Tue, 7 Apr 2020 13:58:05 +0200 Subject: [PATCH 2/4] no specified symmetry defaults to symmetric --- .../tests/test_feast_eigensystem_solver.py | 1 - kratos/python_scripts/eigen_solver_factory.py | 12 ++++++++++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/applications/EigenSolversApplication/tests/test_feast_eigensystem_solver.py b/applications/EigenSolversApplication/tests/test_feast_eigensystem_solver.py index f0b3300cc8d7..bf3f44fb6057 100644 --- a/applications/EigenSolversApplication/tests/test_feast_eigensystem_solver.py +++ b/applications/EigenSolversApplication/tests/test_feast_eigensystem_solver.py @@ -16,7 +16,6 @@ def test_real_symmetric_gev(self): settings = KratosMultiphysics.Parameters('''{ "solver_type": "eigen_feast", - "symmetric": true, "number_of_eigenvalues": 3, "search_lowest_eigenvalues": true, "e_min": 0.0, diff --git a/kratos/python_scripts/eigen_solver_factory.py b/kratos/python_scripts/eigen_solver_factory.py index 7a3e91af6963..6e0452021a8f 100644 --- a/kratos/python_scripts/eigen_solver_factory.py +++ b/kratos/python_scripts/eigen_solver_factory.py @@ -23,7 +23,11 @@ def ConstructSolver(settings): if kratos_utils.CheckIfApplicationsAvailable("EigenSolversApplication"): import KratosMultiphysics.EigenSolversApplication as EiSA if EiSA.HasFEAST(): - if settings.Has("symmetric") and settings["symmetric"].GetBool(): + if settings.Has("symmetric"): + is_symmetric = settings["symmetric"].GetBool() + else: + is_symmetric = True + if is_symmetric: eigen_solver = EiSA.FEASTSymmetricEigensystemSolver(settings) else: eigen_solver = EiSA.FEASTGeneralEigensystemSolver(settings) @@ -36,7 +40,11 @@ def ConstructSolver(settings): if kratos_utils.CheckIfApplicationsAvailable("EigenSolversApplication"): import KratosMultiphysics.EigenSolversApplication as EiSA if EiSA.HasFEAST(): - if settings.Has("symmetric") and settings["symmetric"].GetBool(): + if settings.Has("symmetric"): + is_symmetric = settings["symmetric"].GetBool() + else: + is_symmetric = True + if is_symmetric: eigen_solver = EiSA.ComplexFEASTSymmetricEigensystemSolver(settings) else: eigen_solver = EiSA.ComplexFEASTGeneralEigensystemSolver(settings) From cf1cae7672c10cf6677aa14969f1f09543fe6370 Mon Sep 17 00:00:00 2001 From: qaumann Date: Thu, 16 Apr 2020 15:31:29 +0200 Subject: [PATCH 3/4] add function to sort eigensolver result --- .../custom_solvers/feast_eigensystem_solver.h | 98 +++++++++++++++++++ .../tests/test_feast_eigensystem_solver.py | 14 ++- 2 files changed, 108 insertions(+), 4 deletions(-) diff --git a/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h b/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h index 29f747ad3f3a..75763fe69764 100644 --- a/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h +++ b/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h @@ -86,6 +86,93 @@ 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](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 template + void SortingHelper>::SortEigenvalues(VectorType &rEigenvalues, MatrixType &rEigenvectors) + { + std::vector 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 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(); } @@ -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(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) From 30b50c0084a9f045920e6669c9bc0375ba70e4e4 Mon Sep 17 00:00:00 2001 From: qaumann Date: Thu, 30 Apr 2020 16:51:50 +0200 Subject: [PATCH 4/4] style review comments --- .../custom_solvers/feast_eigensystem_solver.h | 37 +++++++++---------- kratos/python_scripts/eigen_solver_factory.py | 20 ++-------- 2 files changed, 22 insertions(+), 35 deletions(-) diff --git a/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h b/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h index 75763fe69764..bdc5b4f0fd85 100644 --- a/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h +++ b/applications/EigenSolversApplication/custom_solvers/feast_eigensystem_solver.h @@ -104,21 +104,21 @@ namespace { // helpers namespace 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::vector 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];}); + [&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](size_t i1, size_t i2) {return std::abs(rEigenvalues[i1]) < std::abs(rEigenvalues[i2]);}); + [&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](size_t i1, size_t i2) {return rEigenvalues[i1] > rEigenvalues[i2];}); + [&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](size_t i1, size_t i2) {return std::abs(rEigenvalues[i1]) > std::abs(rEigenvalues[i2]);}); + [&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; } @@ -126,7 +126,7 @@ namespace { // helpers namespace VectorType tmp_eigenvalues(rEigenvalues.size()); MatrixType tmp_eigenvectors(rEigenvectors.size1(), rEigenvectors.size2()); - for( size_t i=0; i template void SortingHelper>::SortEigenvalues(VectorType &rEigenvalues, MatrixType &rEigenvectors) { - std::vector idx(rEigenvalues.size()); + std::vector 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]);}); + [&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](size_t i1, size_t i2) {return std::abs(rEigenvalues[i1]) < std::abs(rEigenvalues[i2]);}); + [&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](size_t i1, size_t i2) {return std::imag(rEigenvalues[i1]) < std::imag(rEigenvalues[i2]);}); + [&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](size_t i1, size_t i2) {return std::real(rEigenvalues[i1]) > std::real(rEigenvalues[i2]);}); + [&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](size_t i1, size_t i2) {return std::abs(rEigenvalues[i1]) > std::abs(rEigenvalues[i2]);}); + [&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](size_t i1, size_t i2) {return std::imag(rEigenvalues[i1]) > std::imag(rEigenvalues[i2]);}); + [&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; } @@ -166,7 +165,7 @@ namespace { // helpers namespace VectorType tmp_eigenvalues(rEigenvalues.size()); MatrixType tmp_eigenvectors(rEigenvectors.size1(), rEigenvectors.size2()); - for( size_t i=0; i(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 diff --git a/kratos/python_scripts/eigen_solver_factory.py b/kratos/python_scripts/eigen_solver_factory.py index e5b9c70ddbcd..d6d17df9e977 100644 --- a/kratos/python_scripts/eigen_solver_factory.py +++ b/kratos/python_scripts/eigen_solver_factory.py @@ -23,14 +23,8 @@ def ConstructSolver(settings): if kratos_utils.CheckIfApplicationsAvailable("EigenSolversApplication"): import KratosMultiphysics.EigenSolversApplication as EiSA if EiSA.HasFEAST(): - if settings.Has("symmetric"): - is_symmetric = settings["symmetric"].GetBool() - else: - is_symmetric = True - if is_symmetric: - 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") @@ -40,14 +34,8 @@ def ConstructSolver(settings): if kratos_utils.CheckIfApplicationsAvailable("EigenSolversApplication"): import KratosMultiphysics.EigenSolversApplication as EiSA if EiSA.HasFEAST(): - if settings.Has("symmetric"): - is_symmetric = settings["symmetric"].GetBool() - else: - is_symmetric = True - if is_symmetric: - 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")