update

.gitignore

@@ -16,7 +16,6 @@ dist/
 downloads/
 eggs/
 .eggs/
-lib/
 lib64/
 parts/
 sdist/

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

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
 )

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()

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

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()

example/use_hmatrix_compression.py

@@ -1,10 +1,14 @@
+import logging
+
 import Htool
 import matplotlib.pyplot as plt
 import mpi4py
 import numpy as np
 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
 )

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
 )

src/htool/distributed_operator/utility.hpp

@@ -27,26 +27,26 @@ void declare_distributed_operator_utility(py::module &m, std::string prefix = ""
     default_approximation_class.def(py::init<const VirtualGenerator<CoefficientPrecision> &, const Cluster<CoordinatePrecision> &, const Cluster<CoordinatePrecision> &, htool::underlying_type<CoefficientPrecision>, htool::underlying_type<CoefficientPrecision>, 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<CoefficientPrecision, CoordinatePrecision> &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);
 
     py::class_<LocalDefaultApproximation> default_local_approximation_class(m, default_local_approximation_name.c_str());
     default_local_approximation_class.def(py::init<const VirtualGenerator<CoefficientPrecision> &, const Cluster<CoordinatePrecision> &, const Cluster<CoordinatePrecision> &, htool::underlying_type<CoefficientPrecision>, htool::underlying_type<CoefficientPrecision>, 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);
 
     py::class_<DistributedOperatorFromHMatrix> distributed_operator_from_hmatrix_class(m, distributed_operator_from_hmatrix_name.c_str());
     distributed_operator_from_hmatrix_class.def(py::init<const VirtualGenerator<CoefficientPrecision> &, const Cluster<CoordinatePrecision> &, const Cluster<CoordinatePrecision> &, const HMatrixTreeBuilder<CoefficientPrecision, CoordinatePrecision> &, 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);
 }

src/htool/hmatrix/hmatrix.hpp

@@ -1,6 +1,7 @@
 #ifndef HTOOL_HMATRIX_CPP
 #define HTOOL_HMATRIX_CPP
 
+#include <pybind11/functional.h>
 #include <pybind11/numpy.h>
 #include <pybind11/operators.h>
 #include <pybind11/pybind11.h>
@@ -11,6 +12,7 @@
 #include <htool/hmatrix/hmatrix.hpp>
 #include <htool/hmatrix/hmatrix_distributed_output.hpp>
 #include <htool/hmatrix/hmatrix_output.hpp>
+#include <htool/hmatrix/utils/recompression.hpp>
 
 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<CoefficientPrecision, CoordinatePrecision> &hmatrix) { return htool::get_tree_parameters(hmatrix); });
     py_class.def("get_local_information", [](const HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix) { return htool::get_hmatrix_information(hmatrix); });
     py_class.def("get_distributed_information", [](const HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix, MPI_Comm_wrapper comm) { return htool::get_distributed_hmatrix_information(hmatrix, comm); });
+
+    m.def("recompression", &htool::recompression<CoefficientPrecision, CoordinatePrecision, std::function<void(LowRankMatrix<CoefficientPrecision> &)>>);
+    m.def("recompression", [](HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix) { recompression(hmatrix); });
+    m.def("openmp_recompression", &htool::openmp_recompression<CoefficientPrecision, CoordinatePrecision, std::function<void(LowRankMatrix<CoefficientPrecision> &)>>);
+    m.def("openmp_recompression", [](HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix) { recompression(hmatrix); });
 }
 
 #endif

src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp

@@ -21,19 +21,29 @@ class VirtualLowRankGeneratorPython : public VirtualLowRankGenerator<Coefficient
 
     VirtualLowRankGeneratorPython() {}
 
-    void copy_low_rank_approximation(int M, int N, const int *const rows, const int *const cols, underlying_type<CoefficientPrecision> epsilon, int &rank, Matrix<CoefficientPrecision> &U, Matrix<CoefficientPrecision> &V) const override {
+    bool copy_low_rank_approximation(int M, int N, const int *const rows, const int *const cols, LowRankMatrix<CoefficientPrecision> &lrmat) const override {
+        auto &U = lrmat.get_U();
+        auto &V = lrmat.get_V();
         py::array_t<int> py_rows(std::array<long int, 1>{M}, rows, py::capsule(rows));
         py::array_t<int> py_cols(std::array<long int, 1>{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<CoefficientPrecision> &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<int, py::array::f_style> &rows, const py::array_t<int, py::array::f_style> &cols, underlying_type<CoefficientPrecision> epsilon) const = 0; // LCOV_EXCL_LINE
+    virtual bool build_low_rank_approximation(const py::array_t<int, py::array::f_style> &rows, const py::array_t<int, py::array::f_style> &cols, underlying_type<CoefficientPrecision> epsilon) const = 0; // LCOV_EXCL_LINE
 
     void set_U(py::array_t<CoefficientPrecision, py::array::f_style> U0) {
         m_mats_U.push_back(U0); // no copy here
@@ -47,9 +57,9 @@ class PyVirtualLowRankGenerator : public VirtualLowRankGeneratorPython<Coefficie
     using VirtualLowRankGeneratorPython<CoefficientPrecision>::VirtualLowRankGeneratorPython;
 
     /* Trampoline (need one for each virtual function) */
-    virtual void build_low_rank_approximation(const py::array_t<int, py::array::f_style> &rows, const py::array_t<int, py::array::f_style> &cols, underlying_type<CoefficientPrecision> epsilon) const override {
+    virtual bool build_low_rank_approximation(const py::array_t<int, py::array::f_style> &rows, const py::array_t<int, py::array::f_style> &cols, underlying_type<CoefficientPrecision> epsilon) const override {
         PYBIND11_OVERRIDE_PURE(
-            void,                                                /* Return type */
+            bool,                                                /* Return type */
             VirtualLowRankGeneratorPython<CoefficientPrecision>, /* Parent class */
             build_low_rank_approximation,                        /* Name of function in C++ (must match Python name) */
             rows,