From 8de449c7149951f3b7d2319709f129501b1257cf Mon Sep 17 00:00:00 2001 From: Pierre Marchand Date: Tue, 16 Jul 2024 17:00:12 +0200 Subject: [PATCH] add conversion in user numbering --- include/htool/hmatrix/hmatrix.hpp | 148 ++++++++++++++---- include/htool/misc/petsc.hpp | 7 +- .../test_hmatrix_build_complex_double.cpp | 40 +++-- .../test_hmatrix_build_double.cpp | 36 +++-- .../hmatrix/test_hmatrix_build.hpp | 51 ++++++ 5 files changed, 211 insertions(+), 71 deletions(-) diff --git a/include/htool/hmatrix/hmatrix.hpp b/include/htool/hmatrix/hmatrix.hpp index 74f1ace5..e842a24d 100644 --- a/include/htool/hmatrix/hmatrix.hpp +++ b/include/htool/hmatrix/hmatrix.hpp @@ -46,6 +46,7 @@ class HMatrix : public TreeNode m_leaves{}; mutable std::vector m_leaves_for_symmetry{}; mutable char m_symmetry_type_for_leaves{'N'}; + mutable char m_UPLO_for_leaves{'N'}; // std::vector m_leaves_in_diagonal_block{}; StorageType m_storage_type{StorageType::Hierarchical}; @@ -71,6 +72,7 @@ class HMatrix : public TreeNodeget_symmetry() != 'N') { m_symmetry_type_for_leaves = current_hmatrix->get_symmetry(); + m_UPLO_for_leaves = current_hmatrix->get_UPLO(); } for (const auto &child : current_hmatrix->get_children()) { @@ -251,6 +253,7 @@ class HMatrix : public TreeNode &get_infos() const { return infos; } @@ -1992,37 +1995,85 @@ void copy_to_dense(const HMatrix &hma } } } +} + +template +void copy_to_dense_in_user_numbering(const HMatrix &hmatrix, CoefficientPrecision *ptr) { + + const auto &target_cluster = hmatrix.get_target_cluster(); + const auto &source_cluster = hmatrix.get_source_cluster(); + if (!target_cluster.is_root() && !is_cluster_on_partition(target_cluster)) { + htool::Logger::get_instance().log(LogLevel::ERROR, "Target cluster is neither root nor local, permutation is not stable and copy_to_dense_in_user_numbering cannot be used."); // LCOV_EXCL_LINE + } + + if (!source_cluster.is_root() && !is_cluster_on_partition(source_cluster)) { + htool::Logger::get_instance().log(LogLevel::ERROR, "Source cluster is neither root nor local, permutation is not stable and copy_to_dense_in_user_numbering cannot be used."); // LCOV_EXCL_LINE + } + + if (is_cluster_on_partition(target_cluster) && !(target_cluster.is_permutation_local())) { + htool::Logger::get_instance().log(LogLevel::ERROR, "Target cluster is local, but permutation is not local, copy_to_dense_in_user_numbering cannot be used"); // LCOV_EXCL_LINE} + } + + if (is_cluster_on_partition(source_cluster) && !(source_cluster.is_permutation_local())) { + htool::Logger::get_instance().log(LogLevel::ERROR, "Source cluster is local, but permutation is not local, copy_to_dense_in_user_numbering cannot be used"); // LCOV_EXCL_LINE} + } + int target_offset = target_cluster.get_offset(); + int source_offset = source_cluster.get_offset(); + int target_size = target_cluster.get_size(); + const auto &target_permutation = target_cluster.get_permutation(); + const auto &source_permutation = source_cluster.get_permutation(); + + for (auto leaf : hmatrix.get_leaves()) { + int local_nr = leaf->get_target_cluster().get_size(); + int local_nc = leaf->get_source_cluster().get_size(); + int row_offset = leaf->get_target_cluster().get_offset(); + int col_offset = leaf->get_source_cluster().get_offset(); + + if (leaf->is_dense()) { + for (int k = 0; k < local_nc; k++) { + for (int j = 0; j < local_nr; j++) { + ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size] = (*leaf->get_dense_data())(j, k); + } + } + } else { - // char symmetry_type = hmatrix.get_symmetry_for_leaves(); - // for (auto leaf : hmatrix.get_leaves_for_symmetry()) { - // int local_nr = leaf->get_target_cluster().get_size(); - // int local_nc = leaf->get_source_cluster().get_size(); - // int col_offset = leaf->get_target_cluster().get_offset() - source_offset; - // int row_offset = leaf->get_source_cluster().get_offset() - target_offset; - // if (leaf->is_dense()) { - // for (int j = 0; j < local_nr; j++) { - // for (int k = 0; k < local_nc; k++) { - // ptr[k + row_offset + (j + col_offset) * target_size] = (*leaf->get_dense_data())(j, k); - // } - // } - // } else { - - // Matrix low_rank_to_dense(local_nr, local_nc); - // leaf->get_low_rank_data()->copy_to_dense(low_rank_to_dense.data()); - // for (int j = 0; j < local_nr; j++) { - // for (int k = 0; k < local_nc; k++) { - // ptr[k + row_offset + (j + col_offset) * target_size] = low_rank_to_dense(j, k); - // } - // } - // } - // if (symmetry_type == 'H') { - // for (int k = 0; k < local_nc; k++) { - // for (int j = 0; j < local_nr; j++) { - // ptr[k + row_offset + (j + col_offset) * target_size] = conj_if_complex(ptr[k + row_offset + (j + col_offset) * target_size]); - // } - // } - // } - // } + Matrix low_rank_to_dense(local_nr, local_nc); + leaf->get_low_rank_data()->copy_to_dense(low_rank_to_dense.data()); + for (int k = 0; k < local_nc; k++) { + for (int j = 0; j < local_nr; j++) { + ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size] = low_rank_to_dense(j, k); + } + } + } + } + if (hmatrix.get_symmetry_for_leaves() != 'N') { + for (auto leaf : hmatrix.get_leaves_for_symmetry()) { + int local_nr = leaf->get_target_cluster().get_size(); + int local_nc = leaf->get_source_cluster().get_size(); + int row_offset = leaf->get_target_cluster().get_offset(); + int col_offset = leaf->get_source_cluster().get_offset(); + + if (leaf->get_target_cluster().get_offset() == leaf->get_source_cluster().get_offset() and hmatrix.get_UPLO_for_leaves() == 'L') { + for (int k = 0; k < local_nc; k++) { + for (int j = k + 1; j < local_nr; j++) { + ptr[source_permutation[k + col_offset] - target_offset + (target_permutation[j + row_offset] - source_offset) * target_size] = hmatrix.get_symmetry_for_leaves() == 'S' ? ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size] : conj_if_complex(ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size]); + } + } + } else if (leaf->get_target_cluster().get_offset() == leaf->get_source_cluster().get_offset() and hmatrix.get_UPLO_for_leaves() == 'U') { + for (int k = 0; k < local_nc; k++) { + for (int j = 0; j < k; j++) { + ptr[source_permutation[k + col_offset] - target_offset + (target_permutation[j + row_offset] - source_offset) * target_size] = hmatrix.get_symmetry_for_leaves() == 'S' ? ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size] : conj_if_complex(ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size]); + } + } + } else { + for (int k = 0; k < local_nc; k++) { + for (int j = 0; j < local_nr; j++) { + ptr[source_permutation[k + col_offset] - target_offset + (target_permutation[j + row_offset] - source_offset) * target_size] = hmatrix.get_symmetry_for_leaves() == 'S' ? ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size] : conj_if_complex(ptr[target_permutation[j + row_offset] - target_offset + (source_permutation[k + col_offset] - source_offset) * target_size]); + } + } + } + } + } } template @@ -2054,6 +2105,43 @@ void copy_diagonal(const HMatrix &hma } } +template +void copy_diagonal_in_user_numbering(const HMatrix &hmatrix, CoefficientPrecision *ptr) { + if (hmatrix.get_target_cluster().get_offset() != hmatrix.get_source_cluster().get_offset() || hmatrix.get_target_cluster().get_size() != hmatrix.get_source_cluster().get_size()) { + htool::Logger::get_instance().log(LogLevel::ERROR, "Matrix is not square a priori, get_local_diagonal cannot be used"); // LCOV_EXCL_LINE + } + + if (!(hmatrix.get_target_cluster().is_root() && hmatrix.get_source_cluster().is_root()) && !(is_cluster_on_partition(hmatrix.get_target_cluster()) && is_cluster_on_partition(hmatrix.get_source_cluster()))) { + htool::Logger::get_instance().log(LogLevel::ERROR, "Clusters are neither root nor local, permutations are not stable and copy_diagonal_in_user_numbering cannot be used."); // LCOV_EXCL_LINE + } + + if ((is_cluster_on_partition(hmatrix.get_target_cluster()) && is_cluster_on_partition(hmatrix.get_source_cluster())) && !(hmatrix.get_target_cluster().is_permutation_local() && hmatrix.get_source_cluster().is_permutation_local())) { + htool::Logger::get_instance().log(LogLevel::ERROR, "Clusters are local, but permutations are not local, copy_diagonal_in_user_numbering cannot be used"); // LCOV_EXCL_LINE} + } + + const auto &permutation = hmatrix.get_target_cluster().get_permutation(); + int target_offset = hmatrix.get_target_cluster().get_offset(); + // int source_offset = hmatrix.get_source_cluster().get_offset(); + + for (auto leaf : hmatrix.get_leaves()) { + int local_nr = leaf->get_target_cluster().get_size(); + int local_nc = leaf->get_source_cluster().get_size(); + int offset_i = leaf->get_target_cluster().get_offset(); + int offset_j = leaf->get_source_cluster().get_offset(); + if (leaf->is_dense() && offset_i == offset_j) { + for (int k = 0; k < std::min(local_nr, local_nc); k++) { + ptr[permutation[k + offset_i] - target_offset] = (*leaf->get_dense_data())(k, k); + } + } else if (leaf->is_low_rank() && offset_i == offset_j) { // pretty rare... + Matrix low_rank_to_dense(local_nc, local_nr); + leaf->get_low_rank_data()->copy_to_dense(low_rank_to_dense.data()); + for (int k = 0; k < std::min(local_nr, local_nc); k++) { + ptr[permutation[k + offset_i] - target_offset] = (low_rank_to_dense)(k, k); + } + } + } +} + // template // Matrix HMatrix::get_local_interaction(bool permutation) const { // int local_size_source = cluster_tree_s->get_masteroffset(rankWorld).second; diff --git a/include/htool/misc/petsc.hpp b/include/htool/misc/petsc.hpp index d29fd5a7..6f629ee6 100644 --- a/include/htool/misc/petsc.hpp +++ b/include/htool/misc/petsc.hpp @@ -12,13 +12,14 @@ #include "define.hpp" #include "misc.hpp" -#include "../clustering/bounding_box_1.hpp" -#include "../clustering/pca.hpp" +#include "../clustering/clustering.hpp" +#include "../distributed_operator/distributed_operator.hpp" +#include "../distributed_operator/implementations/partition_from_cluster.hpp" +#include "../hmatrix/hmatrix.hpp" #include "../hmatrix/lrmat/SVD.hpp" #include "../hmatrix/lrmat/fullACA.hpp" #include "../hmatrix/lrmat/lrmat.hpp" #include "../hmatrix/lrmat/sympartialACA.hpp" -#include "../types/hmatrix.hpp" #if defined(__clang__) # pragma clang diagnostic pop diff --git a/tests/functional_tests/hmatrix/hmatrix_build/test_hmatrix_build_complex_double.cpp b/tests/functional_tests/hmatrix/hmatrix_build/test_hmatrix_build_complex_double.cpp index f0b18a75..67afadeb 100644 --- a/tests/functional_tests/hmatrix/hmatrix_build/test_hmatrix_build_complex_double.cpp +++ b/tests/functional_tests/hmatrix/hmatrix_build/test_hmatrix_build_complex_double.cpp @@ -6,27 +6,25 @@ using namespace htool; int main(int argc, char *argv[]) { MPI_Init(&argc, &argv); - bool is_error = false; - const int number_of_rows = 200; - const int number_of_rows_increased = 400; - const int number_of_columns = 200; - const int number_of_columns_increased = 400; - - for (auto use_local_cluster : {true, false}) { - for (auto epsilon : {1e-14, 1e-6}) { - for (auto use_dense_block_generator : {true, false}) { - // Square matrix - is_error = is_error || test_hmatrix_build, GeneratorTestComplexSymmetric>(number_of_rows, number_of_columns, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator); - - is_error = is_error || test_hmatrix_build, GeneratorTestComplexSymmetric>(number_of_rows, number_of_columns, use_local_cluster, 'S', 'U', epsilon, use_dense_block_generator); - is_error = is_error || test_hmatrix_build, GeneratorTestComplexSymmetric>(number_of_rows, number_of_columns, use_local_cluster, 'S', 'L', epsilon, use_dense_block_generator); - - is_error = is_error || test_hmatrix_build, GeneratorTestComplexHermitian>(number_of_rows, number_of_columns, use_local_cluster, 'H', 'U', epsilon, use_dense_block_generator); - is_error = is_error || test_hmatrix_build, GeneratorTestComplexHermitian>(number_of_rows, number_of_columns, use_local_cluster, 'H', 'L', epsilon, use_dense_block_generator); - - // Rectangle matrix - is_error = is_error || test_hmatrix_build, GeneratorTestComplex>(number_of_rows_increased, number_of_columns, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator); - is_error = is_error || test_hmatrix_build, GeneratorTestComplex>(number_of_rows, number_of_columns_increased, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator); + bool is_error = false; + + for (auto nr : {200, 400}) { + for (auto nc : {200, 400}) { + for (auto use_local_cluster : {true, false}) { + for (auto epsilon : {1e-14, 1e-6}) { + for (auto use_dense_block_generator : {true, false}) { + std::cout << nr << " " << nc << " " << use_local_cluster << " " << epsilon << " " << use_dense_block_generator << "\n"; + + is_error = is_error || test_hmatrix_build, GeneratorTestComplexSymmetric>(nr, nc, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator); + if (nr == nc) { + for (auto UPLO : {'U', 'L'}) { + is_error = is_error || test_hmatrix_build, GeneratorTestComplexSymmetric>(nr, nr, use_local_cluster, 'S', UPLO, epsilon, use_dense_block_generator); + + is_error = is_error || test_hmatrix_build, GeneratorTestComplexHermitian>(nr, nr, use_local_cluster, 'H', UPLO, epsilon, use_dense_block_generator); + } + } + } + } } } } diff --git a/tests/functional_tests/hmatrix/hmatrix_build/test_hmatrix_build_double.cpp b/tests/functional_tests/hmatrix/hmatrix_build/test_hmatrix_build_double.cpp index 3aea02fc..dd2c6cdb 100644 --- a/tests/functional_tests/hmatrix/hmatrix_build/test_hmatrix_build_double.cpp +++ b/tests/functional_tests/hmatrix/hmatrix_build/test_hmatrix_build_double.cpp @@ -6,23 +6,25 @@ using namespace htool; int main(int argc, char *argv[]) { MPI_Init(&argc, &argv); - bool is_error = false; - const int number_of_rows = 200; - const int number_of_rows_increased = 400; - const int number_of_columns = 200; - const int number_of_columns_increased = 400; - - for (auto use_local_cluster : {true, false}) { - for (auto epsilon : {1e-14, 1e-6}) { - for (auto use_dense_block_generator : {true, false}) { - // Square matrix - is_error = is_error || test_hmatrix_build(number_of_rows, number_of_columns, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator); - is_error = is_error || test_hmatrix_build(number_of_rows, number_of_columns, use_local_cluster, 'S', 'U', epsilon, use_dense_block_generator); - is_error = is_error || test_hmatrix_build(number_of_rows, number_of_columns, use_local_cluster, 'S', 'L', epsilon, use_dense_block_generator); - - // Rectangle matrix - is_error = is_error || test_hmatrix_build(number_of_rows_increased, number_of_columns, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator); - is_error = is_error || test_hmatrix_build(number_of_rows, number_of_columns_increased, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator); + bool is_error = false; + + for (auto nr : {200, 400}) { + for (auto nc : {200, 400}) { + for (auto use_local_cluster : {true, false}) { + for (auto epsilon : {1e-14, 1e-6}) { + for (auto use_dense_block_generator : {true, false}) { + std::cout << nr << " " << nc << " " << use_local_cluster << " " << epsilon << " " << use_dense_block_generator << "\n"; + + is_error = is_error || test_hmatrix_build(nr, nc, use_local_cluster, 'N', 'N', epsilon, use_dense_block_generator); + + if (nr == nc) { + for (auto UPLO : {'U', 'L'}) { + std::cout << UPLO << "\n"; + is_error = is_error || test_hmatrix_build(nr, nc, use_local_cluster, 'S', UPLO, epsilon, use_dense_block_generator); + } + } + } + } } } } diff --git a/tests/functional_tests/hmatrix/test_hmatrix_build.hpp b/tests/functional_tests/hmatrix/test_hmatrix_build.hpp index 7a26d967..324cad80 100644 --- a/tests/functional_tests/hmatrix/test_hmatrix_build.hpp +++ b/tests/functional_tests/hmatrix/test_hmatrix_build.hpp @@ -193,6 +193,49 @@ bool test_hmatrix_build(int nr, int nc, bool use_local_cluster, char Symmetry, c is_error = is_error || !(frobenius_error < epsilon); std::cout << "Check dense conversion: " << frobenius_error << "\n"; + // Check dense conversion with permutation + Matrix permuted_matrix(hmatrix_target_cluster.get_size(), hmatrix_source_cluster.get_size()), temp_matrix(hmatrix_target_cluster.get_size(), hmatrix_source_cluster.get_size()); + copy_to_dense_in_user_numbering(root_hmatrix, permuted_matrix.data()); + + const auto &source_permutation = hmatrix_source_cluster.get_permutation(); + for (int i = 0; i < hmatrix_target_cluster.get_size(); i++) { + for (int j = 0; j < hmatrix_source_cluster.get_size(); j++) { + temp_matrix(i, j) = permuted_matrix(target_permutation[i + hmatrix_target_cluster.get_offset()] - hmatrix_target_cluster.get_offset(), source_permutation[j + hmatrix_source_cluster.get_offset()] - hmatrix_source_cluster.get_offset()); + } + } + // MPI_Barrier(MPI_COMM_WORLD); + // if (rankWorld == 0) { + // dense_matrix.print(std::cout, ","); + // std::cout << root_hmatrix.get_symmetry() << "\n"; + // } + // MPI_Barrier(MPI_COMM_WORLD); + // MPI_Barrier(MPI_COMM_WORLD); + // if (rankWorld == 1) { + // dense_matrix.print(std::cout, ","); + // std::cout << root_hmatrix.get_symmetry() << "\n"; + // } + // MPI_Barrier(MPI_COMM_WORLD); + // MPI_Barrier(MPI_COMM_WORLD); + // if (rankWorld == 0) { + // save_leaves_with_rank(root_hmatrix, "test_0"); + // Matrix test = temp_matrix - dense_matrix; + // test.print(std::cout, ","); + // std::cout << root_hmatrix.get_symmetry() << "\n"; + // } + // MPI_Barrier(MPI_COMM_WORLD); + // MPI_Barrier(MPI_COMM_WORLD); + // if (rankWorld == 1) { + // save_leaves_with_rank(root_hmatrix, "test_1"); + // Matrix test = temp_matrix - dense_matrix; + // test.print(std::cout, ","); + // std::cout << root_hmatrix.get_symmetry() << "\n"; + // } + // MPI_Barrier(MPI_COMM_WORLD); + frobenius_error = normFrob(temp_matrix - dense_matrix); + frobenius_error /= dense_matrix_norm; + is_error = is_error || !(frobenius_error < epsilon); + std::cout << "Check dense conversion in user numbering: " << frobenius_error << "\n"; + // Check get diagonal conversion if (Symmetry != 'N' || nr == nc) { int local_size = diagonal_hmatrix->get_target_cluster().get_size(); @@ -208,6 +251,14 @@ bool test_hmatrix_build(int nr, int nc, bool use_local_cluster, char Symmetry, c is_error = is_error || !(error_on_diagonal < epsilon); std::cout << "Check get diagonal: " << error_on_diagonal << "\n"; + + vector diagonal_in_user_numbering(local_size), temp_vector(local_size); + copy_diagonal_in_user_numbering(*diagonal_hmatrix, diagonal_in_user_numbering.data()); + user_to_cluster(diagonal_hmatrix->get_target_cluster(), diagonal_in_user_numbering.data(), temp_vector.data()); + error_on_diagonal = norm2(temp_vector - dense_diagonal) / norm2(dense_diagonal); + + is_error = is_error || !(error_on_diagonal < epsilon); + std::cout << "Check get diagonal in user numbering: " << error_on_diagonal << "\n"; } //