diff --git a/CHANGELOG.md b/CHANGELOG.md index a1d7ee41..c238d6fd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,8 @@ All notable changes to this project will be documented in this file. ## Unreleased +## [0.9.0] - 2024-TBA + ### Added - The old implementation of `HMatrix` was mixing the distributed operations and compression via hierarchical matrices. This is fixed by replacing `HMatrix` by: @@ -33,12 +35,13 @@ All notable changes to this project will be documented in this file. - `VirtualLocalOperator` which is the interface local operators must satisfy, - `LocalDenseMatrix` is an example of local operator consisting of a dense matrix `Matrix`, - and `LocalHMatrix` which is an example of local operator consisting of a hierarchical matrix based on `HMatrix` (different from the previous `HMatrix`, see below). - ### Changed - `HMatrix` is now a class representing a hierarchical matrix without distributed-memory parallelism (note that it can still use shared-memory parallelism): - - It inherits from `TreeNode`, and it provides the algebra related to hierchical matrices. + - It inherits from `TreeNode`, and it provides the algebra related to hierchical matrices: + - product with vector and matrix (threaded with OpenMP), + - and with this new version, LU and Cholesky factorization (not threaded yet, WIP). - The algorithms for building the block cluster tree is contained in `HMatrixTreeBuilder`. Users can provide their own "factory". - `VirtualCluster` is removed and the clustering part of the library has been rewritten: - `Cluster` now derives from `TreeNode`, whose template parameter corresponds to the precision of cluster nodes' radius and centre (previously only `double`). diff --git a/tests/functional_tests/multilrmat/CMakeLists.txt b/tests/functional_tests/multilrmat/CMakeLists.txt deleted file mode 100755 index 657f0cbf..00000000 --- a/tests/functional_tests/multilrmat/CMakeLists.txt +++ /dev/null @@ -1,15 +0,0 @@ -#=============================================================================# -#=========================== Executables =====================================# -#=============================================================================# - -#=== lrmat_SVD -set(compressions "partialACA") -# list(APPEND compressions "partialACA") -# list(APPEND compressions "sympartialACA") -# list(APPEND compressions "SVD") -foreach(compression ${compressions}) - add_executable(Test_multi_lrmat_${compression} test_multi_lrmat_${compression}.cpp) - target_link_libraries(Test_multi_lrmat_${compression} htool) - add_dependencies(build-tests Test_multi_lrmat_${compression}) - add_test(Test_multi_lrmat_${compression} Test_multi_lrmat_${compression}) -endforeach() diff --git a/tests/functional_tests/multilrmat/test_multi_lrmat.hpp b/tests/functional_tests/multilrmat/test_multi_lrmat.hpp deleted file mode 100644 index c4c1c722..00000000 --- a/tests/functional_tests/multilrmat/test_multi_lrmat.hpp +++ /dev/null @@ -1,146 +0,0 @@ -#include -#include -#include -#include - -#include -#include -#include -#include - -using namespace std; -using namespace htool; -class MyMultiMatrix : public VirtualMultiGenerator { - const vector &p1; - const vector &p2; - int space_dim; - - public: - MyMultiMatrix(int space_dim0, int nr, int nc, int nm, const vector &p10, const vector &p20) : VirtualMultiGenerator(nr, nc, nm), p1(p10), p2(p20), space_dim(space_dim0) {} - - std::vector get_coefs(const int &i, const int &j) const { - return std::vector{ - (1.) / (4 * M_PI * std::sqrt(std::inner_product(p1.begin() + space_dim * i, p1.begin() + space_dim * i + space_dim, p2.begin() + space_dim * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))), - (2.) / (4 * M_PI * std::sqrt(std::inner_product(p1.begin() + space_dim * i, p1.begin() + space_dim * i + space_dim, p2.begin() + space_dim * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))), - (3.) / (4 * M_PI * std::sqrt(std::inner_product(p1.begin() + space_dim * i, p1.begin() + space_dim * i + space_dim, p2.begin() + space_dim * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))), - (4.) / (4 * M_PI * std::sqrt(std::inner_product(p1.begin() + space_dim * i, p1.begin() + space_dim * i + space_dim, p2.begin() + space_dim * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); }))), - (5.) / (4 * M_PI * std::sqrt(std::inner_product(p1.begin() + space_dim * i, p1.begin() + space_dim * i + space_dim, p2.begin() + space_dim * j, double(0), std::plus(), [](double u, double v) { return (u - v) * (u - v); })))}; - } - - void copy_submatrices(int M, int N, const int *const rows, const int *const cols, int NM, double *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]); - } - } - } - std::vector mult(std::vector &a, int l) const { - std::vector result(this->nr, 0); - for (int i = 0; i < this->nr; i++) { - for (int k = 0; k < this->nc; k++) { - result[i] += this->get_coefs(i, k)[l] * a[k]; - } - } - return result; - } - - double normFrob(int l) const { - 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_coefs(j, k))[l]), 2); - } - } - return sqrt(norm); - } -}; - -template -int test_multi_lrmat(const MyMultiMatrix &A, const MultiLowRankMatrix &Fixed_approximation, const MultiLowRankMatrix &Auto_approximation, const std::vector &permt, const std::vector &perms, std::pair fixed_compression_interval, std::pair auto_compression_interval) { - - bool test = 0; - int nr = permt.size(); - int nc = perms.size(); - - // Random vector - double lower_bound = 0; - double upper_bound = 10000; - std::random_device rd; - std::mt19937 mersenne_engine(rd()); - std::uniform_real_distribution dist(lower_bound, upper_bound); - auto gen = [&dist, &mersenne_engine]() { - return dist(mersenne_engine); - }; - - vector f(nc, 1); - generate(begin(f), end(f), gen); - - // ACA with fixed rank - int reqrank_max = 10; - std::vector> fixed_errors; - for (int k = 0; k < reqrank_max + 1; k++) { - fixed_errors.push_back(Frobenius_absolute_error(Fixed_approximation, A, k)); - } - - // Test rank - test = test || !(Fixed_approximation.rank_of() == reqrank_max); - std::cout << test << " " << reqrank_max << " " << Fixed_approximation.rank_of() << std::endl; - cout << "Compression with fixed rank" << endl; - cout << "> rank : " << Fixed_approximation.rank_of() << endl; - - // Test Frobenius errors - test = test || !(max(fixed_errors.back()) < 1e-7); - cout << "> Errors with Frobenius norm : " << fixed_errors << endl; - - for (int l = 0; l < A.nb_matrix(); l++) { - - // Test compression - test = test || !(fixed_compression_interval.first < Fixed_approximation[l].space_saving() && Fixed_approximation[l].space_saving() < fixed_compression_interval.second); - cout << "> Compression rate : " << Fixed_approximation[l].space_saving() << endl; - - // Test mat vec prod - std::vector out_perm(nr); - std::vector out = Fixed_approximation[l] * f; - for (int i = 0; i < permt.size(); i++) { - out_perm[permt[i]] = out[i]; - } - double error = norm2(A.mult(f, l) - out_perm) / norm2(A.mult(f, l)); - test = test || !(error < Fixed_approximation[l].get_epsilon() * 10); - cout << "> Errors on a mat vec prod : " << error << " " << (Fixed_approximation[l].get_epsilon() * 10) << " " << (error < Fixed_approximation[l].get_epsilon() * 10) << endl; - cout << "test : " << test << endl - << endl; - } - - // ACA automatic building - std::vector> auto_errors; - for (int k = 0; k < Auto_approximation.rank_of() + 1; k++) { - auto_errors.push_back(Frobenius_absolute_error(Auto_approximation, A, k)); - } - cout << "Automatic compression" << endl; - - // Test Frobenius error - test = test || !(max(auto_errors[Auto_approximation.rank_of()]) < Auto_approximation.get_epsilon()); - cout << "> Errors with Frobenius norm: " << auto_errors << endl; - - for (int l = 0; l < A.nb_matrix(); l++) { - - // Test compression rate - test = test || !(auto_compression_interval.first < Auto_approximation[l].space_saving() && Auto_approximation[l].space_saving() < auto_compression_interval.second); - cout << "> Compression rate : " << Auto_approximation[l].space_saving() << endl; - - // Test mat vec prod - std::vector out_perm(nr); - std::vector out = Auto_approximation[l] * f; - for (int i = 0; i < permt.size(); i++) { - out_perm[permt[i]] = out[i]; - } - double error = norm2(A.mult(f, l) - out_perm) / norm2(A.mult(f, l)); - test = test || !(error < Auto_approximation[l].get_epsilon() * 10); - cout << "> Errors on a mat vec prod : " << error << endl; - - cout << "test : " << test << endl - << endl; - } - - return test; -} diff --git a/tests/functional_tests/multilrmat/test_multi_lrmat_partialACA.cpp b/tests/functional_tests/multilrmat/test_multi_lrmat_partialACA.cpp deleted file mode 100644 index 815b8c99..00000000 --- a/tests/functional_tests/multilrmat/test_multi_lrmat_partialACA.cpp +++ /dev/null @@ -1,86 +0,0 @@ -#include -#include -#include - -#include "test_multi_lrmat.hpp" -#include -#include - -using namespace std; -using namespace htool; - -int main(int argc, char *argv[]) { - - // Initialize the MPI environment - MPI_Init(&argc, &argv); - - const int ndistance = 4; - double distance[ndistance]; - distance[0] = 15; - distance[1] = 20; - distance[2] = 30; - distance[3] = 40; - - double epsilon = 0.0001; - - int nr = 500; - int nc = 100; - - std::vector xt(3 * nr); - std::vector xs(3 * nc); - bool test = 0; - - for (int idist = 0; idist < ndistance; idist++) { - - create_disk(3, 0, nr, xt.data()); - create_disk(3, distance[idist], nc, xs.data()); - - Cluster t, s; - - t.build(nr, xt.data()); - s.build(nc, xs.data()); - - MyMultiMatrix A(3, nr, nc, xt, xs); - int nm = A.nb_matrix(); - GeneratorTestDouble A_test(3, nr, nc, xt, xs); - - // partialACA fixed rank - int reqrank_max = 10; - MultiLowRankMatrix A_partialACA_fixed(A_test.get_dimension(), t.get_perm(), s.get_perm(), nm, reqrank_max, epsilon); - LowRankMatrix A_partialACA_fixed_test(A_test.get_dimension(), t.get_perm(), s.get_perm(), reqrank_max, epsilon); - A_partialACA_fixed.build(A, MultipartialACA(), t, xt.data(), s, xs.data()); - ; - A_partialACA_fixed_test.build(A_test, partialACA(), t, xt.data(), s, xs.data()); - ; - - // // ACA automatic building - // MultiLowRankMatrix A_partialACA(A_test.get_dimension(), t.get_perm(), s.get_perm(), nm); - // A_partialACA.set_epsilon(epsilon); - // A_partialACA.build(A, MultipartialACA(), t, xt.data(), s, xs.data()); - // LowRankMatrix A_partialACA_test(A_test.get_dimension(), t.get_perm(), s.get_perm()); - // A_partialACA_test.set_epsilon(epsilon); - // A_partialACA_test.build(A_test, partialACA(), t, xt.data(), s, xs.data()); - // ; - - // // Comparison with lrmat - // std::vector one(nc, 1); - // test = test || !(norm2(A_partialACA_fixed[0] * one - A_partialACA_fixed_test * one) < 1e-10); - // cout << "> Errors for fixed rank compared to lrmat: " << norm2(A_partialACA_fixed[0] * one - A_partialACA_fixed_test * one) << endl; - - // test = test || !(norm2(A_partialACA[0] * one - A_partialACA_test * one) < 1e-10); - // cout << "> Errors for auto rank compared to lrmat: " << norm2(A_partialACA[0] * one - A_partialACA_test * one) << endl; - - // // Test multi lrmat - // std::pair fixed_compression_interval(0.87, 0.89); - // std::pair auto_compression_interval(0.93, 0.96); - - // test = test || (test_multi_lrmat(A, A_partialACA_fixed, A_partialACA, t.get_perm(), s.get_perm(), fixed_compression_interval, auto_compression_interval)); - } - - cout << "test : " << test << endl; - - // Finalize the MPI environment. - MPI_Finalize(); - - return test; -} diff --git a/tests/test_warnings.cpp b/tests/test_warnings.cpp deleted file mode 100644 index a6f5a220..00000000 --- a/tests/test_warnings.cpp +++ /dev/null @@ -1,7 +0,0 @@ -#define OMPI_SKIP_MPICXX 1 -#define MPICH_SKIP_MPICXX 1 -#include - -int main(int, char *[]) { - return 0; -} \ No newline at end of file