diff --git a/.gitignore b/.gitignore index 023bd72..d963fdc 100644 --- a/.gitignore +++ b/.gitignore @@ -16,7 +16,6 @@ dist/ downloads/ eggs/ .eggs/ -lib/ lib64/ parts/ sdist/ diff --git a/example/define_custom_low_rank_generator.py b/example/define_custom_low_rank_generator.py index f4f9123..a8d263b 100644 --- a/example/define_custom_low_rank_generator.py +++ b/example/define_custom_low_rank_generator.py @@ -16,10 +16,16 @@ def build_low_rank_approximation(self, rows, cols, epsilon): norm = np.linalg.norm(submat) svd_norm = 0 - count = len(s) - 1 - while count > 0 and math.sqrt(svd_norm) / norm < epsilon: - svd_norm += s[count] ** 2 - count -= 1 - count = min(count + 1, min(len(rows), len(cols))) - self.set_U(u[:, 0:count] * s[0:count]) - self.set_V(vh[0:count, :]) + truncated_rank = len(s) - 1 + while truncated_rank > 0 and math.sqrt(svd_norm) / norm < epsilon: + svd_norm += s[truncated_rank] ** 2 + truncated_rank -= 1 + truncated_rank += 1 + + if truncated_rank * (len(rows) + len(cols)) > (len(rows) * len(cols)): + print("coucou") + return False # the low rank approximation is not worthwhile + + self.set_U(u[:, 0:truncated_rank] * s[0:truncated_rank]) + self.set_V(vh[0:truncated_rank, :]) + return True diff --git a/example/use_block_jacobi_solver.py b/example/use_block_jacobi_solver.py index eb9c277..a338d5b 100644 --- a/example/use_block_jacobi_solver.py +++ b/example/use_block_jacobi_solver.py @@ -1,4 +1,5 @@ import copy +import logging import Htool import matplotlib.pyplot as plt @@ -7,6 +8,8 @@ from create_geometry import create_random_geometries from define_custom_generators import CustomGenerator +logging.basicConfig(level=logging.INFO) + # Random geometry size = 500 dimension = 3 @@ -41,9 +44,14 @@ "L", mpi4py.MPI.COMM_WORLD, ) +hmatrix = default_approximation.hmatrix +Htool.recompression(hmatrix) # Solver with block Jacobi preconditionner -block_diagonal_hmatrix = copy.deepcopy(default_approximation.block_diagonal_hmatrix) +block_diagonal_hmatrix = copy.deepcopy( + default_approximation.block_diagonal_hmatrix +) # inplace factorization requires deepcopy + default_solver_builder = Htool.DDMSolverBuilder( default_approximation.distributed_operator, block_diagonal_hmatrix ) @@ -64,7 +72,6 @@ # Several ways to display information -hmatrix = default_approximation.hmatrix hmatrix_distributed_information = hmatrix.get_distributed_information( mpi4py.MPI.COMM_WORLD ) diff --git a/example/use_custom_dense_block_generator.py b/example/use_custom_dense_block_generator.py index b7055df..64a6b78 100644 --- a/example/use_custom_dense_block_generator.py +++ b/example/use_custom_dense_block_generator.py @@ -1,9 +1,14 @@ +import logging + import Htool import mpi4py import numpy as np from create_geometry import create_random_geometries from define_custom_dense_blocks_generator import CustomDenseBlocksGenerator from define_custom_generators import CustomGenerator +from matplotlib import pyplot as plt + +logging.basicConfig(level=logging.INFO) # Random geometry size = 500 @@ -46,6 +51,7 @@ distributed_operator = distributed_operator_from_hmatrix.distributed_operator hmatrix = distributed_operator_from_hmatrix.hmatrix +Htool.openmp_recompression(hmatrix) # Test matrix vector product np.random.seed(0) @@ -60,3 +66,35 @@ Y_1 = distributed_operator @ X Y_2 = generator.mat_mat(X) print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(Y_1 - Y_2) / np.linalg.norm(Y_2)) +# Several ways to display information +hmatrix_distributed_information = hmatrix.get_distributed_information( + mpi4py.MPI.COMM_WORLD +) +hmatrix_tree_parameter = hmatrix.get_tree_parameters() +hmatrix_local_information = hmatrix.get_local_information() +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + print(hmatrix_distributed_information) + print(hmatrix_local_information) + print(hmatrix_tree_parameter) + + fig = plt.figure() + + if dimension == 2: + ax1 = fig.add_subplot(2, 2, 1) + ax2 = fig.add_subplot(2, 2, 2) + # ax3 = fig.add_subplot(2, 2, 3) + ax4 = fig.add_subplot(2, 2, 4) + elif dimension == 3: + ax1 = fig.add_subplot(2, 2, 1, projection="3d") + ax2 = fig.add_subplot(2, 2, 2, projection="3d") + # ax3 = fig.add_subplot(2, 2, 3, projection="3d") + ax4 = fig.add_subplot(2, 2, 4) + + ax1.set_title("target cluster at depth 1") + ax2.set_title("target cluster at depth 2") + # ax3.set_title("source cluster at depth 1") + ax4.set_title("Hmatrix on rank 0") + Htool.plot(ax1, target_cluster, points, 1) + Htool.plot(ax2, target_cluster, points, 2) + Htool.plot(ax4, hmatrix) + plt.show() diff --git a/example/use_custom_local_operator.py b/example/use_custom_local_operator.py index d9814fe..b97d031 100644 --- a/example/use_custom_local_operator.py +++ b/example/use_custom_local_operator.py @@ -15,8 +15,6 @@ # Htool parameters -eta = 10 -epsilon = 1e-3 minclustersize = 10 number_of_children = 2 @@ -42,7 +40,7 @@ "N", "N", False, - True, + False, ) # Build distributed operator diff --git a/example/use_custom_low_rank_approximation.py b/example/use_custom_low_rank_approximation.py index 1168828..6565822 100644 --- a/example/use_custom_low_rank_approximation.py +++ b/example/use_custom_low_rank_approximation.py @@ -1,9 +1,14 @@ +import logging + import Htool import mpi4py import numpy as np from create_geometry import create_partitionned_geometries from define_custom_generators import CustomGenerator from define_custom_low_rank_generator import CustomSVD +from matplotlib import pyplot as plt + +logging.basicConfig(level=logging.INFO) # Random geometry nb_rows = 500 @@ -48,7 +53,8 @@ ) distributed_operator = distributed_operator_from_hmatrix.distributed_operator - +hmatrix = distributed_operator_from_hmatrix.hmatrix +Htool.openmp_recompression(hmatrix) # Test matrix vector product np.random.seed(0) @@ -63,3 +69,37 @@ Y_1 = distributed_operator @ X Y_2 = generator.mat_mat(X) print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(Y_1 - Y_2) / np.linalg.norm(Y_2)) + +# Several ways to display information +hmatrix_distributed_information = hmatrix.get_distributed_information( + mpi4py.MPI.COMM_WORLD +) +hmatrix_tree_parameter = hmatrix.get_tree_parameters() +hmatrix_local_information = hmatrix.get_local_information() +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + print(hmatrix_distributed_information) + print(hmatrix_local_information) + print(hmatrix_tree_parameter) + + fig = plt.figure() + + if dimension == 2: + ax1 = fig.add_subplot(2, 2, 1) + ax2 = fig.add_subplot(2, 2, 2) + ax3 = fig.add_subplot(2, 2, 3) + ax4 = fig.add_subplot(2, 2, 4) + elif dimension == 3: + ax1 = fig.add_subplot(2, 2, 1, projection="3d") + ax2 = fig.add_subplot(2, 2, 2, projection="3d") + ax3 = fig.add_subplot(2, 2, 3, projection="3d") + ax4 = fig.add_subplot(2, 2, 4) + + ax1.set_title("target cluster at depth 1") + ax2.set_title("target cluster at depth 2") + ax3.set_title("source cluster at depth 1") + ax4.set_title("Hmatrix on rank 0") + Htool.plot(ax1, target_cluster, target_points, 1) + Htool.plot(ax2, target_cluster, target_points, 2) + Htool.plot(ax3, source_cluster, source_points, 1) + Htool.plot(ax4, hmatrix) + plt.show() diff --git a/example/use_hmatrix_compression.py b/example/use_hmatrix_compression.py index 932a118..60b6fc5 100644 --- a/example/use_hmatrix_compression.py +++ b/example/use_hmatrix_compression.py @@ -1,3 +1,5 @@ +import logging + import Htool import matplotlib.pyplot as plt import mpi4py @@ -5,6 +7,8 @@ from create_geometry import create_partitionned_geometries from define_custom_generators import CustomGenerator +logging.basicConfig(level=logging.INFO) + # Random geometry nb_rows = 500 nb_cols = 500 @@ -47,6 +51,9 @@ ) distributed_operator = default_approximation.distributed_operator +hmatrix = default_approximation.hmatrix +Htool.openmp_recompression(hmatrix) + # Test matrix vector product np.random.seed(0) @@ -63,7 +70,6 @@ print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(Y_1 - Y_2) / np.linalg.norm(Y_2)) # Several ways to display information -hmatrix = default_approximation.hmatrix hmatrix_distributed_information = hmatrix.get_distributed_information( mpi4py.MPI.COMM_WORLD ) diff --git a/example/use_local_hmatrix_compression.py b/example/use_local_hmatrix_compression.py index 60d9e23..10e7387 100644 --- a/example/use_local_hmatrix_compression.py +++ b/example/use_local_hmatrix_compression.py @@ -1,3 +1,5 @@ +import logging + import Htool import matplotlib.pyplot as plt import mpi4py @@ -6,6 +8,8 @@ from define_custom_generators import CustomGenerator from define_custom_local_operator import CustomLocalOperator +logging.basicConfig(level=logging.INFO) + # Random geometry target_size = 500 source_size = 500 @@ -68,6 +72,8 @@ mpi4py.MPI.COMM_WORLD, ) distributed_operator = default_local_approximation.distributed_operator +hmatrix = default_local_approximation.hmatrix +Htool.recompression(hmatrix) # Build off diagonal operators @@ -114,9 +120,9 @@ True, ) -if off_diagonal_nc_1 > 0: +if local_operator_1: distributed_operator.add_local_operator(local_operator_1) -if off_diagonal_nc_2 > 0: +if local_operator_2: distributed_operator.add_local_operator(local_operator_2) # Test matrix vector product @@ -135,7 +141,6 @@ print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(Y_1 - Y_2) / np.linalg.norm(Y_2)) # Several ways to display information -hmatrix = default_local_approximation.hmatrix hmatrix_distributed_information = hmatrix.get_distributed_information( mpi4py.MPI.COMM_WORLD ) diff --git a/lib/htool b/lib/htool index 2b8bb26..c9fd344 160000 --- a/lib/htool +++ b/lib/htool @@ -1 +1 @@ -Subproject commit 2b8bb26695de4051bd9d72b02bc3fc9b368edddd +Subproject commit c9fd34469fb9cb66a7b7aa67bede6c6a6159906c diff --git a/src/htool/distributed_operator/utility.hpp b/src/htool/distributed_operator/utility.hpp index 3eceecb..13f3034 100644 --- a/src/htool/distributed_operator/utility.hpp +++ b/src/htool/distributed_operator/utility.hpp @@ -27,8 +27,8 @@ void declare_distributed_operator_utility(py::module &m, std::string prefix = "" default_approximation_class.def(py::init &, const Cluster &, const Cluster &, htool::underlying_type, htool::underlying_type, char, char, MPI_Comm_wrapper>()); default_approximation_class.def_property_readonly( "distributed_operator", [](const DefaultApproximation &self) { return &self.distributed_operator; }, py::return_value_policy::reference_internal); - default_approximation_class.def_property_readonly( - "hmatrix", [](const DefaultApproximation &self) { return &self.hmatrix; }, py::return_value_policy::reference_internal); + default_approximation_class.def_property( + "hmatrix", [](const DefaultApproximation &self) { return &self.hmatrix; }, [](DefaultApproximation &self, const HMatrix &hmatrix) { self.hmatrix = hmatrix; }); default_approximation_class.def_property_readonly( "block_diagonal_hmatrix", [](const DefaultApproximation &self) { return &*self.block_diagonal_hmatrix; }, py::return_value_policy::reference_internal); @@ -36,8 +36,8 @@ void declare_distributed_operator_utility(py::module &m, std::string prefix = "" default_local_approximation_class.def(py::init &, const Cluster &, const Cluster &, htool::underlying_type, htool::underlying_type, char, char, MPI_Comm_wrapper>()); default_local_approximation_class.def_property_readonly( "distributed_operator", [](const LocalDefaultApproximation &self) { return &self.distributed_operator; }, py::return_value_policy::reference_internal); - default_local_approximation_class.def_property_readonly( - "hmatrix", [](const LocalDefaultApproximation &self) { return &self.hmatrix; }, py::return_value_policy::reference_internal); + default_local_approximation_class.def_readwrite( + "hmatrix", &LocalDefaultApproximation::hmatrix); default_local_approximation_class.def_property_readonly( "block_diagonal_hmatrix", [](const LocalDefaultApproximation &self) { return &*self.block_diagonal_hmatrix; }, py::return_value_policy::reference_internal); @@ -45,8 +45,8 @@ void declare_distributed_operator_utility(py::module &m, std::string prefix = "" distributed_operator_from_hmatrix_class.def(py::init &, const Cluster &, const Cluster &, const HMatrixTreeBuilder &, MPI_Comm_wrapper>()); distributed_operator_from_hmatrix_class.def_property_readonly( "distributed_operator", [](const DistributedOperatorFromHMatrix &self) { return &self.distributed_operator; }, py::return_value_policy::reference_internal); - distributed_operator_from_hmatrix_class.def_property_readonly( - "hmatrix", [](const DistributedOperatorFromHMatrix &self) { return &self.hmatrix; }, py::return_value_policy::reference_internal); + distributed_operator_from_hmatrix_class.def_readwrite( + "hmatrix", &DistributedOperatorFromHMatrix::hmatrix); distributed_operator_from_hmatrix_class.def_property_readonly( "block_diagonal_hmatrix", [](const DistributedOperatorFromHMatrix &self) { return &*self.block_diagonal_hmatrix; }, py::return_value_policy::reference_internal); } diff --git a/src/htool/hmatrix/hmatrix.hpp b/src/htool/hmatrix/hmatrix.hpp index 899f134..d5237cb 100644 --- a/src/htool/hmatrix/hmatrix.hpp +++ b/src/htool/hmatrix/hmatrix.hpp @@ -1,6 +1,7 @@ #ifndef HTOOL_HMATRIX_CPP #define HTOOL_HMATRIX_CPP +#include #include #include #include @@ -11,6 +12,7 @@ #include #include #include +#include namespace py = pybind11; using namespace pybind11::literals; @@ -43,6 +45,11 @@ void declare_HMatrix(py::module &m, const std::string &className) { py_class.def("get_tree_parameters", [](const HMatrix &hmatrix) { return htool::get_tree_parameters(hmatrix); }); py_class.def("get_local_information", [](const HMatrix &hmatrix) { return htool::get_hmatrix_information(hmatrix); }); py_class.def("get_distributed_information", [](const HMatrix &hmatrix, MPI_Comm_wrapper comm) { return htool::get_distributed_hmatrix_information(hmatrix, comm); }); + + m.def("recompression", &htool::recompression &)>>); + m.def("recompression", [](HMatrix &hmatrix) { recompression(hmatrix); }); + m.def("openmp_recompression", &htool::openmp_recompression &)>>); + m.def("openmp_recompression", [](HMatrix &hmatrix) { recompression(hmatrix); }); } #endif diff --git a/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp b/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp index 3bab9b3..c87cb21 100644 --- a/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp +++ b/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp @@ -21,19 +21,29 @@ class VirtualLowRankGeneratorPython : public VirtualLowRankGenerator epsilon, int &rank, Matrix &U, Matrix &V) const override { + bool copy_low_rank_approximation(int M, int N, const int *const rows, const int *const cols, LowRankMatrix &lrmat) const override { + auto &U = lrmat.get_U(); + auto &V = lrmat.get_V(); py::array_t py_rows(std::array{M}, rows, py::capsule(rows)); py::array_t py_cols(std::array{N}, cols, py::capsule(cols)); - build_low_rank_approximation(py_rows, py_cols, epsilon); - U.assign(m_mats_U.back().shape()[0], m_mats_U.back().shape()[1], m_mats_U.back().mutable_data(), false); - V.assign(m_mats_V.back().shape()[0], m_mats_V.back().shape()[1], m_mats_V.back().mutable_data(), false); + bool success = build_low_rank_approximation(py_rows, py_cols, lrmat.get_epsilon()); + if (success) { + U.assign(m_mats_U.back().shape()[0], m_mats_U.back().shape()[1], m_mats_U.back().mutable_data(), false); + V.assign(m_mats_V.back().shape()[0], m_mats_V.back().shape()[1], m_mats_V.back().mutable_data(), false); + } + return success; + } + + bool copy_low_rank_approximation(int M, int N, const int *const rows, const int *const cols, int reqrank, LowRankMatrix &lrmat) const override { + Logger::get_instance().log(LogLevel::ERROR, "copy_low_rank_approximation with required rank is not implemented in the python interface."); + return false; } bool is_htool_owning_data() const override { return false; } // lcov does not see it because of trampoline I assume - virtual void build_low_rank_approximation(const py::array_t &rows, const py::array_t &cols, underlying_type epsilon) const = 0; // LCOV_EXCL_LINE + virtual bool build_low_rank_approximation(const py::array_t &rows, const py::array_t &cols, underlying_type epsilon) const = 0; // LCOV_EXCL_LINE void set_U(py::array_t U0) { m_mats_U.push_back(U0); // no copy here @@ -47,9 +57,9 @@ class PyVirtualLowRankGenerator : public VirtualLowRankGeneratorPython::VirtualLowRankGeneratorPython; /* Trampoline (need one for each virtual function) */ - virtual void build_low_rank_approximation(const py::array_t &rows, const py::array_t &cols, underlying_type epsilon) const override { + virtual bool build_low_rank_approximation(const py::array_t &rows, const py::array_t &cols, underlying_type epsilon) const override { PYBIND11_OVERRIDE_PURE( - void, /* Return type */ + bool, /* Return type */ VirtualLowRankGeneratorPython, /* Parent class */ build_low_rank_approximation, /* Name of function in C++ (must match Python name) */ rows, diff --git a/src/htool/hmatrix/lrmat.hpp b/src/htool/hmatrix/lrmat.hpp new file mode 100644 index 0000000..a074da3 --- /dev/null +++ b/src/htool/hmatrix/lrmat.hpp @@ -0,0 +1,19 @@ +#ifndef HTOOL_LRMAT_CPP +#define HTOOL_LRMAT_CPP + +#include +#include + +namespace py = pybind11; +using namespace htool; + +template +void declare_LowRankMatrix(py::module &m, const std::string &className) { + using Class = LowRankMatrix; + py::class_ py_class(m, className.c_str()); + + py_class.def("nb_rows", &Class::nb_rows); + py_class.def("nb_cols", &Class::nb_cols); + py_class.def("rank", &Class::rank_of, "test"); +} +#endif diff --git a/src/htool/main.cpp b/src/htool/main.cpp index 7f94e5a..e9b688b 100644 --- a/src/htool/main.cpp +++ b/src/htool/main.cpp @@ -15,6 +15,7 @@ #include "hmatrix/interfaces/virtual_dense_blocks_generator.hpp" #include "hmatrix/interfaces/virtual_generator.hpp" #include "hmatrix/interfaces/virtual_low_rank_generator.hpp" +#include "hmatrix/lrmat.hpp" #include "local_operator/local_operator.hpp" #include "local_operator/virtual_local_operator.hpp" @@ -55,6 +56,7 @@ PYBIND11_MODULE(Htool, m) { declare_geometric_splitting(m); declare_hmatrix_builder(m, "HMatrixBuilder"); + declare_LowRankMatrix(m, "LowRankMatrix"); declare_HMatrix(m, "HMatrix"); declare_virtual_generator(m, "VirtualGenerator", "IGenerator"); declare_custom_VirtualLowRankGenerator(m, "VirtualLowRankGenerator"); @@ -77,6 +79,7 @@ PYBIND11_MODULE(Htool, m) { declare_hmatrix_builder, double>(m, "ComplexHMatrixBuilder"); declare_HMatrix, double>(m, "ComplexHMatrix"); + declare_LowRankMatrix>(m, "ComplexLowRankMatrix"); declare_virtual_generator>(m, "ComplexVirtualGenerator", "IComplexGenerator"); declare_custom_VirtualLowRankGenerator>(m, "VirtualComplexLowRankGenerator");