add off diagonal approximation

.github/workflows/CI.yml

@@ -29,27 +29,27 @@ jobs:
           - name: "Ubuntu OpenMPI g++"
             CC_COMPILER: gcc
             CXX_COMPILER: g++
-            container: "ubuntu_openmpi"
+            container: "ubuntu_gcc_openmpi"
             MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe"
             USE_SANITIZER: ""
             CMAKE_BUILD_TYPE: Release
-            DOCKER_OPTIONS: " "
+            DOCKER_OPTIONS: "--cap-add SYS_PTRACE"
             CODE_COVERAGE: "OFF"
 
           - name: "Ubuntu OpenMPI clang++"
             CC_COMPILER: clang
             CXX_COMPILER: clang++
-            container: "ubuntu_openmpi"
+            container: "ubuntu_clang_openmpi"
             MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe"
             USE_SANITIZER: ""
             CMAKE_BUILD_TYPE: Release
-            DOCKER_OPTIONS: " "
+            DOCKER_OPTIONS: "--cap-add SYS_PTRACE"
             CODE_COVERAGE: "OFF"
 
           - name: "Ubuntu MPICH g++"
             CC_COMPILER: gcc
             CXX_COMPILER: g++
-            container: "ubuntu_mpich"
+            container: "ubuntu_gcc_mpich"
             MPIEXEC_PREFLAGS: ""
             CMAKE_BUILD_TYPE: Debug
             DOCKER_OPTIONS: " "
@@ -58,35 +58,37 @@ jobs:
           - name: "Ubuntu MPICH clang++"
             CC_COMPILER: clang
             CXX_COMPILER: clang++
-            container: "ubuntu_mpich"
+            container: "ubuntu_clang_mpich"
             MPIEXEC_PREFLAGS: ""
             CMAKE_BUILD_TYPE: Release
             DOCKER_OPTIONS: " "
             CODE_COVERAGE: "OFF"
 
           # Hangs on github
-          # - name: "Debian OpenMPI g++"
-          #   CC_COMPILER: gcc
-          #   CXX_COMPILER: g++
-          #   container: "debian_openmpi"
-          #   MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe --mca btl_vader_single_copy_mechanism none"
-          #   USE_SANITIZER: ""
-          #   CMAKE_BUILD_TYPE: Debug
-          #   DOCKER_OPTIONS: "--cap-add SYS_PTRACE"
-
-          # - name: "Debian OpenMPI clang++"
-          #   CC_COMPILER: clang
-          #   CXX_COMPILER: clang++
-          #   container: "debian_openmpi"
-          #   MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe --mca btl_vader_single_copy_mechanism none"
-          #   USE_SANITIZER: ""
-          #   CMAKE_BUILD_TYPE: Debug
-          #   DOCKER_OPTIONS: "--cap-add SYS_PTRACE"
+          - name: "Debian OpenMPI g++"
+            CC_COMPILER: gcc
+            CXX_COMPILER: g++
+            container: "debian_gcc_openmpi"
+            MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe"
+            USE_SANITIZER: ""
+            CMAKE_BUILD_TYPE: Release
+            DOCKER_OPTIONS: "--cap-add SYS_PTRACE"
+            CODE_COVERAGE: "OFF"
+
+          - name: "Debian OpenMPI clang++"
+            CC_COMPILER: clang
+            CXX_COMPILER: clang++
+            container: "debian_clang_openmpi"
+            MPIEXEC_PREFLAGS: "--allow-run-as-root --oversubscribe"
+            USE_SANITIZER: ""
+            CMAKE_BUILD_TYPE: Release
+            DOCKER_OPTIONS: "--cap-add SYS_PTRACE"
+            CODE_COVERAGE: "OFF"
 
           - name: "Debian MPICH g++"
             CC_COMPILER: gcc
             CXX_COMPILER: g++
-            container: "debian_mpich"
+            container: "debian_gcc_mpich"
             MPIEXEC_PREFLAGS: ""
             CMAKE_BUILD_TYPE: Release
             DOCKER_OPTIONS: " "
