diff --git a/include/htool/solvers/ddm.hpp b/include/htool/solvers/ddm.hpp index e2a10b40..3a6d843b 100644 --- a/include/htool/solvers/ddm.hpp +++ b/include/htool/solvers/ddm.hpp @@ -487,7 +487,7 @@ DDM make_DDM_solver_w_custom_local template > DDM make_DDM_solver_w_custom_local_solver(const DistributedOperator &distributed_operator, HMatrix &local_hmatrix, Matrix &B, Matrix &C, Matrix &D, const std::vector &neighbors, const std::vector> &intersections) { - int n = local_hmatrix.get_target_cluster().get_size() + C.nb_rows(); + int n = local_hmatrix.get_target_cluster().get_size() + D.nb_rows(); // Timing double mytime, maxtime; diff --git a/include/htool/solvers/geneo/coarse_space_builder.hpp b/include/htool/solvers/geneo/coarse_space_builder.hpp index 56a72798..125d9b72 100644 --- a/include/htool/solvers/geneo/coarse_space_builder.hpp +++ b/include/htool/solvers/geneo/coarse_space_builder.hpp @@ -12,7 +12,7 @@ namespace htool { -template +template > class GeneoCoarseSpaceDenseBuilder : public VirtualCoarseSpaceBuilder { protected: @@ -24,16 +24,42 @@ class GeneoCoarseSpaceDenseBuilder : public VirtualCoarseSpaceBuilder m_geneo_threshold = -1.; - explicit GeneoCoarseSpaceDenseBuilder(int size_wo_overlap, const Matrix &Ai, Matrix &Bi, char symmetry, char uplo, int geneo_nu, htool::underlying_type geneo_threshold) : m_size_wo_overlap(size_wo_overlap), m_size_with_overlap(Ai.nb_cols()), m_DAiD(m_size_with_overlap, m_size_with_overlap), m_Bi(Bi), m_symmetry(symmetry), m_uplo(uplo), m_geneo_nu(geneo_nu), m_geneo_threshold(geneo_threshold) { - for (int i = 0; i < m_size_wo_overlap; i++) { - std::copy_n(Ai.data() + i * m_size_with_overlap, m_size_wo_overlap, &(m_DAiD(0, i))); + explicit GeneoCoarseSpaceDenseBuilder(int size_wo_overlap, int size_with_overlap, const Matrix &Ai, const Matrix &Bi, char symmetry, char uplo, int geneo_nu, htool::underlying_type geneo_threshold) : m_size_wo_overlap(size_wo_overlap), m_size_with_overlap(size_with_overlap), m_DAiD(m_size_with_overlap, m_size_with_overlap), m_Bi(Bi), m_symmetry(symmetry), m_uplo(uplo), m_geneo_nu(geneo_nu), m_geneo_threshold(geneo_threshold) { + if (Ai.nb_cols() == size_with_overlap) { + for (int i = 0; i < m_size_wo_overlap; i++) { + std::copy_n(Ai.data() + i * m_size_with_overlap, m_size_wo_overlap, &(m_DAiD(0, i))); + } + } else if (Ai.nb_cols() == size_wo_overlap) { + for (int i = 0; i < m_size_wo_overlap; i++) { + std::copy_n(Ai.data() + i * m_size_wo_overlap, m_size_wo_overlap, &(m_DAiD(0, i))); + } } } public: - static GeneoCoarseSpaceDenseBuilder GeneoWithNu(int size_wo_overlap, const Matrix &Ai, Matrix &Bi, char symmetry, char uplo, int geneo_nu) { return GeneoCoarseSpaceDenseBuilder{size_wo_overlap, Ai, Bi, symmetry, uplo, geneo_nu, -1}; } + static GeneoCoarseSpaceDenseBuilder GeneoWithNu(int size_wo_overlap, int size_with_overlap, const Matrix &Ai, const Matrix &Bi, char symmetry, char uplo, int geneo_nu) { return GeneoCoarseSpaceDenseBuilder{size_wo_overlap, size_with_overlap, Ai, Bi, symmetry, uplo, geneo_nu, -1}; } + + static GeneoCoarseSpaceDenseBuilder GeneoWithNu(int size_wo_overlap, int size_with_overlap, const HMatrix &Ai, const Matrix &Bi, char symmetry, char uplo, int geneo_nu) { + Matrix temp(Ai.nb_rows(), Ai.nb_cols()); + if (Ai.nb_rows() == size_with_overlap) { + copy_to_dense_in_user_numbering(Ai, temp.data()); + } else if (Ai.nb_rows() == size_wo_overlap) { + copy_to_dense(Ai, temp.data()); + } + return GeneoCoarseSpaceDenseBuilder{size_wo_overlap, size_with_overlap, temp, Bi, symmetry, uplo, geneo_nu, -1}; + } + + static GeneoCoarseSpaceDenseBuilder GeneoWithThreshold(int size_wo_overlap, int size_with_overlap, const Matrix &Ai, const Matrix &Bi, char symmetry, char uplo, htool::underlying_type geneo_threshold) { return GeneoCoarseSpaceDenseBuilder{size_wo_overlap, size_with_overlap, Ai, Bi, symmetry, uplo, 0, geneo_threshold}; } - static GeneoCoarseSpaceDenseBuilder GeneoWithThreshold(int size_wo_overlap, const Matrix &Ai, Matrix &Bi, char symmetry, char uplo, htool::underlying_type geneo_threshold) { return GeneoCoarseSpaceDenseBuilder{size_wo_overlap, Ai, Bi, symmetry, uplo, 0, geneo_threshold}; } + static GeneoCoarseSpaceDenseBuilder GeneoWithThreshold(int size_wo_overlap, int size_with_overlap, const HMatrix &Ai, const Matrix &Bi, char symmetry, char uplo, htool::underlying_type geneo_threshold) { + Matrix temp(Ai.nb_rows(), Ai.nb_cols()); + if (Ai.nb_rows() == size_with_overlap) { + copy_to_dense_in_user_numbering(Ai, temp.data()); + } else if (Ai.nb_rows() == size_wo_overlap) { + copy_to_dense(Ai, temp.data()); + } + return GeneoCoarseSpaceDenseBuilder{size_wo_overlap, size_with_overlap, temp, Bi, symmetry, uplo, 0, geneo_threshold}; + } virtual Matrix build_coarse_space() override { diff --git a/include/htool/solvers/local_solvers/local_hmatrix_plus_overlap_solvers.hpp b/include/htool/solvers/local_solvers/local_hmatrix_plus_overlap_solvers.hpp index 0cbe2abe..803604b9 100644 --- a/include/htool/solvers/local_solvers/local_hmatrix_plus_overlap_solvers.hpp +++ b/include/htool/solvers/local_solvers/local_hmatrix_plus_overlap_solvers.hpp @@ -32,16 +32,20 @@ class LocalHMatrixPlusOverlapSolver : public VirtualLocalSolver 0) { + if (m_local_hmatrix.get_UPLO() == 'L' && m_C.nb_rows() > 0) { triangular_hmatrix_matrix_solve('R', m_local_hmatrix.get_UPLO(), is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(1), m_local_hmatrix, m_C); add_matrix_matrix_product('N', is_complex() ? 'C' : 'T', CoefficientPrecision(-1), m_C, m_C, CoefficientPrecision(1), m_D); cholesky_factorization(m_local_hmatrix.get_UPLO(), m_D); + } else if (m_local_hmatrix.get_UPLO() == 'U' && m_B.nb_cols() > 0) { + triangular_hmatrix_matrix_solve('L', m_local_hmatrix.get_UPLO(), is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(1), m_local_hmatrix, m_B); + add_matrix_matrix_product(is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(-1), m_B, m_B, CoefficientPrecision(1), m_D); + cholesky_factorization(m_local_hmatrix.get_UPLO(), m_D); } } } void solve(CoefficientPrecision *const b, const unsigned short &mu = 1) const { int local_size_wo_overlap = m_local_hmatrix.get_target_cluster().get_size(); - int size_overlap = m_C.nb_rows(); + int size_overlap = m_D.nb_rows(); int local_size_w_overlap = local_size_wo_overlap + size_overlap; Matrix b1(local_size_wo_overlap, mu), b2(size_overlap, mu); @@ -50,30 +54,7 @@ class LocalHMatrixPlusOverlapSolver : public VirtualLocalSolver 0) { - add_matrix_matrix_product('N', 'N', CoefficientPrecision(-1), m_C, b1, CoefficientPrecision(1), b2); - triangular_matrix_matrix_solve('L', 'L', 'N', 'U', CoefficientPrecision(1), m_D, b2); - - triangular_matrix_matrix_solve('L', 'U', 'N', 'N', CoefficientPrecision(1), m_D, b2); - add_matrix_matrix_product('N', 'N', CoefficientPrecision(-1), m_B, b2, CoefficientPrecision(1), b1); - } - triangular_hmatrix_matrix_solve('L', 'U', 'N', 'N', CoefficientPrecision(1), m_local_hmatrix, b1); - - } else if (m_local_hmatrix.get_symmetry() == 'S' || m_local_hmatrix.get_symmetry() == 'H') { - triangular_hmatrix_matrix_solve('L', 'L', 'N', 'N', CoefficientPrecision(1), m_local_hmatrix, b1); - - if (m_C.nb_rows() > 0) { - add_matrix_matrix_product('N', 'N', CoefficientPrecision(-1), m_C, b1, CoefficientPrecision(1), b2); - triangular_matrix_matrix_solve('L', 'L', 'N', 'U', CoefficientPrecision(1), m_D, b2); - - triangular_matrix_matrix_solve('L', 'L', is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(1), m_D, b2); - add_matrix_matrix_product('N', 'N', CoefficientPrecision(-1), m_C, b2, CoefficientPrecision(1), b1); - } - triangular_hmatrix_matrix_solve('L', 'L', is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(1), m_local_hmatrix, b1); - } + this->solve(b1, b2); for (int i = 0; i < mu; i++) { std::copy_n(b1.data() + i * local_size_wo_overlap, local_size_wo_overlap, b + i * local_size_w_overlap); @@ -83,7 +64,7 @@ class LocalHMatrixPlusOverlapSolver : public VirtualLocalSolver b1(local_size_wo_overlap, mu), b2(size_overlap, mu); @@ -91,34 +72,48 @@ class LocalHMatrixPlusOverlapSolver : public VirtualLocalSolversolve(b1, b2); + + for (int i = 0; i < mu; i++) { + std::copy_n(b1.data() + i * local_size_wo_overlap, local_size_wo_overlap, x + i * local_size_w_overlap); + std::copy_n(b2.data() + i * size_overlap, size_overlap, x + i * local_size_w_overlap + local_size_wo_overlap); + } + } + + private: + void solve(Matrix &b1, Matrix &b2) const { if (m_local_hmatrix.get_symmetry() == 'N') { triangular_hmatrix_matrix_solve('L', 'L', 'N', 'U', CoefficientPrecision(1), m_local_hmatrix, b1); - if (m_C.nb_rows() > 0) { add_matrix_matrix_product('N', 'N', CoefficientPrecision(-1), m_C, b1, CoefficientPrecision(1), b2); triangular_matrix_matrix_solve('L', 'L', 'N', 'U', CoefficientPrecision(1), m_D, b2); + triangular_matrix_matrix_solve('L', 'U', 'N', 'N', CoefficientPrecision(1), m_D, b2); add_matrix_matrix_product('N', 'N', CoefficientPrecision(-1), m_B, b2, CoefficientPrecision(1), b1); } triangular_hmatrix_matrix_solve('L', 'U', 'N', 'N', CoefficientPrecision(1), m_local_hmatrix, b1); - - } else if (m_local_hmatrix.get_symmetry() == 'S' || m_local_hmatrix.get_symmetry() == 'H') { + } else if ((m_local_hmatrix.get_symmetry() == 'S' || m_local_hmatrix.get_symmetry() == 'H') && m_local_hmatrix.get_UPLO() == 'L') { triangular_hmatrix_matrix_solve('L', 'L', 'N', 'N', CoefficientPrecision(1), m_local_hmatrix, b1); + if (m_C.nb_rows() > 0) { add_matrix_matrix_product('N', 'N', CoefficientPrecision(-1), m_C, b1, CoefficientPrecision(1), b2); triangular_matrix_matrix_solve('L', 'L', 'N', 'U', CoefficientPrecision(1), m_D, b2); - triangular_matrix_matrix_solve('L', 'L', is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(1), m_D, b2); add_matrix_matrix_product(is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(-1), m_C, b2, CoefficientPrecision(1), b1); } triangular_hmatrix_matrix_solve('L', 'L', is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(1), m_local_hmatrix, b1); - } - // std::cout << "TEST " << get_max(b1) << " " << get_max(b2) << "\n"; - for (int i = 0; i < mu; i++) { - std::copy_n(b1.data() + i * local_size_wo_overlap, local_size_wo_overlap, x + i * local_size_w_overlap); - std::copy_n(b2.data() + i * size_overlap, size_overlap, x + i * local_size_w_overlap + local_size_wo_overlap); + } else if ((m_local_hmatrix.get_symmetry() == 'S' || m_local_hmatrix.get_symmetry() == 'H') && m_local_hmatrix.get_UPLO() == 'U') { + triangular_hmatrix_matrix_solve('L', 'U', is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(1), m_local_hmatrix, b1); + + if (m_B.nb_cols() > 0) { + add_matrix_matrix_product(is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(-1), m_B, b1, CoefficientPrecision(1), b2); + triangular_matrix_matrix_solve('L', 'U', 'N', 'U', CoefficientPrecision(1), m_D, b2); + add_matrix_matrix_product('N', 'N', CoefficientPrecision(-1), m_B, b2, CoefficientPrecision(1), b1); + } + + triangular_hmatrix_matrix_solve('L', 'U', 'N', 'N', CoefficientPrecision(1), m_local_hmatrix, b1); } } }; diff --git a/include/htool/solvers/utility.hpp b/include/htool/solvers/utility.hpp index 222edb0d..78dd87b0 100644 --- a/include/htool/solvers/utility.hpp +++ b/include/htool/solvers/utility.hpp @@ -219,16 +219,18 @@ class DDMSolverBuilder { // Overlap D.resize(local_size_with_overlap - local_size_without_overlap, local_size_with_overlap - local_size_without_overlap); - C.resize(local_size_with_overlap - local_size_without_overlap, local_size_without_overlap); - - generator0.copy_submatrix(local_size_with_overlap - local_size_without_overlap, local_size_without_overlap, overlap_num.data(), inside_num.data(), C.data()); generator0.copy_submatrix(local_size_with_overlap - local_size_without_overlap, local_size_with_overlap - local_size_without_overlap, overlap_num.data(), overlap_num.data(), D.data()); bool sym = (distributed_operator0.get_symmetry_type() == 'S' || (distributed_operator0.get_symmetry_type() == 'H' && is_complex())) ? true : false; - if (!sym) { + + if (!sym or block_diagonal_hmatrix0.get_UPLO() == 'U') { B.resize(local_size_without_overlap, local_size_with_overlap - local_size_without_overlap); generator0.copy_submatrix(local_size_without_overlap, local_size_with_overlap - local_size_without_overlap, inside_num.data(), overlap_num.data(), B.data()); + } else if (!sym or block_diagonal_hmatrix0.get_UPLO() == 'L') { + C.resize(local_size_with_overlap - local_size_without_overlap, local_size_without_overlap); + + generator0.copy_submatrix(local_size_with_overlap - local_size_without_overlap, local_size_without_overlap, overlap_num.data(), inside_num.data(), C.data()); } return blocks_in_overlap0; }; @@ -250,7 +252,7 @@ class DDMSolverBuilder { std::unique_ptr> local_cluster; - std::function(const VirtualGeneratorInUserNumbering &, underlying_type, CoordinatePrecision, char)> initialize_local_hmatrix = [this](const VirtualGeneratorInUserNumbering &generator0, underlying_type epsilon0, CoordinatePrecision eta0, char symmetry0) { + std::function(const VirtualGeneratorInUserNumbering &, underlying_type, CoordinatePrecision, char, char)> initialize_local_hmatrix = [this](const VirtualGeneratorInUserNumbering &generator0, underlying_type epsilon0, CoordinatePrecision eta0, char symmetry0, char UPLO0) { struct LocalGeneratorInUserNumbering : public VirtualGeneratorInUserNumbering { const std::vector &m_target_local_to_global_numbering; const std::vector &m_source_local_to_global_numbering; @@ -274,7 +276,7 @@ class DDMSolverBuilder { LocalGeneratorInUserNumbering local_generator(generator0, local_to_global_numbering, local_to_global_numbering); // Local HMatrix - HMatrixTreeBuilder local_hmatrix_builder(*local_cluster, *local_cluster, epsilon0, eta0, symmetry0, symmetry0 != 'N' ? 'L' : 'N', -1, -1, -1); + HMatrixTreeBuilder local_hmatrix_builder(*local_cluster, *local_cluster, epsilon0, eta0, symmetry0, UPLO0, -1, -1, -1); return local_hmatrix_builder.build(local_generator); }; @@ -291,7 +293,7 @@ class DDMSolverBuilder { DDMSolverBuilder(DistributedOperator &distributed_operator, HMatrix &block_diagonal_hmatrix, const VirtualGeneratorInUserNumbering &generator, const std::vector &ovr_subdomain_to_global, const std::vector &cluster_to_ovr_subdomain, const std::vector &neighbors, const std::vector> &intersections) : m_local_numbering(ovr_subdomain_to_global, cluster_to_ovr_subdomain, intersections), local_to_global_numbering(m_local_numbering.local_to_global_numbering), blocks_in_overlap(initialize_blocks_in_overlap(distributed_operator, block_diagonal_hmatrix, generator)), solver(make_DDM_solver_w_custom_local_solver(distributed_operator, block_diagonal_hmatrix, blocks_in_overlap[0], blocks_in_overlap[1], blocks_in_overlap[2], neighbors, m_local_numbering.intersections)) {} // DDM building local hmatrix with overlap - DDMSolverBuilder(DistributedOperator &distributed_operator, const std::vector &ovr_subdomain_to_global, const std::vector &cluster_to_ovr_subdomain, const std::vector &neighbors, const std::vector> &intersections, const VirtualGeneratorInUserNumbering &generator, int spatial_dimension, const CoordinatePrecision *global_geometry, underlying_type epsilon, CoordinatePrecision eta) : m_local_numbering(ovr_subdomain_to_global, cluster_to_ovr_subdomain, intersections), local_to_global_numbering(m_local_numbering.local_to_global_numbering), local_cluster(std::make_unique>(initialize_local_cluster(spatial_dimension, global_geometry))), local_hmatrix(std::make_unique>(initialize_local_hmatrix(generator, epsilon, eta, distributed_operator.get_symmetry_type()))), solver(make_DDM_solver_w_custom_local_solver(distributed_operator, *local_hmatrix, neighbors, m_local_numbering.intersections, true)) {} + DDMSolverBuilder(DistributedOperator &distributed_operator, const std::vector &ovr_subdomain_to_global, const std::vector &cluster_to_ovr_subdomain, const std::vector &neighbors, const std::vector> &intersections, const VirtualGeneratorInUserNumbering &generator, int spatial_dimension, const CoordinatePrecision *global_geometry, underlying_type epsilon, CoordinatePrecision eta) : m_local_numbering(ovr_subdomain_to_global, cluster_to_ovr_subdomain, intersections), local_to_global_numbering(m_local_numbering.local_to_global_numbering), local_cluster(std::make_unique>(initialize_local_cluster(spatial_dimension, global_geometry))), local_hmatrix(std::make_unique>(initialize_local_hmatrix(generator, epsilon, eta, distributed_operator.get_symmetry_type(), distributed_operator.get_storage_type()))), solver(make_DDM_solver_w_custom_local_solver(distributed_operator, *local_hmatrix, neighbors, m_local_numbering.intersections, true)) {} }; } // namespace htool diff --git a/tests/functional_tests/solvers/test_solver.cpp b/tests/functional_tests/solvers/test_solver.cpp index b9a53ad2..6508ee04 100644 --- a/tests/functional_tests/solvers/test_solver.cpp +++ b/tests/functional_tests/solvers/test_solver.cpp @@ -33,19 +33,28 @@ int main(int argc, char *argv[]) { } for (auto symmetry : symmetries) { - is_error = is_error || test_solver_wo_overlap, double>>(argc, argv, nb_rhs, symmetry, datapath_final); + is_error = is_error || test_solver_wo_overlap, double>>(argc, argv, nb_rhs, symmetry, symmetry == 'N' ? 'N' : 'L', datapath_final); is_error = is_error || test_solver_ddm_adding_overlap>>(argc, argv, nb_rhs, data_symmetry, symmetry, datapath_final); - is_error = is_error || test_solver_ddm>>(argc, argv, nb_rhs, data_symmetry, symmetry, datapath_final); - is_error = is_error || test_solver_wo_overlap, double>>(argc, argv, nb_rhs, symmetry, datapath_final); - test_solver_ddm_adding_overlap>>(argc, argv, nb_rhs, data_symmetry, symmetry, datapath_final); - is_error = is_error || test_solver_ddm>>(argc, argv, nb_rhs, data_symmetry, symmetry, datapath_final); + is_error = is_error || test_solver_ddm>>(argc, argv, nb_rhs, data_symmetry, symmetry, symmetry == 'N' ? 'N' : 'L', datapath_final); + + std::vector storage; + if (symmetry == 'N') { + storage.push_back('N'); + } else { + storage.push_back('U'); + storage.push_back('L'); + } + for (auto UPLO : storage) { + is_error = is_error || test_solver_wo_overlap, double>>(argc, argv, nb_rhs, symmetry, UPLO, datapath_final); + test_solver_ddm_adding_overlap(argc, argv, nb_rhs, data_symmetry, symmetry, UPLO, datapath_final); + is_error = is_error || test_solver_ddm>>(argc, argv, nb_rhs, data_symmetry, symmetry, UPLO, datapath_final); + } } } } // Finalize the MPI environment. MPI_Finalize(); - if (is_error) { return 1; } diff --git a/tests/functional_tests/solvers/test_solver_ddm.hpp b/tests/functional_tests/solvers/test_solver_ddm.hpp index 4a03d9a2..2af540b5 100644 --- a/tests/functional_tests/solvers/test_solver_ddm.hpp +++ b/tests/functional_tests/solvers/test_solver_ddm.hpp @@ -24,7 +24,7 @@ using namespace std; using namespace htool; template -int test_solver_ddm(int argc, char *argv[], int mu, char data_symmetry, char symmetric, std::string datapath) { +int test_solver_ddm(int argc, char *argv[], int mu, char data_symmetry, char symmetric, char UPLO, std::string datapath) { // Get the number of processes int size; @@ -50,11 +50,6 @@ int test_solver_ddm(int argc, char *argv[], int mu, char data_symmetry, char sym // epsilon = 1e-10; double eta = 10; - char UPLO = 'N'; - if (symmetric != 'N') { - UPLO = 'L'; - } - // Clustering if (rank == 0) std::cout << "Creating cluster tree" << std::endl; @@ -147,7 +142,8 @@ int test_solver_ddm(int argc, char *argv[], int mu, char data_symmetry, char sym opt.parse("-hpddm_geneo_nu 2"); Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); - auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithNu(local_size_wo_overlap, default_ddm_solver.block_diagonal_dense_matrix, Ki, symmetric, UPLO, 2); + int local_size_with_overlap = Ki.nb_cols(); + auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithNu(local_size_wo_overlap, local_size_with_overlap, *default_ddm_solver.local_hmatrix, Ki, symmetric, UPLO, 2); // geneo_coarse_space_dense_builder.set_geneo_nu(2); GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); ddm_with_overlap.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); @@ -224,7 +220,8 @@ int test_solver_ddm(int argc, char *argv[], int mu, char data_symmetry, char sym auto &ddm_with_overlap_threshold = default_ddm_solver_with_threshold.solver; Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); - auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithThreshold(local_size_wo_overlap, default_ddm_solver.block_diagonal_dense_matrix, Ki, symmetric, UPLO, 100); + int local_size_with_overlap = Ki.nb_cols(); + auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithThreshold(local_size_wo_overlap, local_size_with_overlap, *default_ddm_solver.local_hmatrix, Ki, symmetric, UPLO, 100); // geneo_coarse_space_dense_builder.set_geneo_threshold(100); GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); ddm_with_overlap_threshold.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); @@ -246,7 +243,7 @@ int test_solver_ddm(int argc, char *argv[], int mu, char data_symmetry, char sym } template <> -int test_solver_ddm>>(int argc, char *argv[], int mu, char, char symmetric, std::string datapath) { +int test_solver_ddm>>(int argc, char *argv[], int mu, char data_symmetry, char symmetric, char UPLO, std::string datapath) { // Get the number of processes int size; @@ -272,11 +269,6 @@ int test_solver_ddm>>(int argc, char *argv[], i // epsilon = 1e-10; double eta = 10; - char UPLO = 'N'; - if (symmetric != 'N') { - UPLO = 'L'; - } - // Clustering if (rank == 0) std::cout << "Creating cluster tree" << std::endl; @@ -364,16 +356,17 @@ int test_solver_ddm>>(int argc, char *argv[], i MPI_Barrier(MPI_COMM_WORLD); opt.parse("-hpddm_schwarz_method asm "); Matrix> Ki; - // if (data_symmetry == 'S' && size > 1) { - // opt.remove("geneo_threshold"); - // opt.parse("-hpddm_geneo_nu 2"); - // Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); - // int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); - // auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithNu(local_size_wo_overlap, default_ddm_solver.block_diagonal_dense_matrix, Ki, symmetric, UPLO, 2); - // // geneo_coarse_space_dense_builder.set_geneo_nu(2); - // GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); - // ddm_with_overlap.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); - // } + if (data_symmetry == 'S' && size > 1) { + opt.remove("geneo_threshold"); + opt.parse("-hpddm_geneo_nu 2"); + Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); + int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); + int local_size_with_overlap = Ki.nb_cols(); + auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithNu(local_size_wo_overlap, local_size_with_overlap, *default_ddm_solver.local_hmatrix, Ki, symmetric, UPLO, 2); + // geneo_coarse_space_dense_builder.set_geneo_nu(2); + GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); + ddm_with_overlap.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); + } ddm_with_overlap.facto_one_level(); ddm_with_overlap.solve(f_global.data(), x_global.data(), mu); @@ -405,64 +398,65 @@ int test_solver_ddm>>(int argc, char *argv[], i x_global = 0; // DDM two level ASM with overlap - // if (data_symmetry == 'S' && size > 1) { - // if (rank == 0) - // std::cout << "ASM two level with overlap:" << std::endl; - // MPI_Barrier(MPI_COMM_WORLD); - // opt.parse("-hpddm_schwarz_method asm -hpddm_schwarz_coarse_correction additive"); - // ddm_with_overlap.solve(f_global.data(), x_global.data(), mu); - // ddm_with_overlap.print_infos(); - // error2 = normFrob(f_global - A * x_global) / normFrob(f_global); - // if (rank == 0) { - // cout << "error: " << error2 << endl; - // } - - // test = test || !(error2 < tol); - - // x_global = 0; - - // if (rank == 0) - // std::cout << "RAS two level with overlap:" << std::endl; - // MPI_Barrier(MPI_COMM_WORLD); - // opt.parse("-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction additive"); - // ddm_with_overlap.solve(f_global.data(), x_global.data(), mu); - // ddm_with_overlap.print_infos(); - // ddm_with_overlap.clean(); - // error2 = normFrob(f_global - A * x_global) / normFrob(f_global); - // if (rank == 0) { - // cout << "error: " << error2 << endl; - // } - - // test = test || !(error2 < tol); - - // x_global = 0; - - // // DDM solver with threshold - // if (rank == 0) - // std::cout << "RAS two level with overlap and threshold:" << std::endl; - // opt.parse("-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction additive"); - // // DDM> ddm_with_overlap_threshold(Generator, &Operator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections); - // solver_builder default_ddm_solver_with_threshold(Operator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections, Generator, 3, geometry.data(), epsilon, eta); - // auto &ddm_with_overlap_threshold = default_ddm_solver_with_threshold.solver; - // Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); - // int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); - // auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithThreshold(local_size_wo_overlap, default_ddm_solver.block_diagonal_dense_matrix, Ki, symmetric, UPLO, 100); - // // geneo_coarse_space_dense_builder.set_geneo_threshold(100); - // GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); - // ddm_with_overlap_threshold.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); - // ddm_with_overlap_threshold.facto_one_level(); - // ddm_with_overlap_threshold.solve(f_global.data(), x_global.data(), mu); - // ddm_with_overlap_threshold.print_infos(); - // ddm_with_overlap_threshold.clean(); - // error2 = normFrob(f_global - A * x_global) / normFrob(f_global); - // if (rank == 0) { - // cout << "error: " << error2 << endl; - // } - - // test = test || !(error2 < tol); - - // x_global = 0; - // } + if (data_symmetry == 'S' && size > 1) { + if (rank == 0) + std::cout << "ASM two level with overlap:" << std::endl; + MPI_Barrier(MPI_COMM_WORLD); + opt.parse("-hpddm_schwarz_method asm -hpddm_schwarz_coarse_correction additive"); + ddm_with_overlap.solve(f_global.data(), x_global.data(), mu); + ddm_with_overlap.print_infos(); + error2 = normFrob(f_global - A * x_global) / normFrob(f_global); + if (rank == 0) { + cout << "error: " << error2 << endl; + } + + test = test || !(error2 < tol); + + x_global = 0; + + if (rank == 0) + std::cout << "RAS two level with overlap:" << std::endl; + MPI_Barrier(MPI_COMM_WORLD); + opt.parse("-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction additive"); + ddm_with_overlap.solve(f_global.data(), x_global.data(), mu); + ddm_with_overlap.print_infos(); + ddm_with_overlap.clean(); + error2 = normFrob(f_global - A * x_global) / normFrob(f_global); + if (rank == 0) { + cout << "error: " << error2 << endl; + } + + test = test || !(error2 < tol); + + x_global = 0; + + // DDM solver with threshold + if (rank == 0) + std::cout << "RAS two level with overlap and threshold:" << std::endl; + opt.parse("-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction additive"); + // DDM> ddm_with_overlap_threshold(Generator, &Operator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections); + DDMSolverBuilder> default_ddm_solver_with_threshold(Operator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections, Generator, 3, geometry.data(), epsilon, eta); + auto &ddm_with_overlap_threshold = default_ddm_solver_with_threshold.solver; + Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); + int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); + int local_size_with_overlap = Ki.nb_cols(); + auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithThreshold(local_size_wo_overlap, local_size_with_overlap, *default_ddm_solver.local_hmatrix, Ki, symmetric, UPLO, 100); + // geneo_coarse_space_dense_builder.set_geneo_threshold(100); + GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); + ddm_with_overlap_threshold.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); + ddm_with_overlap_threshold.facto_one_level(); + ddm_with_overlap_threshold.solve(f_global.data(), x_global.data(), mu); + ddm_with_overlap_threshold.print_infos(); + ddm_with_overlap_threshold.clean(); + error2 = normFrob(f_global - A * x_global) / normFrob(f_global); + if (rank == 0) { + cout << "error: " << error2 << endl; + } + + test = test || !(error2 < tol); + + x_global = 0; + } return test; } diff --git a/tests/functional_tests/solvers/test_solver_ddm_adding_overlap.hpp b/tests/functional_tests/solvers/test_solver_ddm_adding_overlap.hpp index 76389dc6..2dba05c5 100644 --- a/tests/functional_tests/solvers/test_solver_ddm_adding_overlap.hpp +++ b/tests/functional_tests/solvers/test_solver_ddm_adding_overlap.hpp @@ -147,7 +147,8 @@ int test_solver_ddm_adding_overlap(int argc, char *argv[], int mu, char data_sym opt.parse("-hpddm_geneo_nu 2"); Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); - auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithNu(local_size_wo_overlap, default_ddm_solver.block_diagonal_dense_matrix, Ki, symmetric, UPLO, 2); + int local_size_with_overlap = Ki.nb_cols(); + auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithNu(local_size_wo_overlap, local_size_with_overlap, default_ddm_solver.block_diagonal_dense_matrix, Ki, symmetric, UPLO, 2); // geneo_coarse_space_dense_builder.set_geneo_nu(4); GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); ddm_with_overlap.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); @@ -225,7 +226,8 @@ int test_solver_ddm_adding_overlap(int argc, char *argv[], int mu, char data_sym Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); - auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithThreshold(local_size_wo_overlap, default_ddm_solver.block_diagonal_dense_matrix, Ki, symmetric, UPLO, 100); + int local_size_with_overlap = Ki.nb_cols(); + auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithThreshold(local_size_wo_overlap, local_size_with_overlap, default_ddm_solver.block_diagonal_dense_matrix, Ki, symmetric, UPLO, 100); // geneo_coarse_space_dense_builder.set_geneo_threshold(100); GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); ddm_with_overlap_threshold.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); @@ -246,8 +248,7 @@ int test_solver_ddm_adding_overlap(int argc, char *argv[], int mu, char data_sym return test; } -template <> -int test_solver_ddm_adding_overlap>>(int argc, char *argv[], int mu, char, char symmetric, std::string datapath) { +int test_solver_ddm_adding_overlap(int argc, char *argv[], int mu, char data_symmetry, char symmetric, char UPLO, std::string datapath) { // Get the number of processes int size; @@ -272,11 +273,6 @@ int test_solver_ddm_adding_overlap>>(int a double epsilon = tol; double eta = 10; - char UPLO = 'N'; - if (symmetric != 'N') { - UPLO = 'L'; - } - // Clustering if (rank == 0) std::cout << "Creating cluster tree" << std::endl; @@ -363,16 +359,17 @@ int test_solver_ddm_adding_overlap>>(int a MPI_Barrier(MPI_COMM_WORLD); opt.parse("-hpddm_schwarz_method asm "); Matrix> Ki; - // if (data_symmetry == 'S' && size > 1) { - // opt.remove("geneo_threshold"); - // opt.parse("-hpddm_geneo_nu 2"); - // Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); - // int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); - // auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithNu(local_size_wo_overlap, default_ddm_solver.block_diagonal_dense_matrix, Ki, symmetric, UPLO, 2); - // // geneo_coarse_space_dense_builder.set_geneo_nu(4); - // GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); - // ddm_with_overlap.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); - // } + if (data_symmetry == 'S' && size > 1) { + opt.remove("geneo_threshold"); + opt.parse("-hpddm_geneo_nu 2"); + Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); + int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); + int local_size_with_overlap = Ki.nb_cols(); + auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithNu(local_size_wo_overlap, local_size_with_overlap, *default_build.block_diagonal_hmatrix, Ki, symmetric, UPLO, 2); + // geneo_coarse_space_dense_builder.set_geneo_nu(4); + GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); + ddm_with_overlap.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); + } ddm_with_overlap.facto_one_level(); ddm_with_overlap.solve(f_global.data(), x_global.data(), mu); @@ -404,65 +401,67 @@ int test_solver_ddm_adding_overlap>>(int a x_global = 0; // DDM two level ASM with overlap - // if (data_symmetry == 'S' && size > 1) { - // if (rank == 0) - // std::cout << "ASM two level with overlap:" << std::endl; - // MPI_Barrier(MPI_COMM_WORLD); - // opt.parse("-hpddm_schwarz_method asm -hpddm_schwarz_coarse_correction additive"); - // ddm_with_overlap.solve(f_global.data(), x_global.data(), mu); - // ddm_with_overlap.print_infos(); - // error2 = normFrob(f_global - A * x_global) / normFrob(f_global); - // if (rank == 0) { - // cout << "error: " << error2 << endl; - // } - - // test = test || !(error2 < tol); - - // x_global = 0; - - // if (rank == 0) - // std::cout << "RAS two level with overlap:" << std::endl; - // MPI_Barrier(MPI_COMM_WORLD); - // opt.parse("-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction additive"); - // ddm_with_overlap.solve(f_global.data(), x_global.data(), mu); - // ddm_with_overlap.print_infos(); - // ddm_with_overlap.clean(); - // error2 = normFrob(f_global - A * x_global) / normFrob(f_global); - // if (rank == 0) { - // cout << "error: " << error2 << endl; - // } - - // test = test || !(error2 < tol); - - // x_global = 0; - - // // DDM solver with threshold - // if (rank == 0) - // std::cout << "RAS two level with overlap and threshold:" << std::endl; - // opt.parse("-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction additive -hpddm_geneo_threshold 100"); - // // DDM> ddm_with_overlap_threshold(Generator, &Operator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections); - // solver_builder default_ddm_solver_with_threshold(Operator, local_block_diagonal_hmatrix, Generator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections); - // auto &ddm_with_overlap_threshold = default_ddm_solver_with_threshold.solver; - - // Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); - // int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); - // auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithThreshold(local_size_wo_overlap, default_ddm_solver.block_diagonal_dense_matrix, Ki, symmetric, UPLO, 100); - // // geneo_coarse_space_dense_builder.set_geneo_threshold(100); - // GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); - // ddm_with_overlap_threshold.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); - // ddm_with_overlap_threshold.facto_one_level(); - // ddm_with_overlap_threshold.solve(f_global.data(), x_global.data(), mu); - // ddm_with_overlap_threshold.print_infos(); - // ddm_with_overlap_threshold.clean(); - // error2 = normFrob(f_global - A * x_global) / normFrob(f_global); - // if (rank == 0) { - // cout << "error: " << error2 << endl; - // } - - // test = test || !(error2 < tol); - - // x_global = 0; - // } + if (data_symmetry == 'S' && size > 1) { + if (rank == 0) + std::cout << "ASM two level with overlap:" << std::endl; + MPI_Barrier(MPI_COMM_WORLD); + opt.parse("-hpddm_schwarz_method asm -hpddm_schwarz_coarse_correction additive"); + ddm_with_overlap.solve(f_global.data(), x_global.data(), mu); + ddm_with_overlap.print_infos(); + error2 = normFrob(f_global - A * x_global) / normFrob(f_global); + if (rank == 0) { + cout << "error: " << error2 << endl; + } + + test = test || !(error2 < tol); + + x_global = 0; + + if (rank == 0) + std::cout << "RAS two level with overlap:" << std::endl; + MPI_Barrier(MPI_COMM_WORLD); + opt.parse("-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction additive"); + ddm_with_overlap.solve(f_global.data(), x_global.data(), mu); + ddm_with_overlap.print_infos(); + ddm_with_overlap.clean(); + error2 = normFrob(f_global - A * x_global) / normFrob(f_global); + if (rank == 0) { + cout << "error: " << error2 << endl; + } + + test = test || !(error2 < tol); + + x_global = 0; + + // DDM solver with threshold + if (rank == 0) + std::cout << "RAS two level with overlap and threshold:" << std::endl; + opt.parse("-hpddm_schwarz_method ras -hpddm_schwarz_coarse_correction additive -hpddm_geneo_threshold 100"); + + auto local_block_diagonal_hmatrix_bis = *default_build.block_diagonal_hmatrix; + DDMSolverBuilder> default_ddm_solver_with_threshold(Operator, local_block_diagonal_hmatrix_bis, Generator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections); + auto &ddm_with_overlap_threshold = default_ddm_solver_with_threshold.solver; + + Ki.bytes_to_matrix(datapath + "/Ki_" + NbrToStr(size) + "_" + NbrToStr(rank) + ".bin"); + int local_size_wo_overlap = Operator.get_target_partition().get_size_of_partition(rank); + int local_size_with_overlap = Ki.nb_cols(); + auto geneo_coarse_space_dense_builder = GeneoCoarseSpaceDenseBuilder>::GeneoWithThreshold(local_size_wo_overlap, local_size_with_overlap, *default_build.block_diagonal_hmatrix, Ki, symmetric, UPLO, 100); + // geneo_coarse_space_dense_builder.set_geneo_threshold(100); + GeneoCoarseOperatorBuilder> geneo_coarse_operator_builder(Operator); + ddm_with_overlap_threshold.build_coarse_space(geneo_coarse_space_dense_builder, geneo_coarse_operator_builder); + ddm_with_overlap_threshold.facto_one_level(); + ddm_with_overlap_threshold.solve(f_global.data(), x_global.data(), mu); + ddm_with_overlap_threshold.print_infos(); + ddm_with_overlap_threshold.clean(); + error2 = normFrob(f_global - A * x_global) / normFrob(f_global); + if (rank == 0) { + cout << "error: " << error2 << endl; + } + + test = test || !(error2 < tol); + + x_global = 0; + } return test; } diff --git a/tests/functional_tests/solvers/test_solver_wo_overlap.hpp b/tests/functional_tests/solvers/test_solver_wo_overlap.hpp index a42283a6..cb8fbd90 100644 --- a/tests/functional_tests/solvers/test_solver_wo_overlap.hpp +++ b/tests/functional_tests/solvers/test_solver_wo_overlap.hpp @@ -26,7 +26,7 @@ using namespace std; using namespace htool; template -int test_solver_wo_overlap(int argc, char *argv[], int mu, char symmetric, std::string datapath) { +int test_solver_wo_overlap(int argc, char *argv[], int mu, char symmetric, char UPLO, std::string datapath) { // Get the number of processes int size; @@ -51,11 +51,6 @@ int test_solver_wo_overlap(int argc, char *argv[], int mu, char symmetric, std:: double epsilon = tol; double eta = 10; - char UPLO = 'N'; - if (symmetric != 'N') { - UPLO = 'L'; - } - // Clustering if (rank == 0) std::cout << "Creating cluster tree" << std::endl;