From d3563fcb8e595760e3c02a9db1b39883d26c47f6 Mon Sep 17 00:00:00 2001 From: Pierre Marchand Date: Thu, 18 Jul 2024 17:36:22 +0200 Subject: [PATCH] change in generator with permutation + clean up --- .../htool/distributed_operator/utility.hpp | 7 + include/htool/hmatrix/hmatrix.hpp | 4 - .../hmatrix/interfaces/virtual_generator.hpp | 24 +- .../hmatrix/tree_builder/tree_builder.hpp | 3 + .../htool/local_operators/local_operator.hpp | 3 - include/htool/testing/generate_test_case.hpp | 87 +++-- include/htool/testing/generator_test.hpp | 313 +++++++++++------- .../test_distributed_operator.hpp | 199 ++--------- ...ibuted_operator_product_complex_double.cpp | 4 +- ...st_distributed_operator_product_double.cpp | 4 +- .../test_hmatrix_addition_complex_double.cpp | 4 +- .../test_hmatrix_addition_double.cpp | 4 +- ...t_hmatrix_factorization_complex_double.cpp | 5 +- .../test_hmatrix_factorization_double.cpp | 5 +- .../test_hmatrix_product_complex_double.cpp | 1 + ...matrix_triangular_solve_complex_double.cpp | 4 +- .../test_hmatrix_triangular_solve_double.cpp | 4 +- .../lrmat_build/test_lrmat_build_SVD.cpp | 3 +- .../lrmat_build/test_lrmat_build_fullACA.cpp | 3 +- .../test_lrmat_build_partialACA.cpp | 3 +- .../test_lrmat_build_sympartialACA.cpp | 3 +- .../hmatrix/lrmat/test_lrmat_build.hpp | 2 +- .../hmatrix/test_hmatrix_build.hpp | 66 +--- .../hmatrix/test_hmatrix_factorization.hpp | 2 +- .../hmatrix/test_hmatrix_hmatrix_product.hpp | 1 + 25 files changed, 333 insertions(+), 425 deletions(-) diff --git a/include/htool/distributed_operator/utility.hpp b/include/htool/distributed_operator/utility.hpp index 6665c7de..81802af3 100644 --- a/include/htool/distributed_operator/utility.hpp +++ b/include/htool/distributed_operator/utility.hpp @@ -30,6 +30,9 @@ class DistributedOperatorFromHMatrix { distributed_operator.add_local_operator(&local_hmatrix); block_diagonal_hmatrix = hmatrix.get_sub_hmatrix(hmatrix_builder.get_target_cluster().get_cluster_on_partition(get_rankWorld(communicator)), hmatrix_builder.get_source_cluster().get_cluster_on_partition(get_rankWorld(communicator))); } + + DistributedOperatorFromHMatrix(const VirtualGeneratorInUserNumbering &generator, const Cluster &target_cluster, const Cluster &source_cluster, const HMatrixTreeBuilder &hmatrix_builder, MPI_Comm communicator) : DistributedOperatorFromHMatrix(GeneratorWithPermutation(generator, target_cluster.get_permutation().data(), source_cluster.get_permutation().data()), target_cluster, source_cluster, hmatrix_builder, communicator) { + } }; template > @@ -49,6 +52,8 @@ class DefaultApproximationBuilder { const HMatrix *block_diagonal_hmatrix{nullptr}; DefaultApproximationBuilder(const VirtualGenerator &generator, const Cluster &target_cluster, const Cluster &source_cluster, htool::underlying_type epsilon, htool::underlying_type eta, char symmetry, char UPLO, MPI_Comm communicator) : distributed_operator_builder(generator, target_cluster, source_cluster, HMatrixTreeBuilder(target_cluster, source_cluster, epsilon, eta, symmetry, UPLO, -1, get_rankWorld(communicator), get_rankWorld(communicator)), communicator), hmatrix(distributed_operator_builder.hmatrix), distributed_operator(distributed_operator_builder.distributed_operator), block_diagonal_hmatrix(distributed_operator_builder.block_diagonal_hmatrix) {} + + DefaultApproximationBuilder(const VirtualGeneratorInUserNumbering &generator, const Cluster &target_cluster, const Cluster &source_cluster, htool::underlying_type epsilon, htool::underlying_type eta, char symmetry, char UPLO, MPI_Comm communicator) : DefaultApproximationBuilder(GeneratorWithPermutation(generator, target_cluster.get_permutation().data(), source_cluster.get_permutation().data()), target_cluster, source_cluster, epsilon, eta, symmetry, UPLO, communicator) {} }; template > @@ -69,6 +74,8 @@ class DefaultLocalApproximationBuilder { public: DefaultLocalApproximationBuilder(const VirtualGenerator &generator, const Cluster &target_cluster, const Cluster &source_cluster, htool::underlying_type epsilon, htool::underlying_type eta, char symmetry, char UPLO, MPI_Comm communicator) : distributed_operator_builder(generator, target_cluster, source_cluster, HMatrixTreeBuilder(target_cluster.get_cluster_on_partition(get_rankWorld(communicator)), source_cluster.get_cluster_on_partition(get_rankWorld(communicator)), epsilon, eta, symmetry, UPLO, -1, get_rankWorld(communicator), get_rankWorld(communicator)), communicator), hmatrix(distributed_operator_builder.hmatrix), distributed_operator(distributed_operator_builder.distributed_operator), block_diagonal_hmatrix(distributed_operator_builder.block_diagonal_hmatrix) {} + + DefaultLocalApproximationBuilder(const VirtualGeneratorInUserNumbering &generator, const Cluster &target_cluster, const Cluster &source_cluster, htool::underlying_type epsilon, htool::underlying_type eta, char symmetry, char UPLO, MPI_Comm communicator) : DefaultLocalApproximationBuilder(GeneratorWithPermutation(generator, target_cluster.get_permutation().data(), source_cluster.get_permutation().data()), target_cluster, source_cluster, epsilon, eta, symmetry, UPLO, communicator) {} }; } // namespace htool diff --git a/include/htool/hmatrix/hmatrix.hpp b/include/htool/hmatrix/hmatrix.hpp index e842a24d..ea9d7cf7 100644 --- a/include/htool/hmatrix/hmatrix.hpp +++ b/include/htool/hmatrix/hmatrix.hpp @@ -14,7 +14,6 @@ #include "lrmat/lrmat.hpp" #include -#include namespace htool { // Class @@ -989,9 +988,6 @@ void HMatrix::threaded_hierarchical_a set_leaves_in_cache(); - int rankWorld; - MPI_Comm_rank(MPI_COMM_WORLD, &rankWorld); - int out_size(m_target_cluster->get_size()); auto get_output_cluster{&HMatrix::get_target_cluster}; auto get_input_cluster{&HMatrix::get_source_cluster}; diff --git a/include/htool/hmatrix/interfaces/virtual_generator.hpp b/include/htool/hmatrix/interfaces/virtual_generator.hpp index 5803c900..af5c1668 100644 --- a/include/htool/hmatrix/interfaces/virtual_generator.hpp +++ b/include/htool/hmatrix/interfaces/virtual_generator.hpp @@ -11,7 +11,6 @@ class VirtualGenerator { public: virtual void copy_submatrix(int M, int N, int row_offset, int col_offset, CoefficientPrecision *ptr) const = 0; - // virtual void copy_submatrix(int M, int N, const int *rows, const int *cols, CoefficientPrecision *ptr) const = 0; VirtualGenerator() {} VirtualGenerator(const VirtualGenerator &) = default; @@ -22,21 +21,34 @@ class VirtualGenerator { }; template -class VirtualGeneratorWithPermutation : public VirtualGenerator { +class VirtualGeneratorInUserNumbering { + + public: + virtual void copy_submatrix(int M, int N, const int *rows, const int *cols, CoefficientPrecision *ptr) const = 0; + + VirtualGeneratorInUserNumbering() {} + VirtualGeneratorInUserNumbering(const VirtualGeneratorInUserNumbering &) = default; + VirtualGeneratorInUserNumbering &operator=(const VirtualGeneratorInUserNumbering &) = default; + VirtualGeneratorInUserNumbering(VirtualGeneratorInUserNumbering &&) = default; + VirtualGeneratorInUserNumbering &operator=(VirtualGeneratorInUserNumbering &&) = default; + virtual ~VirtualGeneratorInUserNumbering() {} +}; + +template +class GeneratorWithPermutation : public VirtualGenerator { protected: + const VirtualGeneratorInUserNumbering &m_generator_in_user_numbering; const int *m_target_permutation; const int *m_source_permutation; public: - VirtualGeneratorWithPermutation(const int *target_permutation, const int *source_permutation) : m_target_permutation(target_permutation), m_source_permutation(source_permutation) { + GeneratorWithPermutation(const VirtualGeneratorInUserNumbering &generator_in_user_numbering, const int *target_permutation, const int *source_permutation) : m_generator_in_user_numbering(generator_in_user_numbering), m_target_permutation(target_permutation), m_source_permutation(source_permutation) { } virtual void copy_submatrix(int M, int N, int row_offset, int col_offset, CoefficientPrecision *ptr) const override { - copy_submatrix_from_user_numbering(M, N, m_target_permutation + row_offset, m_source_permutation + col_offset, ptr); + m_generator_in_user_numbering.copy_submatrix(M, N, m_target_permutation + row_offset, m_source_permutation + col_offset, ptr); } - - virtual void copy_submatrix_from_user_numbering(int M, int N, const int *rows, const int *cols, CoefficientPrecision *ptr) const = 0; }; } // namespace htool diff --git a/include/htool/hmatrix/tree_builder/tree_builder.hpp b/include/htool/hmatrix/tree_builder/tree_builder.hpp index a2be060a..5afc3402 100644 --- a/include/htool/hmatrix/tree_builder/tree_builder.hpp +++ b/include/htool/hmatrix/tree_builder/tree_builder.hpp @@ -123,6 +123,9 @@ class HMatrixTreeBuilder { // Build HMatrixType build(const VirtualGenerator &generator) const; + HMatrixType build(const VirtualGeneratorInUserNumbering &generator) const { + return this->build(GeneratorWithPermutation(generator, m_target_root_cluster.get_permutation().data(), m_source_root_cluster.get_permutation().data())); + } // Setters void set_low_rank_generator(std::shared_ptr> ptr) { m_low_rank_generator = ptr; } diff --git a/include/htool/local_operators/local_operator.hpp b/include/htool/local_operators/local_operator.hpp index 4aa149ec..051b5ef3 100644 --- a/include/htool/local_operators/local_operator.hpp +++ b/include/htool/local_operators/local_operator.hpp @@ -35,9 +35,6 @@ class LocalOperator : public VirtualLocalOperator { virtual void local_sub_matrix_vector_product(const CoefficientPrecision *const in, CoefficientPrecision *const out, int mu, int offset, int size) const { std::vector temp(m_source_cluster.get_size() * mu, 0); - int rankWorld; - MPI_Comm_rank(MPI_COMM_WORLD, &rankWorld); - std::copy_n(in, size * mu, temp.data() + offset * mu); add_matrix_product_global_to_local(1, temp.data(), 0, out, mu); } diff --git a/include/htool/testing/generate_test_case.hpp b/include/htool/testing/generate_test_case.hpp index e15b5401..a44ab0a9 100644 --- a/include/htool/testing/generate_test_case.hpp +++ b/include/htool/testing/generate_test_case.hpp @@ -9,6 +9,10 @@ namespace htool { template class TestCaseProduct { + std::unique_ptr operator_in_user_numbering_A = nullptr; + std::unique_ptr operator_in_user_numbering_B = nullptr; + std::unique_ptr operator_in_user_numbering_C = nullptr; + public: char transa; char transb; @@ -24,9 +28,9 @@ class TestCaseProduct { const Cluster> *root_cluster_B_output = nullptr; const Cluster> *root_cluster_C_input = nullptr; const Cluster> *root_cluster_C_output = nullptr; - std::unique_ptr operator_A = nullptr; - std::unique_ptr operator_B = nullptr; - std::unique_ptr operator_C = nullptr; + std::unique_ptr> operator_A = nullptr; + std::unique_ptr> operator_B = nullptr; + std::unique_ptr> operator_C = nullptr; int ni_A; int no_A; int ni_B; @@ -79,26 +83,30 @@ class TestCaseProduct { // Operators if (transa == 'N') { - operator_A = std::make_unique(3, no_A, ni_A, x1, x2, *root_cluster_1, *root_cluster_2, true, true); - root_cluster_A_output = root_cluster_1.get(); - root_cluster_A_input = root_cluster_2.get(); + operator_in_user_numbering_A = std::make_unique(3, x1, x2); + root_cluster_A_output = root_cluster_1.get(); + root_cluster_A_input = root_cluster_2.get(); } else { - operator_A = std::make_unique(3, no_A, ni_A, x2, x1, *root_cluster_2, *root_cluster_1, true, true); - root_cluster_A_output = root_cluster_2.get(); - root_cluster_A_input = root_cluster_1.get(); + operator_in_user_numbering_A = std::make_unique(3, x2, x1); + root_cluster_A_output = root_cluster_2.get(); + root_cluster_A_input = root_cluster_1.get(); } if (transb == 'N') { - operator_B = std::make_unique(3, no_B, ni_B, x2, x3, *root_cluster_2, *root_cluster_3, true, true); - root_cluster_B_output = root_cluster_2.get(); - root_cluster_B_input = root_cluster_3.get(); + operator_in_user_numbering_B = std::make_unique(3, x2, x3); + root_cluster_B_output = root_cluster_2.get(); + root_cluster_B_input = root_cluster_3.get(); } else { - operator_B = std::make_unique(3, no_B, ni_B, x3, x2, *root_cluster_3, *root_cluster_2, true, true); - root_cluster_B_output = root_cluster_3.get(); - root_cluster_B_input = root_cluster_2.get(); + operator_in_user_numbering_B = std::make_unique(3, x3, x2); + root_cluster_B_output = root_cluster_3.get(); + root_cluster_B_input = root_cluster_2.get(); } - operator_C = std::make_unique(3, no_C, ni_C, x1, x3, *root_cluster_1, *root_cluster_3, true, true); - root_cluster_C_input = root_cluster_3.get(); - root_cluster_C_output = root_cluster_1.get(); + operator_in_user_numbering_C = std::make_unique(3, x1, x3); + root_cluster_C_input = root_cluster_3.get(); + root_cluster_C_output = root_cluster_1.get(); + + operator_A = std::make_unique>(*operator_in_user_numbering_A, root_cluster_A_output->get_permutation().data(), root_cluster_A_input->get_permutation().data()); + operator_B = std::make_unique>(*operator_in_user_numbering_B, root_cluster_B_output->get_permutation().data(), root_cluster_B_input->get_permutation().data()); + operator_C = std::make_unique>(*operator_in_user_numbering_C, root_cluster_C_output->get_permutation().data(), root_cluster_C_input->get_permutation().data()); } }; @@ -257,6 +265,10 @@ class TestCaseSymmetricRankUpdate { template class TestCaseSolve { + private: + std::unique_ptr operator_in_user_numbering_A = nullptr; + std::unique_ptr operator_in_user_numbering_X = nullptr; + public: char side; char trans; @@ -268,8 +280,8 @@ class TestCaseSolve { const Cluster> *root_cluster_A_output = nullptr; const Cluster> *root_cluster_X_input = nullptr; const Cluster> *root_cluster_X_output = nullptr; - std::unique_ptr operator_A = nullptr; - std::unique_ptr operator_X = nullptr; + std::unique_ptr> operator_A = nullptr; + std::unique_ptr> operator_X = nullptr; int ni_A; int no_A; int ni_X; @@ -306,24 +318,30 @@ class TestCaseSolve { } // Operators - operator_A = std::make_unique(3, no_A, ni_A, x1, x1, *root_cluster_1, *root_cluster_1, true, true); - root_cluster_A_output = root_cluster_1.get(); - root_cluster_A_input = root_cluster_1.get(); + operator_in_user_numbering_A = std::make_unique(3, x1, x1); + root_cluster_A_output = root_cluster_1.get(); + root_cluster_A_input = root_cluster_1.get(); if (side == 'L') { - operator_X = std::make_unique(3, no_X, ni_X, x1, x2, *root_cluster_1, *root_cluster_2, true, true); - root_cluster_X_output = root_cluster_1.get(); - root_cluster_X_input = root_cluster_2.get(); + operator_in_user_numbering_X = std::make_unique(3, x1, x2); + root_cluster_X_output = root_cluster_1.get(); + root_cluster_X_input = root_cluster_2.get(); } else { - operator_X = std::make_unique(3, no_X, ni_X, x2, x1, *root_cluster_2, *root_cluster_1, true, true); - root_cluster_X_output = root_cluster_2.get(); - root_cluster_X_input = root_cluster_1.get(); + operator_in_user_numbering_X = std::make_unique(3, x2, x1); + root_cluster_X_output = root_cluster_2.get(); + root_cluster_X_input = root_cluster_1.get(); } + + operator_A = std::make_unique>(*operator_in_user_numbering_A, root_cluster_A_output->get_permutation().data(), root_cluster_A_input->get_permutation().data()); + operator_X = std::make_unique>(*operator_in_user_numbering_X, root_cluster_X_output->get_permutation().data(), root_cluster_X_input->get_permutation().data()); } }; template class TestCaseAddition { + std::unique_ptr operator_in_user_numbering_A = nullptr; + std::unique_ptr operator_in_user_numbering_B = nullptr; + public: std::vector> x1; std::vector> x2; @@ -334,8 +352,8 @@ class TestCaseAddition { const Cluster> *root_cluster_A_output = nullptr; const Cluster> *root_cluster_B_input = nullptr; const Cluster> *root_cluster_B_output = nullptr; - std::unique_ptr operator_A = nullptr; - std::unique_ptr operator_B = nullptr; + std::unique_ptr> operator_A = nullptr; + std::unique_ptr> operator_B = nullptr; int ni_A; int no_A; int ni_B; @@ -371,9 +389,10 @@ class TestCaseAddition { } // Operators - operator_A = std::make_unique(3, no_A, ni_A, x1, x2, *root_cluster_1, *root_cluster_2, true, true); - root_cluster_A_output = root_cluster_1.get(); - root_cluster_A_input = root_cluster_2.get(); + operator_in_user_numbering_A = std::make_unique(3, x1, x2); + root_cluster_A_output = root_cluster_1.get(); + root_cluster_A_input = root_cluster_2.get(); + operator_A = std::make_unique>(*operator_in_user_numbering_A, root_cluster_A_output->get_permutation().data(), root_cluster_A_input->get_permutation().data()); // Sub lrmat two level deep std::random_device rd; diff --git a/include/htool/testing/generator_test.hpp b/include/htool/testing/generator_test.hpp index cab6d0bb..bb6e9dfd 100644 --- a/include/htool/testing/generator_test.hpp +++ b/include/htool/testing/generator_test.hpp @@ -10,138 +10,158 @@ namespace htool { template > -class GeneratorTest : public VirtualGenerator { +class GeneratorTestWithPermutation : public VirtualGeneratorInUserNumbering { protected: - int m_target_size, m_source_size; - const std::vector &p1; - const std::vector &p2; - const Cluster &m_target_cluster; - const Cluster &m_source_cluster; - bool m_use_target_permutation{true}; - bool m_use_source_permutation{true}; - int space_dim; + int m_space_dimension; + const std::vector &m_target_points; + const std::vector &m_source_points; public: - explicit GeneratorTest(int space_dim0, int nr0, int nc0, const std::vector> &p10, const std::vector> &p20, const Cluster &target_cluster, const Cluster &source_cluster, bool use_target_permutation, bool use_source_permutation) : m_target_size(nr0), m_source_size(nc0), p1(p10), p2(p20), m_target_cluster(target_cluster), m_source_cluster(source_cluster), m_use_target_permutation(use_target_permutation), m_use_source_permutation(use_source_permutation), space_dim(space_dim0) {} + GeneratorTestWithPermutation(int space_dim, const std::vector &target_points, const std::vector &source_points) : m_space_dimension(space_dim), m_target_points(target_points), m_source_points(source_points) {} virtual CoefficientPrecision get_coef(const int &i, const int &j) const = 0; - - void copy_submatrix(int M, int N, int row_offset, int col_offset, CoefficientPrecision *ptr) const override { - if (m_use_target_permutation && m_use_source_permutation) { - const auto &target_permutation = m_target_cluster.get_permutation(); - const auto &source_permutation = m_source_cluster.get_permutation(); - for (int i = 0; i < M; i++) { - for (int j = 0; j < N; j++) { - ptr[i + M * j] = this->get_coef(target_permutation[i + row_offset], source_permutation[j + col_offset]); - } - } - } else if (m_use_target_permutation) { - const auto &target_permutation = m_target_cluster.get_permutation(); - for (int i = 0; i < M; i++) { - for (int j = 0; j < N; j++) { - ptr[i + M * j] = this->get_coef(target_permutation[i + row_offset], j + col_offset); - } - } - } else if (m_use_source_permutation) { - const auto &source_permutation = m_source_cluster.get_permutation(); - for (int i = 0; i < M; i++) { - for (int j = 0; j < N; j++) { - ptr[i + M * j] = this->get_coef(i + row_offset, source_permutation[j + col_offset]); - } - } - } else { - for (int i = 0; i < M; i++) { - for (int j = 0; j < N; j++) { - ptr[i + M * j] = this->get_coef(i + row_offset, j + col_offset); - } + CoefficientPrecision operator()(int i, int j) { return get_coef(i, j); } + void copy_submatrix(int M, int N, const int *rows, const int *cols, CoefficientPrecision *ptr) const override { + for (int i = 0; i < M; i++) { + for (int j = 0; j < N; j++) { + ptr[i + M * j] = this->get_coef(rows[i], cols[j]); } } } +}; +// template > +// class GeneratorTest : public VirtualGenerator { +// protected: +// int m_target_size, m_source_size; +// const std::vector &m_target_points; +// const std::vector &m_source_points; +// const Cluster &m_target_cluster; +// const Cluster &m_source_cluster; +// bool m_use_target_permutation{true}; +// bool m_use_source_permutation{true}; +// int space_dim; - void set_use_target_permutation(bool use_target_permutation) { m_use_target_permutation = use_target_permutation; } - void set_use_source_permutation(bool use_source_permutation) { m_use_source_permutation = use_source_permutation; } +// public: +// explicit GeneratorTest(int space_dim0, int nr0, int nc0, const std::vector> &m_target_points0, const std::vector> &m_source_points0, const Cluster &target_cluster, const Cluster &source_cluster, bool use_target_permutation, bool use_source_permutation) : m_target_size(nr0), m_source_size(nc0), m_target_points(m_target_points0), m_source_points(m_source_points0), m_target_cluster(target_cluster), m_source_cluster(source_cluster), m_use_target_permutation(use_target_permutation), m_use_source_permutation(use_source_permutation), space_dim(space_dim0) {} - CoefficientPrecision operator()(int i, int j) { return get_coef(i, j); } +// virtual CoefficientPrecision get_coef(const int &i, const int &j) const = 0; - std::vector operator*(std::vector &a) const { - std::vector result(m_target_size, 0); - for (int i = 0; i < m_target_size; i++) { - for (int k = 0; k < m_source_size; k++) { - result[i] += this->get_coef(i, k) * a[k]; - } - } - return result; - } - double normFrob() { - double norm = 0; - for (int j = 0; j < this->nb_rows(); j++) { - for (int k = 0; k < this->nb_cols(); k++) { - norm = norm + std::pow(std::abs(this->get_coef(j, k)), 2); - } - } - return sqrt(norm); - } +// void copy_submatrix(int M, int N, int row_offset, int col_offset, CoefficientPrecision *ptr) const override { +// if (m_use_target_permutation && m_use_source_permutation) { +// const auto &target_permutation = m_target_cluster.get_permutation(); +// const auto &source_permutation = m_source_cluster.get_permutation(); +// for (int i = 0; i < M; i++) { +// for (int j = 0; j < N; j++) { +// ptr[i + M * j] = this->get_coef(target_permutation[i + row_offset], source_permutation[j + col_offset]); +// } +// } +// } else if (m_use_target_permutation) { +// const auto &target_permutation = m_target_cluster.get_permutation(); +// for (int i = 0; i < M; i++) { +// for (int j = 0; j < N; j++) { +// ptr[i + M * j] = this->get_coef(target_permutation[i + row_offset], j + col_offset); +// } +// } +// } else if (m_use_source_permutation) { +// const auto &source_permutation = m_source_cluster.get_permutation(); +// for (int i = 0; i < M; i++) { +// for (int j = 0; j < N; j++) { +// ptr[i + M * j] = this->get_coef(i + row_offset, source_permutation[j + col_offset]); +// } +// } +// } else { +// for (int i = 0; i < M; i++) { +// for (int j = 0; j < N; j++) { +// ptr[i + M * j] = this->get_coef(i + row_offset, j + col_offset); +// } +// } +// } +// } - void mvprod(const CoefficientPrecision *const in, CoefficientPrecision *const out, const int &mu) const { - int nr = m_target_size; - int nc = m_source_size; - for (int i = 0; i < nr * mu; i++) { - out[i] = 0; - } - for (int m = 0; m < mu; m++) { - for (int i = 0; i < nr; i++) { - for (int j = 0; j < nc; j++) { - out[nr * m + i] += this->get_coef(i, j) * in[j + m * nc]; - } - } - } - } +// void set_use_target_permutation(bool use_target_permutation) { m_use_target_permutation = use_target_permutation; } +// void set_use_source_permutation(bool use_source_permutation) { m_use_source_permutation = use_source_permutation; } - void mvprod_transp(const CoefficientPrecision *const in, CoefficientPrecision *const out, const int &mu) const { - int nc = m_target_size; - int nr = m_source_size; - for (int i = 0; i < nr * mu; i++) { - out[i] = 0; - } - for (int m = 0; m < mu; m++) { - for (int i = 0; i < nr; i++) { - for (int j = 0; j < nc; j++) { - out[nr * m + i] += this->get_coef(j, i) * in[j + m * nc]; - } - } - } - } +// CoefficientPrecision operator()(int i, int j) { return get_coef(i, j); } - void mvprod_conj(const CoefficientPrecision *const in, CoefficientPrecision *const out, const int &mu) const { - int nc = m_target_size; - int nr = m_source_size; - for (int i = 0; i < nr * mu; i++) { - out[i] = 0; - } - for (int m = 0; m < mu; m++) { - for (int i = 0; i < nr; i++) { - for (int j = 0; j < nc; j++) { - out[nr * m + i] += std::conj(this->get_coef(j, i) * in[j + m * nc]); - } - } - } - } -}; +// std::vector operator*(std::vector &a) const { +// std::vector result(m_target_size, 0); +// for (int i = 0; i < m_target_size; i++) { +// for (int k = 0; k < m_source_size; k++) { +// result[i] += this->get_coef(i, k) * a[k]; +// } +// } +// return result; +// } +// double normFrob() { +// double norm = 0; +// for (int j = 0; j < this->nb_rows(); j++) { +// for (int k = 0; k < this->nb_cols(); k++) { +// norm = norm + std::pow(std::abs(this->get_coef(j, k)), 2); +// } +// } +// return sqrt(norm); +// } + +// void mvprod(const CoefficientPrecision *const in, CoefficientPrecision *const out, const int &mu) const { +// int nr = m_target_size; +// int nc = m_source_size; +// for (int i = 0; i < nr * mu; i++) { +// out[i] = 0; +// } +// for (int m = 0; m < mu; m++) { +// for (int i = 0; i < nr; i++) { +// for (int j = 0; j < nc; j++) { +// out[nr * m + i] += this->get_coef(i, j) * in[j + m * nc]; +// } +// } +// } +// } + +// void mvprod_transp(const CoefficientPrecision *const in, CoefficientPrecision *const out, const int &mu) const { +// int nc = m_target_size; +// int nr = m_source_size; +// for (int i = 0; i < nr * mu; i++) { +// out[i] = 0; +// } +// for (int m = 0; m < mu; m++) { +// for (int i = 0; i < nr; i++) { +// for (int j = 0; j < nc; j++) { +// out[nr * m + i] += this->get_coef(j, i) * in[j + m * nc]; +// } +// } +// } +// } -class GeneratorTestDouble : public GeneratorTest { +// void mvprod_conj(const CoefficientPrecision *const in, CoefficientPrecision *const out, const int &mu) const { +// int nc = m_target_size; +// int nr = m_source_size; +// for (int i = 0; i < nr * mu; i++) { +// out[i] = 0; +// } +// for (int m = 0; m < mu; m++) { +// for (int i = 0; i < nr; i++) { +// for (int j = 0; j < nc; j++) { +// out[nr * m + i] += std::conj(this->get_coef(j, i) * in[j + m * nc]); +// } +// } +// } +// } +// }; + +class GeneratorTestDouble : public GeneratorTestWithPermutation { public: - using GeneratorTest::GeneratorTest; + using GeneratorTestWithPermutation::GeneratorTestWithPermutation; double get_coef(const int &i, const int &j) const override { - return 1. / (4 * M_PI * std::sqrt(std::inner_product(p1.begin() + this->space_dim * i, this->p1.begin() + this->space_dim * i + this->space_dim, p2.begin() + this->space_dim * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))); + return 1. / (4 * M_PI * std::sqrt(std::inner_product(m_target_points.begin() + m_space_dimension * i, this->m_target_points.begin() + m_space_dimension * i + m_space_dimension, m_source_points.begin() + m_space_dimension * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))); } }; -class GeneratorTestComplex : public GeneratorTest> { +class GeneratorTestComplex : public GeneratorTestWithPermutation> { public: - using GeneratorTest::GeneratorTest; + using GeneratorTestWithPermutation::GeneratorTestWithPermutation; std::complex get_coef(const int &i, const int &j) const override { - return (1. + std::complex(0, 1)) / (4 * M_PI * std::sqrt(std::inner_product(p1.begin() + this->space_dim * i, this->p1.begin() + this->space_dim * i + this->space_dim, p2.begin() + this->space_dim * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))); + return (1. + std::complex(0, 1)) / (4 * M_PI * std::sqrt(std::inner_product(m_target_points.begin() + m_space_dimension * i, this->m_target_points.begin() + m_space_dimension * i + m_space_dimension, m_source_points.begin() + m_space_dimension * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))); } }; @@ -153,41 +173,41 @@ double sign(double x) { return 0; } -class GeneratorTestDoubleSymmetric : public GeneratorTest { +class GeneratorTestDoubleSymmetric : public GeneratorTestWithPermutation { public: - using GeneratorTest::GeneratorTest; + using GeneratorTestWithPermutation::GeneratorTestWithPermutation; double get_coef(const int &i, const int &j) const override { - return 1. / (1e-5 + 4 * M_PI * std::sqrt(std::inner_product(p1.begin() + this->space_dim * i, this->p1.begin() + this->space_dim * i + this->space_dim, p2.begin() + this->space_dim * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))); + return 1. / (1e-5 + 4 * M_PI * std::sqrt(std::inner_product(m_target_points.begin() + m_space_dimension * i, this->m_target_points.begin() + m_space_dimension * i + m_space_dimension, m_source_points.begin() + m_space_dimension * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))); } }; -class GeneratorTestComplexSymmetric : public GeneratorTest> { +class GeneratorTestComplexSymmetric : public GeneratorTestWithPermutation> { public: - using GeneratorTest::GeneratorTest; + using GeneratorTestWithPermutation::GeneratorTestWithPermutation; std::complex get_coef(const int &i, const int &j) const override { - return (1. + std::complex(0, 1)) / (1e-5 + 4 * M_PI * std::sqrt(std::inner_product(p1.begin() + this->space_dim * i, this->p1.begin() + this->space_dim * i + this->space_dim, p2.begin() + this->space_dim * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))); + return (1. + std::complex(0, 1)) / (1e-5 + 4 * M_PI * std::sqrt(std::inner_product(m_target_points.begin() + m_space_dimension * i, this->m_target_points.begin() + m_space_dimension * i + m_space_dimension, m_source_points.begin() + m_space_dimension * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))); } }; -class GeneratorTestComplexHermitian : public GeneratorTest> { +class GeneratorTestComplexHermitian : public GeneratorTestWithPermutation> { public: - using GeneratorTest::GeneratorTest; + using GeneratorTestWithPermutation::GeneratorTestWithPermutation; std::complex get_coef(const int &i, const int &j) const override { - return (1. + sign(p1[this->space_dim * i] - p2[this->space_dim * j]) * std::complex(0, 1)) / (1e-5 + 4 * M_PI * std::sqrt(std::inner_product(p1.begin() + this->space_dim * i, this->p1.begin() + this->space_dim * i + this->space_dim, p2.begin() + this->space_dim * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))); + return (1. + sign(m_target_points[m_space_dimension * i] - m_source_points[m_space_dimension * j]) * std::complex(0, 1)) / (1e-5 + 4 * M_PI * std::sqrt(std::inner_product(m_target_points.begin() + m_space_dimension * i, this->m_target_points.begin() + m_space_dimension * i + m_space_dimension, m_source_points.begin() + m_space_dimension * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))); } }; template -class GeneratorFromMatrix : public VirtualGeneratorWithPermutation { +class GeneratorInUserNumberingFromMatrix : public VirtualGeneratorInUserNumbering { + public: const Matrix &A; - public: - explicit GeneratorFromMatrix(const Matrix &A0, const std::vector &target_permutation, const std::vector &source_permutation) : VirtualGeneratorWithPermutation(target_permutation.data(), source_permutation.data()), A(A0) {} + GeneratorInUserNumberingFromMatrix(const Matrix &A0) : A(A0) {} - void copy_submatrix_from_user_numbering(int M, int N, const int *const rows, const int *const cols, T *ptr) const override { + virtual void copy_submatrix(int M, int N, const int *rows, const int *cols, T *ptr) const override { for (int i = 0; i < M; i++) { for (int j = 0; j < N; j++) { ptr[i + M * j] = A(rows[i], cols[j]); @@ -196,16 +216,38 @@ class GeneratorFromMatrix : public VirtualGeneratorWithPermutation { } }; +// class GeneratorFromMatrix : public GeneratorWithPermutation { + +// public: +// explicit GeneratorFromMatrix(const Matrix &A0, const std::vector &target_permutation, const std::vector &source_permutation) : GeneratorWithPermutation(GeneratorInUserNumberingFromMatrix(A0), target_permutation.data(), source_permutation.data()) {} + +// protected: +// class GeneratorInUserNumberingFromMatrix : VirtualGeneratorInUserNumbering { +// public: +// const Matrix &A; + +// GeneratorInUserNumberingFromMatrix(const Matrix &A0) : A(A0) {} + +// virtual void copy_submatrix(int M, int N, const int *rows, const int *cols, T *ptr) { +// for (int i = 0; i < M; i++) { +// for (int j = 0; j < N; j++) { +// ptr[i + M * j] = A(rows[i], cols[j]); +// } +// } +// } +// }; +// }; + template -class LocalGeneratorFromMatrix : public VirtualGeneratorWithPermutation { +class LocalGeneratorInUserNumberingFromMatrix : public VirtualGeneratorInUserNumbering { const Matrix &m_A; const std::vector &m_target_local_to_global_numbering; const std::vector &m_source_local_to_global_numbering; public: - explicit LocalGeneratorFromMatrix(const Matrix &A, const std::vector &target_permutation, const std::vector &source_permutation, const std::vector &target_local_to_global_numbering, const std::vector &source_local_to_global_numbering) : VirtualGeneratorWithPermutation(target_permutation.data(), source_permutation.data()), m_A(A), m_target_local_to_global_numbering(target_local_to_global_numbering), m_source_local_to_global_numbering(source_local_to_global_numbering) {} + LocalGeneratorInUserNumberingFromMatrix(const Matrix &A, const std::vector &target_local_to_global_numbering, const std::vector &source_local_to_global_numbering) : m_A(A), m_target_local_to_global_numbering(target_local_to_global_numbering), m_source_local_to_global_numbering(source_local_to_global_numbering) {} - void copy_submatrix_from_user_numbering(int M, int N, const int *const rows, const int *const cols, T *ptr) const override { + void copy_submatrix(int M, int N, const int *const rows, const int *const cols, T *ptr) const override { for (int i = 0; i < M; i++) { for (int j = 0; j < N; j++) { ptr[i + M * j] = m_A(m_target_local_to_global_numbering[rows[i]], m_source_local_to_global_numbering[cols[j]]); @@ -213,6 +255,27 @@ class LocalGeneratorFromMatrix : public VirtualGeneratorWithPermutation { } } }; +// class LocalGeneratorFromMatrix : public GeneratorWithPermutation { +// const std::vector &m_target_local_to_global_numbering; +// const std::vector &m_source_local_to_global_numbering; + +// public: +// explicit LocalGeneratorFromMatrix(const Matrix &A, const std::vector &target_permutation, const std::vector &source_permutation, const std::vector &target_local_to_global_numbering, const std::vector &source_local_to_global_numbering) : GeneratorWithPermutation(LocalGeneratorInUserNumberingFromMatrix(A), target_permutation.data(), source_permutation.data()), m_target_local_to_global_numbering(target_local_to_global_numbering), m_source_local_to_global_numbering(source_local_to_global_numbering) {} + +// private: +// class LocalGeneratorInUserNumberingFromMatrix : VirtualGeneratorInUserNumbering { +// public: +// const Matrix &m_A; +// LocalGeneratorInUserNumberingFromMatrix(const Matrix &A) : m_A(A) {} +// void copy_submatrix_from(int M, int N, const int *const rows, const int *const cols, T *ptr) const override { +// for (int i = 0; i < M; i++) { +// for (int j = 0; j < N; j++) { +// ptr[i + M * j] = m_A(m_target_local_to_global_numbering[rows[i]], m_source_local_to_global_numbering[cols[j]]); +// } +// } +// } +// }; +// }; } // namespace htool #endif diff --git a/tests/functional_tests/distributed_operator/test_distributed_operator.hpp b/tests/functional_tests/distributed_operator/test_distributed_operator.hpp index 44d568b9..f5f42dd8 100644 --- a/tests/functional_tests/distributed_operator/test_distributed_operator.hpp +++ b/tests/functional_tests/distributed_operator/test_distributed_operator.hpp @@ -187,16 +187,15 @@ int test_vector_product(GeneratorTestType generator, const DistributedOperator -auto add_off_diagonal_operator(ClusterTreeBuilder> &recursive_build, DistributedOperator &distributed_operator, const Cluster> &target_root_cluster, const std::vector p1, const std::vector p1_permuted, const Cluster> &source_root_cluster, const std::vector p2_permuted, bool use_permutation, DataType data_type, htool::underlying_type epsilon, htool::underlying_type eta) { +auto add_off_diagonal_operator(ClusterTreeBuilder> &recursive_build, DistributedOperator &distributed_operator, const Cluster> &target_root_cluster, const std::vector p1, const Cluster> &source_root_cluster, const std::vector p2_permuted, DataType data_type, htool::underlying_type epsilon, htool::underlying_type eta) { struct Holder { std::unique_ptr>> off_diagonal_cluster; - // std::unique_ptr>> local_off_diagonal_cluster_tree_1, local_off_diagonal_cluster_tree_2; std::unique_ptr generator_off_diagonal; std::unique_ptr> off_diagonal_matrix_1, off_diagonal_matrix_2; std::unique_ptr>> off_diagonal_hmatrix_1, off_diagonal_hmatrix_2; std::unique_ptr>> local_off_diagonal_operator_1, local_off_diagonal_operator_2; - Holder(ClusterTreeBuilder> &recursive_build, const Cluster> &target_root_cluster, const std::vector p1, const std::vector p1_permuted, const Cluster> &source_root_cluster, const std::vector p2_permuted, bool use_permutation, DataType data_type, htool::underlying_type epsilon, htool::underlying_type eta) { + Holder(ClusterTreeBuilder> &recursive_build, const Cluster> &target_root_cluster, const std::vector &p1, const Cluster> &source_root_cluster, const std::vector &p2_permuted, DataType data_type, htool::underlying_type epsilon, htool::underlying_type eta) { // Local clusters int rankWorld; MPI_Comm_rank(MPI_COMM_WORLD, &rankWorld); @@ -204,8 +203,8 @@ auto add_off_diagonal_operator(ClusterTreeBuilder> &re const Cluster> &local_source_root_cluster = source_root_cluster.get_cluster_on_partition(rankWorld); // Sizes - int nc = source_root_cluster.get_size(); - int nr = target_root_cluster.get_size(); + int nc = source_root_cluster.get_size(); + // int nr = target_root_cluster.get_size(); int nc_local = local_source_root_cluster.get_size(); int off_diagonal_nc_1{local_source_root_cluster.get_offset()}; int off_diagonal_nc_2{source_root_cluster.get_size() - local_source_root_cluster.get_size() - local_source_root_cluster.get_offset()}; @@ -219,17 +218,20 @@ auto add_off_diagonal_operator(ClusterTreeBuilder> &re off_diagonal_partition.push_back(off_diagonal_nc_1 + nc_local); off_diagonal_partition.push_back(off_diagonal_nc_2); - // recursive_build.set_partition(2, off_diagonal_partition.data()); - off_diagonal_cluster = make_unique>>(recursive_build.create_cluster_tree(nc, 3, p2_permuted.data(), 2, 2, off_diagonal_partition.data())); + off_diagonal_cluster = make_unique>>(recursive_build.create_cluster_tree(nc, 3, p2_permuted.data(), 2, 2, off_diagonal_partition.data())); + generator_off_diagonal = std::make_unique(3, p1, p2_permuted); - // Generators - if (use_permutation) { - generator_off_diagonal = std::make_unique(3, nr, nc, p1, p2_permuted, local_target_root_cluster, *off_diagonal_cluster, true, true); - } else { - generator_off_diagonal = std::make_unique(3, nr, nc, p1_permuted, p2_permuted, local_target_root_cluster, *off_diagonal_cluster, true, true); - generator_off_diagonal->set_use_target_permutation(false); - generator_off_diagonal->set_use_source_permutation(true); - } + // // recursive_build.set_partition(2, off_diagonal_partition.data()); + // off_diagonal_cluster = make_unique>>(recursive_build.create_cluster_tree(nc, 3, p2_permuted.data(), 2, 2, off_diagonal_partition.data())); + + // // Generators + // if (use_permutation) { + // generator_off_diagonal = std::make_unique(3, nr, nc, p1, p2_permuted, local_target_root_cluster, *off_diagonal_cluster, true, true); + // } else { + // generator_off_diagonal = std::make_unique(3, nr, nc, p1_permuted, p2_permuted, local_target_root_cluster, *off_diagonal_cluster, true, true); + // generator_off_diagonal->set_use_target_permutation(false); + // generator_off_diagonal->set_use_source_permutation(true); + // } // Off diagonal LocalDenseMatrix const Cluster> *local_off_diagonal_cluster_tree_1 = &off_diagonal_cluster->get_cluster_on_partition(0); @@ -238,8 +240,11 @@ auto add_off_diagonal_operator(ClusterTreeBuilder> &re if (data_type == DataType::Matrix) { off_diagonal_matrix_1 = std::make_unique>(local_target_root_cluster.get_size(), local_off_diagonal_cluster_tree_1->get_size()); off_diagonal_matrix_2 = std::make_unique>(local_target_root_cluster.get_size(), local_off_diagonal_cluster_tree_2->get_size()); - generator_off_diagonal->copy_submatrix(local_target_root_cluster.get_size(), local_off_diagonal_cluster_tree_1->get_size(), local_target_root_cluster.get_offset(), local_off_diagonal_cluster_tree_1->get_offset(), off_diagonal_matrix_1->data()); - generator_off_diagonal->copy_submatrix(local_target_root_cluster.get_size(), local_off_diagonal_cluster_tree_2->get_size(), local_target_root_cluster.get_offset(), local_off_diagonal_cluster_tree_2->get_offset(), off_diagonal_matrix_2->data()); + + generator_off_diagonal->copy_submatrix(local_target_root_cluster.get_size(), local_off_diagonal_cluster_tree_1->get_size(), local_target_root_cluster.get_permutation().data() + local_target_root_cluster.get_offset(), local_off_diagonal_cluster_tree_1->get_permutation().data() + local_off_diagonal_cluster_tree_1->get_offset(), off_diagonal_matrix_1->data()); + + generator_off_diagonal->copy_submatrix(local_target_root_cluster.get_size(), local_off_diagonal_cluster_tree_2->get_size(), local_target_root_cluster.get_permutation().data() + local_target_root_cluster.get_offset(), local_off_diagonal_cluster_tree_2->get_permutation().data() + local_off_diagonal_cluster_tree_2->get_offset(), off_diagonal_matrix_2->data()); + local_off_diagonal_operator_1 = make_unique>>(*off_diagonal_matrix_1, local_target_root_cluster, *local_off_diagonal_cluster_tree_1, 'N', 'N', false, true); local_off_diagonal_operator_2 = make_unique>>(*off_diagonal_matrix_2, local_target_root_cluster, *local_off_diagonal_cluster_tree_2, 'N', 'N', false, true); @@ -256,7 +261,7 @@ auto add_off_diagonal_operator(ClusterTreeBuilder> &re } }; - Holder holder(recursive_build, target_root_cluster, p1, p1_permuted, source_root_cluster, p2_permuted, use_permutation, data_type, epsilon, eta); + Holder holder(recursive_build, target_root_cluster, p1, source_root_cluster, p2_permuted, data_type, epsilon, eta); if (holder.off_diagonal_cluster->get_cluster_on_partition(0).get_size() > 0) distributed_operator.add_local_operator(holder.local_off_diagonal_operator_1.get()); if (holder.off_diagonal_cluster->get_cluster_on_partition(1).get_size() > 0) @@ -264,146 +269,6 @@ auto add_off_diagonal_operator(ClusterTreeBuilder> &re return holder; } -template -bool test_custom_distributed_operator(int nr, int nc, int mu, bool use_permutation, char Symmetry, char UPLO, char op, bool off_diagonal_approximation, DataType data_type, htool::underlying_type epsilon = 1e-14) { - - // Get the number of processes - int sizeWorld; - MPI_Comm_size(MPI_COMM_WORLD, &sizeWorld); - - // Get the rankWorld of the process - int rankWorld; - MPI_Comm_rank(MPI_COMM_WORLD, &rankWorld); - - srand(1); - bool test = 0; - htool::underlying_type eta = 10; - - // Geometry - double z1 = 1; - vector p1(3 * nr), p1_permuted, off_diagonal_p1; - vector p2(Symmetry == 'N' ? 3 * nc : 1), p2_permuted, off_diagonal_p2; - create_disk(3, z1, nr, p1.data()); - int size_numbering = nr / (sizeWorld); - int count_size = 0; - std::vector MasterOffset_target, MasterOffset_source; - for (int p = 0; p < sizeWorld - 1; p++) { - MasterOffset_target.push_back(count_size); - MasterOffset_target.push_back(size_numbering); - - count_size += size_numbering; - } - MasterOffset_target.push_back(count_size); - MasterOffset_target.push_back(nr - count_size); - - size_numbering = nc / sizeWorld; - count_size = 0; - - ClusterTreeBuilder> recursive_build; - std::shared_ptr>> source_root_cluster; - std::shared_ptr>> target_root_cluster = make_shared>>(recursive_build.create_cluster_tree(nr, 3, p1.data(), 2, sizeWorld, MasterOffset_target.data())); - - if (Symmetry == 'N') { - double z2 = 1 + 0.1; - create_disk(3, z2, nc, p2.data()); - - for (int p = 0; p < sizeWorld - 1; p++) { - MasterOffset_source.push_back(count_size); - MasterOffset_source.push_back(size_numbering); - - count_size += size_numbering; - } - MasterOffset_source.push_back(count_size); - MasterOffset_source.push_back(nc - count_size); - - // recursive_build.set_partition(sizeWorld, MasterOffset_source.data()); - - source_root_cluster = make_shared>>(recursive_build.create_cluster_tree(nc, 3, p2.data(), 2, sizeWorld, MasterOffset_source.data())); - - } else { - - MasterOffset_source = MasterOffset_target; - source_root_cluster = target_root_cluster; - p2 = p1; - } - - // Permutation on geometry - p1_permuted.resize(3 * nr); - const auto &target_permutation = target_root_cluster->get_permutation(); - for (int i = 0; i < target_permutation.size(); i++) { - p1_permuted[i * 3 + 0] = p1[target_permutation[i] * 3 + 0]; - p1_permuted[i * 3 + 1] = p1[target_permutation[i] * 3 + 1]; - p1_permuted[i * 3 + 2] = p1[target_permutation[i] * 3 + 2]; - } - p2_permuted.resize(3 * nc); - if (Symmetry == 'N') { - const auto &source_permutation = source_root_cluster->get_permutation(); - for (int i = 0; i < source_permutation.size(); i++) { - p2_permuted[i * 3 + 0] = p2[source_permutation[i] * 3 + 0]; - p2_permuted[i * 3 + 1] = p2[source_permutation[i] * 3 + 1]; - p2_permuted[i * 3 + 2] = p2[source_permutation[i] * 3 + 2]; - } - } else { - p2_permuted = p1_permuted; - } - - // Generator - GeneratorTestType generator(3, nr, nc, p1, p2, *target_root_cluster, *source_root_cluster, true, true); - GeneratorTestType generator_permuted(3, nr, nc, p1_permuted, p2_permuted, *target_root_cluster, *source_root_cluster, false, false); - - // Diagonal LocalDenseMatrix - std::shared_ptr>> local_operator; - const Cluster> *actual_target_cluster = &target_root_cluster->get_cluster_on_partition(rankWorld); - const Cluster> *actual_source_cluster = off_diagonal_approximation ? &source_root_cluster->get_cluster_on_partition(rankWorld) : source_root_cluster.get(); - - char symmetry = (sizeWorld == 1 || off_diagonal_approximation) ? Symmetry : 'N'; - char uplo = (sizeWorld == 1 || off_diagonal_approximation) ? UPLO : 'N'; - - std::unique_ptr>> hmatrix; - std::unique_ptr> matrix; - - if (data_type == DataType::Matrix) { - matrix = std::make_unique>(actual_target_cluster->get_size(), actual_source_cluster->get_size()); - if (symmetry == 'N') { - generator_permuted.copy_submatrix(matrix->nb_rows(), matrix->nb_cols(), actual_target_cluster->get_offset(), actual_source_cluster->get_offset(), matrix->data()); - } else if ((symmetry == 'S' || symmetry == 'H') && uplo == 'L') { - for (int i = 0; i < matrix->nb_rows(); i++) { - for (int j = 0; j < i + 1; j++) { - generator_permuted.copy_submatrix(1, 1, i + actual_target_cluster->get_offset(), j + actual_source_cluster->get_offset(), matrix->data() + i + j * matrix->nb_rows()); - } - } - } else if ((symmetry == 'S' || symmetry == 'H') && uplo == 'U') { - for (int j = 0; j < matrix->nb_cols(); j++) { - for (int i = 0; i < j + 1; i++) { - generator_permuted.copy_submatrix(1, 1, i + actual_target_cluster->get_offset(), j + actual_source_cluster->get_offset(), matrix->data() + i + j * matrix->nb_rows()); - } - } - } - local_operator = std::make_unique>>(*matrix, *actual_target_cluster, *actual_source_cluster, symmetry, uplo); - } else if (data_type == DataType::HMatrix) { - HMatrixTreeBuilder> hmatrix_builder(*actual_target_cluster, *actual_source_cluster, epsilon, eta, symmetry, uplo, -1, -1, rankWorld); - hmatrix = std::make_unique>>(hmatrix_builder.build(generator_permuted)); - local_operator = std::make_unique>>(*hmatrix, *actual_target_cluster, *actual_source_cluster, symmetry, uplo); - } - - // Distributed operator - PartitionFromCluster> target_partition(*target_root_cluster); - PartitionFromCluster> source_partition(*source_root_cluster); - DistributedOperator distributed_operator(target_partition, source_partition, Symmetry, UPLO, MPI_COMM_WORLD); - distributed_operator.add_local_operator(local_operator.get()); - distributed_operator.use_permutation() = use_permutation; - - if (off_diagonal_approximation) { - auto dependencies = add_off_diagonal_operator(recursive_build, distributed_operator, *target_root_cluster, p1, p1_permuted, *source_root_cluster, p2_permuted, use_permutation, data_type, epsilon, eta); - - test = test_vector_product(generator, distributed_operator, *target_root_cluster, MasterOffset_target, *source_root_cluster, MasterOffset_source, mu, op, use_permutation, epsilon); - } else { - test = test_vector_product(generator, distributed_operator, *target_root_cluster, MasterOffset_target, *source_root_cluster, MasterOffset_source, mu, op, use_permutation, epsilon); - } - - return test; -} - template bool test_default_distributed_operator(int nr, int nc, int mu, bool use_permutation, char Symmetry, char UPLO, char op, bool off_diagonal_approximation, DataType data_type, htool::underlying_type epsilon = 1e-14) { @@ -488,23 +353,22 @@ bool test_default_distributed_operator(int nr, int nc, int mu, bool use_permutat } // Generator - GeneratorTestType generator(3, nr, nc, p1, p2, *target_root_cluster, *source_root_cluster, true, true); - GeneratorTestType generator_permuted(3, nr, nc, p1_permuted, p2_permuted, *target_root_cluster, *source_root_cluster, false, false); + GeneratorTestType generator_in_user_numbering(3, p1, p2); if (off_diagonal_approximation) { - DefaultLocalApproximationBuilder> distributed_operator_holder(generator_permuted, *target_root_cluster, *source_root_cluster, epsilon, eta, Symmetry, UPLO, MPI_COMM_WORLD); + DefaultLocalApproximationBuilder> distributed_operator_holder(generator_in_user_numbering, *target_root_cluster, *source_root_cluster, epsilon, eta, Symmetry, UPLO, MPI_COMM_WORLD); DistributedOperator &distributed_operator = distributed_operator_holder.distributed_operator; distributed_operator.use_permutation() = use_permutation; - auto dependencies = add_off_diagonal_operator(recursive_build, distributed_operator, *target_root_cluster, p1, p1_permuted, *source_root_cluster, p2_permuted, use_permutation, data_type, epsilon, eta); + auto dependencies = add_off_diagonal_operator(recursive_build, distributed_operator, *target_root_cluster, p1, *source_root_cluster, p2_permuted, data_type, epsilon, eta); - test = test_vector_product(generator, distributed_operator, *target_root_cluster, MasterOffset_target, *source_root_cluster, MasterOffset_source, mu, op, use_permutation, epsilon); + test = test_vector_product(generator_in_user_numbering, distributed_operator, *target_root_cluster, MasterOffset_target, *source_root_cluster, MasterOffset_source, mu, op, use_permutation, epsilon); } else { - DefaultApproximationBuilder> distributed_operator_holder(generator_permuted, *target_root_cluster, *source_root_cluster, epsilon, eta, Symmetry, UPLO, MPI_COMM_WORLD); + DefaultApproximationBuilder> distributed_operator_holder(generator_in_user_numbering, *target_root_cluster, *source_root_cluster, epsilon, eta, Symmetry, UPLO, MPI_COMM_WORLD); DistributedOperator &distributed_operator = distributed_operator_holder.distributed_operator; distributed_operator.use_permutation() = use_permutation; - test = test_vector_product(generator, distributed_operator, *target_root_cluster, MasterOffset_target, *source_root_cluster, MasterOffset_source, mu, op, use_permutation, epsilon); + test = test_vector_product(generator_in_user_numbering, distributed_operator, *target_root_cluster, MasterOffset_target, *source_root_cluster, MasterOffset_source, mu, op, use_permutation, epsilon); } return test; @@ -512,8 +376,5 @@ bool test_default_distributed_operator(int nr, int nc, int mu, bool use_permutat template bool test_distributed_operator(int nr, int nc, int mu, bool use_permutation, char Symmetry, char UPLO, char op, bool off_diagonal_approximation, DataType data_type, htool::underlying_type epsilon) { - if (data_type == DataType::DefaultHMatrix) - return test_default_distributed_operator(nr, nc, mu, use_permutation, Symmetry, UPLO, op, off_diagonal_approximation, data_type, epsilon); - else - return test_custom_distributed_operator(nr, nc, mu, use_permutation, Symmetry, UPLO, op, off_diagonal_approximation, data_type, epsilon); + return test_default_distributed_operator(nr, nc, mu, use_permutation, Symmetry, UPLO, op, off_diagonal_approximation, data_type, epsilon); } diff --git a/tests/functional_tests/distributed_operator/test_distributed_operator_product_complex_double.cpp b/tests/functional_tests/distributed_operator/test_distributed_operator_product_complex_double.cpp index 24581b02..8417003d 100644 --- a/tests/functional_tests/distributed_operator/test_distributed_operator_product_complex_double.cpp +++ b/tests/functional_tests/distributed_operator/test_distributed_operator_product_complex_double.cpp @@ -16,9 +16,9 @@ int main(int argc, char *argv[]) { for (auto offdiagonal_approximation : {true, false}) { for (auto use_permutation : {true, false}) { for (auto operation : {'N', 'T'}) { - for (auto data_type : {DataType::Matrix, DataType::HMatrix, DataType::DefaultHMatrix}) { + for (auto data_type : {DataType::Matrix, DataType::HMatrix}) { std::vector tolerances{1e-14}; - if (data_type == DataType::HMatrix || data_type == DataType::DefaultHMatrix) { + if (data_type == DataType::HMatrix) { tolerances.push_back(1e-3); } for (auto epsilon : tolerances) { diff --git a/tests/functional_tests/distributed_operator/test_distributed_operator_product_double.cpp b/tests/functional_tests/distributed_operator/test_distributed_operator_product_double.cpp index 76f97ce6..8b7f4a20 100644 --- a/tests/functional_tests/distributed_operator/test_distributed_operator_product_double.cpp +++ b/tests/functional_tests/distributed_operator/test_distributed_operator_product_double.cpp @@ -16,9 +16,9 @@ int main(int argc, char *argv[]) { for (auto offdiagonal_approximation : {true, false}) { for (auto use_permutation : {true, false}) { for (auto operation : {'N', 'T'}) { - for (auto data_type : {DataType::Matrix, DataType::HMatrix, DataType::DefaultHMatrix}) { + for (auto data_type : {DataType::Matrix, DataType::HMatrix}) { std::vector tolerances{1e-14}; - if (data_type == DataType::HMatrix || data_type == DataType::DefaultHMatrix) { + if (data_type == DataType::HMatrix) { tolerances.push_back(1e-3); } for (auto epsilon : tolerances) { diff --git a/tests/functional_tests/hmatrix/hmatrix_addition/test_hmatrix_addition_complex_double.cpp b/tests/functional_tests/hmatrix/hmatrix_addition/test_hmatrix_addition_complex_double.cpp index 9dfd7673..b68598ef 100644 --- a/tests/functional_tests/hmatrix/hmatrix_addition/test_hmatrix_addition_complex_double.cpp +++ b/tests/functional_tests/hmatrix/hmatrix_addition/test_hmatrix_addition_complex_double.cpp @@ -5,9 +5,8 @@ using namespace std; using namespace htool; -int main(int argc, char *argv[]) { +int main(int, char *[]) { - MPI_Init(&argc, &argv); bool is_error = false; const double margin = 10; @@ -19,7 +18,6 @@ int main(int argc, char *argv[]) { } } } - MPI_Finalize(); if (is_error) { return 1; } diff --git a/tests/functional_tests/hmatrix/hmatrix_addition/test_hmatrix_addition_double.cpp b/tests/functional_tests/hmatrix/hmatrix_addition/test_hmatrix_addition_double.cpp index 26028f6a..5657099f 100644 --- a/tests/functional_tests/hmatrix/hmatrix_addition/test_hmatrix_addition_double.cpp +++ b/tests/functional_tests/hmatrix/hmatrix_addition/test_hmatrix_addition_double.cpp @@ -5,9 +5,8 @@ using namespace std; using namespace htool; -int main(int argc, char *argv[]) { +int main(int, char *[]) { - MPI_Init(&argc, &argv); bool is_error = false; const double margin = 10; @@ -19,7 +18,6 @@ int main(int argc, char *argv[]) { } } } - MPI_Finalize(); if (is_error) { return 1; } diff --git a/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_complex_double.cpp b/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_complex_double.cpp index c8489a33..b0552133 100644 --- a/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_complex_double.cpp +++ b/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_complex_double.cpp @@ -7,9 +7,8 @@ using namespace std; using namespace htool; -int main(int argc, char *argv[]) { +int main(int, char *[]) { - MPI_Init(&argc, &argv); bool is_error = false; const int n1 = 500; const double margin = 1; @@ -20,7 +19,7 @@ int main(int argc, char *argv[]) { is_error = is_error || test_hmatrix_lu, GeneratorTestComplexHermitian>(trans, n1, n2, epsilon, margin); } } - MPI_Finalize(); + if (is_error) { return 1; } diff --git a/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_double.cpp b/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_double.cpp index 614b4488..f68037f4 100644 --- a/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_double.cpp +++ b/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_double.cpp @@ -7,9 +7,8 @@ using namespace std; using namespace htool; -int main(int argc, char *argv[]) { +int main(int, char *[]) { - MPI_Init(&argc, &argv); bool is_error = false; const int n1 = 1000; const double margin = 1; @@ -23,7 +22,7 @@ int main(int argc, char *argv[]) { // is_error = is_error || test_hmatrix_cholesky(UPLO, n1, n2, epsilon, margin); // } } - MPI_Finalize(); + if (is_error) { return 1; } diff --git a/tests/functional_tests/hmatrix/hmatrix_product/test_hmatrix_product_complex_double.cpp b/tests/functional_tests/hmatrix/hmatrix_product/test_hmatrix_product_complex_double.cpp index 2b081370..fe2fd0e3 100644 --- a/tests/functional_tests/hmatrix/hmatrix_product/test_hmatrix_product_complex_double.cpp +++ b/tests/functional_tests/hmatrix/hmatrix_product/test_hmatrix_product_complex_double.cpp @@ -1,5 +1,6 @@ #include "../test_hmatrix_product.hpp" #include +#include using namespace std; using namespace htool; diff --git a/tests/functional_tests/hmatrix/hmatrix_triangular_solve/test_hmatrix_triangular_solve_complex_double.cpp b/tests/functional_tests/hmatrix/hmatrix_triangular_solve/test_hmatrix_triangular_solve_complex_double.cpp index 53a3fddd..bf2a0d45 100644 --- a/tests/functional_tests/hmatrix/hmatrix_triangular_solve/test_hmatrix_triangular_solve_complex_double.cpp +++ b/tests/functional_tests/hmatrix/hmatrix_triangular_solve/test_hmatrix_triangular_solve_complex_double.cpp @@ -7,9 +7,8 @@ using namespace std; using namespace htool; -int main(int argc, char *argv[]) { +int main(int, char *[]) { - MPI_Init(&argc, &argv); bool is_error = false; const int n1 = 400; const double margin = 1; @@ -25,7 +24,6 @@ int main(int argc, char *argv[]) { } } } - MPI_Finalize(); if (is_error) { return 1; } diff --git a/tests/functional_tests/hmatrix/hmatrix_triangular_solve/test_hmatrix_triangular_solve_double.cpp b/tests/functional_tests/hmatrix/hmatrix_triangular_solve/test_hmatrix_triangular_solve_double.cpp index 744a36b2..fdce0cb2 100644 --- a/tests/functional_tests/hmatrix/hmatrix_triangular_solve/test_hmatrix_triangular_solve_double.cpp +++ b/tests/functional_tests/hmatrix/hmatrix_triangular_solve/test_hmatrix_triangular_solve_double.cpp @@ -7,9 +7,8 @@ using namespace std; using namespace htool; -int main(int argc, char *argv[]) { +int main(int, char *[]) { - MPI_Init(&argc, &argv); bool is_error = false; const int n1 = 400; const double margin = 1; @@ -24,7 +23,6 @@ int main(int argc, char *argv[]) { } } } - MPI_Finalize(); if (is_error) { return 1; } diff --git a/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_SVD.cpp b/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_SVD.cpp index e1774ec4..688fc1b1 100644 --- a/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_SVD.cpp +++ b/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_SVD.cpp @@ -41,7 +41,8 @@ int main(int argc, char *argv[]) { Cluster t = recursive_build_strategy.create_cluster_tree(nr, 3, xt.data(), 2, 2); Cluster s = recursive_build_strategy.create_cluster_tree(nc, 3, xs.data(), 2, 2); - GeneratorTestDouble A(3, nr, nc, xt, xs, t, s, true, true); + GeneratorTestDouble A_in_user_numbering(3, xt, xs); + GeneratorWithPermutation A(A_in_user_numbering, t.get_permutation().data(), s.get_permutation().data()); // SVD fixed rank int reqrank_max = 10; diff --git a/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_fullACA.cpp b/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_fullACA.cpp index 550d8058..4e1c0e53 100644 --- a/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_fullACA.cpp +++ b/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_fullACA.cpp @@ -40,7 +40,8 @@ int main(int argc, char *argv[]) { Cluster target_cluster = recursive_build_strategy.create_cluster_tree(nr, 3, xt.data(), 2, 2); Cluster source_cluster = recursive_build_strategy.create_cluster_tree(nc, 3, xt.data(), 2, 2); - GeneratorTestDouble A(3, nr, nc, xt, xs, target_cluster, source_cluster, true, true); + GeneratorTestDouble A_in_user_numbering(3, xt, xs); + GeneratorWithPermutation A(A_in_user_numbering, target_cluster.get_permutation().data(), source_cluster.get_permutation().data()); // fullACA fixed rank int reqrank_max = 10; diff --git a/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_partialACA.cpp b/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_partialACA.cpp index c0154197..cd144363 100644 --- a/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_partialACA.cpp +++ b/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_partialACA.cpp @@ -41,7 +41,8 @@ int main(int argc, char *argv[]) { Cluster t = recursive_build_strategy.create_cluster_tree(nr, 3, xt.data(), 2, 2); Cluster s = recursive_build_strategy.create_cluster_tree(nc, 3, xt.data(), 2, 2); - GeneratorTestDouble A(3, nr, nc, xt, xs, t, s, true, true); + GeneratorTestDouble A_in_user_numbering(3, xt, xs); + GeneratorWithPermutation A(A_in_user_numbering, t.get_permutation().data(), s.get_permutation().data()); // partialACA fixed rank int reqrank_max = 10; diff --git a/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_sympartialACA.cpp b/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_sympartialACA.cpp index 3d7fe6fd..cda07a26 100644 --- a/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_sympartialACA.cpp +++ b/tests/functional_tests/hmatrix/lrmat/lrmat_build/test_lrmat_build_sympartialACA.cpp @@ -42,7 +42,8 @@ int main(int argc, char *argv[]) { Cluster t = recursive_build_strategy.create_cluster_tree(nr, 3, xt.data(), 2, 2); Cluster s = recursive_build_strategy.create_cluster_tree(nc, 3, xt.data(), 2, 2); - GeneratorTestDouble A(3, nr, nc, xt, xs, t, s, true, true); + GeneratorTestDouble A_in_user_numbering(3, xt, xs); + GeneratorWithPermutation A(A_in_user_numbering, t.get_permutation().data(), s.get_permutation().data()); // sympartialACA fixed rank int reqrank_max = 10; diff --git a/tests/functional_tests/hmatrix/lrmat/test_lrmat_build.hpp b/tests/functional_tests/hmatrix/lrmat/test_lrmat_build.hpp index a2a9997a..c856a182 100644 --- a/tests/functional_tests/hmatrix/lrmat/test_lrmat_build.hpp +++ b/tests/functional_tests/hmatrix/lrmat/test_lrmat_build.hpp @@ -14,7 +14,7 @@ using namespace std; using namespace htool; template -bool test_lrmat(const Cluster &target_cluster, const Cluster &source_cluster, const GeneratorTestDouble &A, const LowRankMatrix &Fixed_approximation, const LowRankMatrix &Auto_approximation, std::pair fixed_compression_interval, std::pair auto_compression_interval) { +bool test_lrmat(const Cluster &target_cluster, const Cluster &source_cluster, const VirtualGenerator &A, const LowRankMatrix &Fixed_approximation, const LowRankMatrix &Auto_approximation, std::pair fixed_compression_interval, std::pair auto_compression_interval) { bool test = 0; diff --git a/tests/functional_tests/hmatrix/test_hmatrix_build.hpp b/tests/functional_tests/hmatrix/test_hmatrix_build.hpp index 324cad80..001874ed 100644 --- a/tests/functional_tests/hmatrix/test_hmatrix_build.hpp +++ b/tests/functional_tests/hmatrix/test_hmatrix_build.hpp @@ -14,7 +14,7 @@ using namespace std; using namespace htool; -template +template bool test_hmatrix_build(int nr, int nc, bool use_local_cluster, char Symmetry, char UPLO, htool::underlying_type epsilon, bool use_dense_blocks_generator) { // Get the number of processes @@ -39,12 +39,12 @@ bool test_hmatrix_build(int nr, int nc, bool use_local_cluster, char Symmetry, c test_partition(3, nr, p1, sizeWorld, partition); // Clustering - ClusterTreeBuilder> recursive_build_strategy; + ClusterTreeBuilder> cluster_tree_builder; // recursive_build_strategy.set_partition(partition); // recursive_build_strategy.set_minclustersize(2); std::shared_ptr>> source_root_cluster; - std::shared_ptr>> target_root_cluster = make_shared>>(recursive_build_strategy.create_cluster_tree(nr, 3, p1.data(), 2, sizeWorld, partition.data())); + std::shared_ptr>> target_root_cluster = make_shared>>(cluster_tree_builder.create_cluster_tree(nr, 3, p1.data(), 2, sizeWorld, partition.data())); if (Symmetry == 'N' && nr != nc) { // Geometry @@ -57,34 +57,14 @@ bool test_hmatrix_build(int nr, int nc, bool use_local_cluster, char Symmetry, c // Clustering // source_recursive_build_strategy.set_minclustersize(2); - source_root_cluster = make_shared>>(recursive_build_strategy.create_cluster_tree(nc, 3, p2.data(), 2, sizeWorld, partition.data())); + source_root_cluster = make_shared>>(cluster_tree_builder.create_cluster_tree(nc, 3, p2.data(), 2, sizeWorld, partition.data())); } else { source_root_cluster = target_root_cluster; p2 = p1; } - // Permutation on geometry - p1_permuted.resize(3 * nr); - const auto &target_permutation = target_root_cluster->get_permutation(); - for (int i = 0; i < target_permutation.size(); i++) { - p1_permuted[i * 3 + 0] = p1[target_permutation[i] * 3 + 0]; - p1_permuted[i * 3 + 1] = p1[target_permutation[i] * 3 + 1]; - p1_permuted[i * 3 + 2] = p1[target_permutation[i] * 3 + 2]; - } - p2_permuted.resize(3 * nc); - if (Symmetry == 'N' && nr != nc) { - const auto &source_permutation = source_root_cluster->get_permutation(); - for (int i = 0; i < source_permutation.size(); i++) { - p2_permuted[i * 3 + 0] = p2[source_permutation[i] * 3 + 0]; - p2_permuted[i * 3 + 1] = p2[source_permutation[i] * 3 + 1]; - p2_permuted[i * 3 + 2] = p2[source_permutation[i] * 3 + 2]; - } - } else { - p2_permuted = p1_permuted; - } - - // Generator - GeneratorTestType generator(3, nr, nc, p1_permuted, p2_permuted, *target_root_cluster, *source_root_cluster, false, false); + GeneratorTestTypeInUserNumbering generator(3, p1, p2); + GeneratorWithPermutation generator_with_permutation(generator, target_root_cluster->get_permutation().data(), source_root_cluster->get_permutation().data()); // HMatrix double eta = 10; @@ -98,7 +78,7 @@ bool test_hmatrix_build(int nr, int nc, bool use_local_cluster, char Symmetry, c std::shared_ptr> dense_blocks_generator; if (use_dense_blocks_generator) { - dense_blocks_generator = std::make_shared>(generator); + dense_blocks_generator = std::make_shared>(generator_with_permutation); } hmatrix_tree_builder->set_dense_blocks_generator(dense_blocks_generator); @@ -141,7 +121,7 @@ bool test_hmatrix_build(int nr, int nc, bool use_local_cluster, char Symmetry, c Matrix dense_matrix; dense_matrix.assign(hmatrix_target_cluster.get_size(), hmatrix_source_cluster.get_size(), dense_data.data(), false); - generator.copy_submatrix(hmatrix_target_cluster.get_size(), hmatrix_source_cluster.get_size(), hmatrix_target_cluster.get_offset(), hmatrix_source_cluster.get_offset(), dense_data.data()); + generator_with_permutation.copy_submatrix(hmatrix_target_cluster.get_size(), hmatrix_source_cluster.get_size(), hmatrix_target_cluster.get_offset(), hmatrix_source_cluster.get_offset(), dense_data.data()); // if (rankWorld == 0) { // std::cout << "dense :\n"; @@ -197,40 +177,14 @@ bool test_hmatrix_build(int nr, int nc, bool use_local_cluster, char Symmetry, c 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 &target_permutation = target_root_cluster->get_permutation(); 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); diff --git a/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp b/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp index 2329189b..771b9089 100644 --- a/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp +++ b/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp @@ -47,7 +47,7 @@ bool test_hmatrix_lu(char trans, int n1, int n2, htool::underlying_type epsil template bool test_hmatrix_cholesky(char UPLO, int n1, int n2, htool::underlying_type epsilon, htool::underlying_type margin) { bool is_error = false; - double eta = 10; + double eta = -100; htool::underlying_type error; // Setup test case diff --git a/tests/functional_tests/hmatrix/test_hmatrix_hmatrix_product.hpp b/tests/functional_tests/hmatrix/test_hmatrix_hmatrix_product.hpp index 639bdf50..6cceda92 100644 --- a/tests/functional_tests/hmatrix/test_hmatrix_hmatrix_product.hpp +++ b/tests/functional_tests/hmatrix/test_hmatrix_hmatrix_product.hpp @@ -11,6 +11,7 @@ #include #include #include +#include using namespace std; using namespace htool;