@@ -95,7 +97,7 @@ jobs:
           - name: "Debian MPICH clang++"
             CC_COMPILER: clang
             CXX_COMPILER: clang++
-            container: "debian_mpich"
+            container: "debian_clang_mpich"
             MPIEXEC_PREFLAGS: ""
             CMAKE_BUILD_TYPE: Release
             DOCKER_OPTIONS: " "
@@ -115,7 +117,7 @@ jobs:
     if: "!contains(github.event.head_commit.message, '[ci skip]')"
     steps:
       - name: Checkout htool-python
-        uses: actions/checkout@v2
+        uses: actions/checkout@v3
         with:
           submodules: "true"
 
@@ -141,7 +143,9 @@ jobs:
 
       - name: Check format
         if: matrix.CODE_COVERAGE == 'ON'
+        # https://github.com/actions/checkout/pull/762
         run: |
+          git config --global --add safe.directory `pwd` 
           cd build
           make format
           make cmake-format

src/htool/cluster.hpp

@@ -16,7 +16,6 @@ using namespace htool;
 
 template <typename ClusterType>
 void declare_Cluster(py::module &m, const std::string &className) {
-
     py::class_<ClusterType, std::shared_ptr<ClusterType>, VirtualCluster> py_class(m, className.c_str());
     py_class.def(py::init<int>());
     py_class.def(
@@ -55,6 +54,7 @@ void declare_Cluster(py::module &m, const std::string &className) {
         }
         return py::array_t<int, py::array::f_style>(shape, masteroffset_raw.data());
     });
+    py_class.def("get_global_perm", overload_cast_<int>()(&ClusterType::get_global_perm, py::const_));
     py_class.def(
         "display", [](ClusterType &self, py::array_t<double, py::array::f_style | py::array::forcecast> x, int depth, bool show, MPI_Comm_wrapper comm) {
             int rankWorld;
@@ -74,7 +74,7 @@ void declare_Cluster(py::module &m, const std::string &className) {
                 // Permuted geometric points
                 for (int i = 0; i < size; ++i) {
                     for (int p = 0; p < space_dim; p++) {
-                        output[i + size * p] = x.at(p, root->get_perm(i));
+                        output[i + size * p] = x.at(p, root->get_global_perm(i));
                     }
                 }
 
@@ -143,4 +143,4 @@ void declare_Cluster(py::module &m, const std::string &className) {
     py_class.def("set_minclustersize", &ClusterType::set_minclustersize);
 }
 
-#endif
+#endif

src/htool/ddm_solver.hpp

@@ -98,4 +98,4 @@ void declare_DDM(py::module &m, const std::string &className) {
     py_class.def("get_infos", &Class::get_infos);
 }
 
-#endif
+#endif

src/htool/hmatrix.hpp

@@ -9,6 +9,7 @@
 #include "lrmat_generator.hpp"
 #include "matrix.hpp"
 #include "misc.hpp"
+#include "off_diagonal_approximation.hpp"
 #include "wrapper_mpi.hpp"
 #include <htool/htool.hpp>
 
@@ -48,6 +49,9 @@ void declare_HMatrix(py::module &m, const std::string &baseclassName, const std:
     py_class.def("set_compression", [](Class &self, std::shared_ptr<VirtualLowRankGeneratorCpp<T>> mat) {
         self.set_compression(mat);
     });
+    py_class.def("set_off_diagonal_approximation", [](Class &self, std::shared_ptr<VirtualOffDiagonalApproximation<T>> mat) {
+        self.set_off_diagonal_approximation(mat);
+    });
 
     // Getters
     py_class.def_property_readonly("shape", [](const Class &self) {
@@ -57,6 +61,14 @@ void declare_HMatrix(py::module &m, const std::string &baseclassName, const std:
     py_class.def("get_perm_s", overload_cast_<>()(&Class::get_perms, py::const_));
     py_class.def("get_MasterOffset_t", overload_cast_<>()(&Class::get_MasterOffset_t, py::const_));
     py_class.def("get_MasterOffset_s", overload_cast_<>()(&Class::get_MasterOffset_s, py::const_));
+    py_class.def("get_off_diagonal_geometries", [](const Class &self, const py::array_t<double, py::array::f_style> &xt, const py::array_t<double, py::array::f_style> &xs) {
+        int new_nc, new_nr, spatial_dimension;
+        self.get_off_diagonal_size(new_nr, new_nc);
+        spatial_dimension = self.get_target_cluster()->get_space_dim();
+        py::array_t<double, py::array::f_style> new_xs(std::array<long int, 2>{spatial_dimension, new_nc}), new_xt(std::array<long int, 2>{spatial_dimension, new_nr});
+        self.get_off_diagonal_geometries(xt.data(), xs.data(), new_xt.mutable_data(), new_xs.mutable_data());
+        return std::tuple<py::array_t<double, py::array::f_style>, py::array_t<double, py::array::f_style>>(new_xt, new_xs);
+    });
 
     // Linear algebra
     py_class.def("__mul__", [](const Class &self, std::vector<T> b) {
@@ -102,50 +114,9 @@ void declare_HMatrix(py::module &m, const std::string &baseclassName, const std:
     // Plot pattern
     py_class.def(
         "display", [](const Class &self, bool show = true) {
-            const std::vector<LowRankMatrix<T> *> &lrmats = self.get_MyFarFieldMats();
-            const std::vector<SubMatrix<T> *> &dmats      = self.get_MyNearFieldMats();
-
-            int nb        = dmats.size() + lrmats.size();
-            int sizeworld = self.get_sizeworld();
-            int rankworld = self.get_rankworld();
-
-            int nbworld[sizeworld];
-            MPI_Allgather(&nb, 1, MPI_INT, nbworld, 1, MPI_INT, self.get_comm());
-            int nbg = 0;
-            for (int i = 0; i < sizeworld; i++) {
-                nbg += nbworld[i];
-            }
-
-            std::vector<int> buf(5 * nbg, 0);
-
-            for (int i = 0; i < dmats.size(); i++) {
-                const SubMatrix<T> &l = *(dmats[i]);
-                buf[5 * i]            = l.get_offset_i();
-                buf[5 * i + 1]        = l.nb_rows();
-                buf[5 * i + 2]        = l.get_offset_j();
-                buf[5 * i + 3]        = l.nb_cols();
-                buf[5 * i + 4]        = -1;
-            }
-
-            for (int i = 0; i < lrmats.size(); i++) {
-                const LowRankMatrix<T> &l       = *(lrmats[i]);
-                buf[5 * (dmats.size() + i)]     = l.get_offset_i();
-                buf[5 * (dmats.size() + i) + 1] = l.nb_rows();
-                buf[5 * (dmats.size() + i) + 2] = l.get_offset_j();
-                buf[5 * (dmats.size() + i) + 3] = l.nb_cols();
-                buf[5 * (dmats.size() + i) + 4] = l.rank_of();
-            }
-
-            int displs[sizeworld];
-            int recvcounts[sizeworld];
-            displs[0] = 0;
-
-            for (int i = 0; i < sizeworld; i++) {
-                recvcounts[i] = 5 * nbworld[i];
-                if (i > 0)
-                    displs[i] = displs[i - 1] + recvcounts[i - 1];
-            }
-            MPI_Gatherv(rankworld == 0 ? MPI_IN_PLACE : buf.data(), recvcounts[rankworld], MPI_INT, buf.data(), recvcounts, displs, MPI_INT, 0, self.get_comm());
+            int rankworld        = self.get_rankworld();
+            std::vector<int> buf = self.get_output();
+            int nbg              = buf.size() / 5;
 
             if (rankworld == 0) {
                 // Import
@@ -158,6 +129,7 @@ void declare_HMatrix(py::module &m, const std::string &baseclassName, const std:
                 int nr = self.nb_rows();
                 int nc = self.nb_cols();
                 py::array_t<int> matrix({nr, nc});
+                matrix.attr("fill")(0);
                 py::array_t<bool> mask_matrix({nr, nc});
                 mask_matrix.attr("fill")(false);
 
@@ -246,7 +218,7 @@ void declare_HMatrix(py::module &m, const std::string &baseclassName, const std:
                 // Permuted geometric points
                 for (int i = 0; i < size; ++i) {
                     for (int p = 0; p < space_dim; p++) {
-                        output[i + size * p] = points_target.at(p, root->get_perm(i));
+                        output[i + size * p] = points_target.at(p, root->get_global_perm(i));
                     }
                 }
 

src/htool/lrmat_generator.hpp

@@ -15,27 +15,28 @@ using namespace htool;
 template <typename T>
 class VirtualLowRankGeneratorCpp : public VirtualLowRankGenerator<T> {
 
-    py::array_t<T, py::array::f_style> mat_U, mat_V;
     int rank;
+    mutable std::vector<py::array_t<T, py::array::f_style>> mats_U; // owned by Python
+    mutable std::vector<py::array_t<T, py::array::f_style>> mats_V; // owned by Python
 
   public:
     using VirtualLowRankGenerator<T>::VirtualLowRankGenerator;
 
     void copy_low_rank_approximation(double epsilon, int M, int N, const int *const rows, const int *const cols, int &rank0, T **U, T **V, const VirtualGenerator<T> &A, const VirtualCluster &t, const double *const xt, const VirtualCluster &s, const double *const xs) const override {
 
         build_low_rank_approximation(epsilon, rank0, A, std::vector<int>(rows, rows + M), std::vector<int>(cols, cols + N));
-        *U    = new T[M * rank];
-        *V    = new T[N * rank];
+        *U    = mats_U.back().mutable_data();
+        *V    = mats_V.back().mutable_data();
         rank0 = rank;
-        std::copy_n(mat_U.data(), mat_U.size(), *U);
-        std::copy_n(mat_V.data(), mat_V.size(), *V);
     }
 
+    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(double epsilon, int rank, const VirtualGenerator<T> &A, const std::vector<int> &J, const std::vector<int> &K) const = 0; // LCOV_EXCL_LINE
 
-    void set_U(py::array_t<T, py::array::f_style> U0) { mat_U = U0; }
-    void set_V(py::array_t<T, py::array::f_style> V0) { mat_V = V0; }
+    void set_U(py::array_t<T, py::array::f_style> U0) { mats_U.push_back(U0); }
+    void set_V(py::array_t<T, py::array::f_style> V0) { mats_V.push_back(V0); }
     void set_rank(int rank0) { rank = rank0; }
 };
 

src/htool/main.cpp

@@ -4,6 +4,7 @@
 #include "hmatrix.hpp"
 #include "lrmat_generator.hpp"
 #include "matrix.hpp"
+#include "off_diagonal_approximation.hpp"
 #include "wrapper_mpi.hpp"
 
 PYBIND11_MODULE(Htool, m) {
@@ -27,6 +28,9 @@ PYBIND11_MODULE(Htool, m) {
     declare_HMatrix<double>(m, "HMatrixVirtual", "HMatrix");
     declare_HMatrix<std::complex<double>>(m, "ComplexHMatrixVirtual", "ComplexHMatrix");
 
+    declare_VirtualOffDiagonalApproximation<double>(m, "VirtualOffDiagonalApproximation", "CustomOffDiagonalApproximation", "HMatrixOffDiagonalApproximation");
+    declare_VirtualOffDiagonalApproximation<std::complex<double>>(m, "ComplexVirtualOffDiagonalApproximation", "ComplexCustomOffDiagonalApproximation", "ComplexHMatrixOffDiagonalApproximation");
+
     declare_DDM<double>(m, "DDM");
     declare_DDM<std::complex<double>>(m, "ComplexDDM");
 }