improve coverage

Author: PierreMarchand20

lib/htool

@@ -2,20 +2,25 @@
 #define HTOOL_PYTHON_GENEO_COARSE_SPACE_BUILDER_HPP
 
 #include <htool/hmatrix/linalg/add_hmatrix_vector_product.hpp>
-#include <htool/solvers/geneo/coarse_space_builder.hpp>
+#include <htool/solvers/interfaces/virtual_coarse_space_builder.hpp>
 #include <pybind11/pybind11.h>
 namespace py = pybind11;
 
-template <typename CoefficientPrecision, typename CoordinatePrecision = CoefficientPrecision>
-class VirtualGeneoCoarseSpaceBuilderPython : public VirtualGeneoCoarseSpaceBuilder<CoefficientPrecision, CoordinatePrecision> {
+template <typename CoefficientPrecision, typename CoordinatePrecision = htool::underlying_type<CoefficientPrecision>>
+class VirtualGeneoCoarseSpaceBuilderPython : public VirtualCoarseSpaceBuilder<CoefficientPrecision> {
     py::array_t<CoefficientPrecision, py::array::f_style> m_coarse_space;
 
+    int m_size_wo_overlap;
+    int m_size_with_overlap;
+    const HMatrix<CoefficientPrecision, CoordinatePrecision> &m_local_hmatrix;
+    int m_geneo_nu                                                 = 2;
+    htool::underlying_type<CoefficientPrecision> m_geneo_threshold = -1.;
+
   public:
-    char get_symmetry() { return this->m_symmetry; }
     int get_geneo_nu() { return this->m_geneo_nu; }
     htool::underlying_type<CoefficientPrecision> get_geneo_threshold() { return this->m_geneo_threshold; }
 
-    explicit VirtualGeneoCoarseSpaceBuilderPython(int size_wo_overlap, int size_with_overlap, const HMatrix<CoefficientPrecision> &Ai, char symmetry, char uplo, int geneo_nu, htool::underlying_type<CoefficientPrecision> geneo_threshold) : VirtualGeneoCoarseSpaceBuilder<CoefficientPrecision>(size_wo_overlap, size_with_overlap, Ai, symmetry, uplo, geneo_nu, geneo_threshold) {}
+    explicit VirtualGeneoCoarseSpaceBuilderPython(int size_wo_overlap, int size_with_overlap, const HMatrix<CoefficientPrecision> &Ai, int geneo_nu, htool::underlying_type<CoefficientPrecision> geneo_threshold) : m_size_wo_overlap(size_wo_overlap), m_size_with_overlap(size_with_overlap), m_local_hmatrix(Ai), m_geneo_nu(geneo_nu), m_geneo_threshold(geneo_threshold) {}
 
     Matrix<CoefficientPrecision> build_coarse_space() override {
         std::function<py::array_t<CoefficientPrecision>(py::array_t<CoefficientPrecision>)> hmatrix_callback;
@@ -74,28 +79,23 @@ void declare_virtual_geneo_coarse_space_builder(py::module &m, const std::string
 
     using Class = VirtualGeneoCoarseSpaceBuilderPython<CoefficientPrecision>;
     py::class_<Class, PyVirtualGeneoCoarseSpaceBuilder<CoefficientPrecision>, VirtualCoarseSpaceBuilder<CoefficientPrecision>> py_class(m, className.c_str());
-    py_class.def(py::init([](int size_wo_overlap, int size_with_overlap, const HMatrix<CoefficientPrecision> &Ai, char symmetry, char uplo, int geneo_nu) {
-                     return PyVirtualGeneoCoarseSpaceBuilder<CoefficientPrecision>(size_wo_overlap, size_with_overlap, Ai, symmetry, uplo, geneo_nu, -1);
+    py_class.def(py::init([](int size_wo_overlap, int size_with_overlap, const HMatrix<CoefficientPrecision> &Ai, int geneo_nu) {
+                     return PyVirtualGeneoCoarseSpaceBuilder<CoefficientPrecision>(size_wo_overlap, size_with_overlap, Ai, geneo_nu, -1);
                  }),
                  py::arg("size_wo_overlap"),
                  py::arg("size_with_overlap"),
                  py::arg("Ai"),
-                 py::arg("symmetry"),
-                 py::arg("uplo"),
                  py::kw_only(),
                  py::arg("geneo_nu"));
-    py_class.def(py::init([](int size_wo_overlap, int size_with_overlap, const HMatrix<CoefficientPrecision> &Ai, char symmetry, char uplo, double geneo_threshold) {
-                     return PyVirtualGeneoCoarseSpaceBuilder<CoefficientPrecision>(size_wo_overlap, size_with_overlap, Ai, symmetry, uplo, 0, geneo_threshold);
+    py_class.def(py::init([](int size_wo_overlap, int size_with_overlap, const HMatrix<CoefficientPrecision> &Ai, double geneo_threshold) {
+                     return PyVirtualGeneoCoarseSpaceBuilder<CoefficientPrecision>(size_wo_overlap, size_with_overlap, Ai, 0, geneo_threshold);
                  }),
                  py::arg("size_wo_overlap"),
                  py::arg("size_with_overlap"),
                  py::arg("Ai"),
-                 py::arg("symmetry"),
-                 py::arg("uplo"),
                  py::kw_only(),
                  py::arg("geneo_threshold"));
     py_class.def("set_coarse_space", &Class::set_coarse_space);
-    py_class.def_property_readonly("symmetry", &Class::get_symmetry);
     py_class.def_property_readonly("geneo_nu", &Class::get_geneo_nu);
     py_class.def_property_readonly("geneo_threshold", &Class::get_geneo_threshold);
 }

src/htool/solver/geneo/coarse_space_builder.hpp

@@ -6,10 +6,13 @@
 import numpy as np
 import pytest
 from conftest import ComplexGeneratorFromMatrix, GeneratorFromMatrix
-from scipy.linalg import eig, eigh
+from scipy.linalg import eig, eigh, ldl, solve_triangular
+from scipy.sparse.linalg import LinearOperator, eigsh
 
 
-class CustomGeneoCoarseSpaceDenseBuilder(Htool.VirtualGeneoCoarseSpaceDenseBuilder):
+class CustomDenseGeneoCoarseSpaceDenseBuilder(
+    Htool.VirtualGeneoCoarseSpaceDenseBuilder
+):
     def compute_coarse_space(self, Ai, Bi):
         coarse_space = None
         if self.symmetry == "S" or self.symmetry == "H":
@@ -30,6 +33,102 @@ def compute_coarse_space(self, Ai, Bi):
         self.set_coarse_space(coarse_space)
 
 
+class CustomGeneoCoarseSpaceBuilder(Htool.VirtualGeneoCoarseSpaceBuilder):
+    class DirichletHMatrix(LinearOperator):
+        def __init__(self, hmatrix_callback, size, dtype):
+            self.hmatrix_callback = hmatrix_callback
+            self.shape = (size, size)
+            self.dtype = dtype
+
+        def _matvec(self, x):
+            return self.hmatrix_callback(x)
+
+    class DenseNeumannMatrix(LinearOperator):
+        def __init__(self, Bi):
+            self.L, self.D, self.perm = ldl(Bi)
+            self.shape = Bi.shape
+            self.dtype = Bi.dtype
+
+        def _matvec(self, x):
+            tmp = solve_triangular(self.L[self.perm, :], x[self.perm], lower=True)
+            for index, coef in enumerate(self.D.diagonal()):
+                #     if coef>1e-4:
+                tmp[index] /= coef
+            #     else :
+            #         tmp[index]=0
+            res = np.empty_like(x)
+            res[self.perm] = solve_triangular(
+                self.L[self.perm, :], tmp, lower=True, trans="T"
+            )
+            return res
+
+    def __init__(
+        self,
+        size_wo_overlap,
+        size_with_overlap,
+        Ai,
+        Bi,
+        geneo_nu=2,
+        geneo_threshold=None,
+    ):
+        if not geneo_threshold:
+            Htool.VirtualGeneoCoarseSpaceBuilder.__init__(
+                self,
+                size_wo_overlap,
+                size_with_overlap,
+                Ai,
+                geneo_nu=geneo_nu,
+            )
+        else:
+            Htool.VirtualGeneoCoarseSpaceBuilder.__init__(
+                self,
+                size_wo_overlap,
+                size_with_overlap,
+                Ai,
+                geneo_threshold=geneo_threshold,
+            )
+        self.local_size = size_with_overlap
+        self.Bi = Bi
+        self.Ai = Ai
+
+    def compute_coarse_space(self, hmatrix_callback):
+        self.coarse_space = None
+        self.eigenvalues = None
+
+        if self.geneo_threshold > 0:
+            [w, v] = eigsh(
+                self.Bi,
+                10,
+                M=self.DirichletHMatrix(
+                    hmatrix_callback, self.local_size, np.dtype("float64")
+                ),
+                OPinv=self.DenseNeumannMatrix(self.Bi),
+                return_eigenvectors=True,
+                sigma=0,
+            )
+        else:
+            [w, v] = eigsh(
+                self.Bi,
+                self.geneo_nu,
+                M=self.DirichletHMatrix(
+                    hmatrix_callback, self.local_size, np.dtype("float64")
+                ),
+                OPinv=self.DenseNeumannMatrix(self.Bi),
+                return_eigenvectors=True,
+                sigma=0,
+            )
+        idx = w.argsort()
+        if self.geneo_threshold > 0:
+            nb_eig = (w < 1.0 / self.geneo_threshold).sum()
+            self.coarse_space = v[:, idx[0:nb_eig]].real
+            self.eigenvalues = 1.0 / w[idx[:]]
+        else:
+            self.coarse_space = v[:, idx[0 : self.geneo_nu]]
+            self.eigenvalues = 1.0 / w[idx[:]]
+
+        self.set_coarse_space(self.coarse_space)
+
+
 @pytest.mark.parametrize("epsilon", [1e-6])
 @pytest.mark.parametrize("eta", [10])
 @pytest.mark.parametrize("tol", [1e-6])
@@ -96,6 +195,70 @@ def compute_coarse_space(self, Ai, Bi):
             "additive",
             "geneo_threshold",
         ),
+        (
+            1,
+            "S",
+            "DDMWithHMatrixPlusOverlap",
+            "asm",
+            "additive",
+            "custom_dense_geneo_nu",
+        ),
+        (
+            1,
+            "S",
+            "DDMWithHMatrixPlusOverlap",
+            "ras",
+            "additive",
+            "custom_dense_geneo_nu",
+        ),
+        (
+            10,
+            "S",
+            "DDMWithHMatrixPlusOverlap",
+            "asm",
+            "additive",
+            "custom_dense_geneo_nu",
+        ),
+        (
+            10,
+            "S",
+            "DDMWithHMatrixPlusOverlap",
+            "ras",
+            "additive",
+            "custom_dense_geneo_nu",
+        ),
+        (
+            1,
+            "S",
+            "DDMWithHMatrixPlusOverlap",
+            "asm",
+            "additive",
+            "custom_dense_geneo_threshold",
+        ),
+        (
+            1,
+            "S",
+            "DDMWithHMatrixPlusOverlap",
+            "ras",
+            "additive",
+            "custom_dense_geneo_threshold",
+        ),
+        (
+            10,
+            "S",
+            "DDMWithHMatrixPlusOverlap",
+            "asm",
+            "additive",
+            "custom_dense_geneo_threshold",
+        ),
+        (
+            10,
+            "S",
+            "DDMWithHMatrixPlusOverlap",
+            "ras",
+            "additive",
+            "custom_dense_geneo_threshold",
+        ),
         (1, "S", "DDMWithHMatrixPlusOverlap", "asm", "additive", "custom_geneo_nu"),
         (1, "S", "DDMWithHMatrixPlusOverlap", "ras", "additive", "custom_geneo_nu"),
         (
@@ -156,6 +319,14 @@ def compute_coarse_space(self, Ai, Bi):
         (1, "S", "DDMWithHMatrix", "ras", "additive", "geneo_threshold"),
         (10, "S", "DDMWithHMatrix", "asm", "additive", "geneo_threshold"),
         (10, "S", "DDMWithHMatrix", "ras", "additive", "geneo_threshold"),
+        (1, "S", "DDMWithHMatrix", "asm", "additive", "custom_dense_geneo_nu"),
+        (1, "S", "DDMWithHMatrix", "ras", "additive", "custom_dense_geneo_nu"),
+        (10, "S", "DDMWithHMatrix", "asm", "additive", "custom_dense_geneo_nu"),
+        (10, "S", "DDMWithHMatrix", "ras", "additive", "custom_dense_geneo_nu"),
+        (1, "S", "DDMWithHMatrix", "asm", "additive", "custom_dense_geneo_threshold"),
+        (1, "S", "DDMWithHMatrix", "ras", "additive", "custom_dense_geneo_threshold"),
+        (10, "S", "DDMWithHMatrix", "asm", "additive", "custom_dense_geneo_threshold"),
+        (10, "S", "DDMWithHMatrix", "ras", "additive", "custom_dense_geneo_threshold"),
         (1, "S", "DDMWithHMatrix", "asm", "additive", "custom_geneo_nu"),
         (1, "S", "DDMWithHMatrix", "ras", "additive", "custom_geneo_nu"),
         (10, "S", "DDMWithHMatrix", "asm", "additive", "custom_geneo_nu"),
@@ -383,7 +554,24 @@ def test_ddm_solver(
                 geneo_threshold=100,
             )
         elif geneo_type == "custom_geneo_nu":
-            geneo_space_operator_builder = CustomGeneoCoarseSpaceDenseBuilder(
+            print("ICI")
+            geneo_space_operator_builder = CustomGeneoCoarseSpaceBuilder(
+                local_size_wo_overlap,
+                local_size_with_overlap,
+                default_approximation.block_diagonal_hmatrix,
+                local_neumann_matrix,
+                geneo_nu=2,
+            )
+        elif geneo_type == "custom_geneo_threshold":
+            geneo_space_operator_builder = CustomGeneoCoarseSpaceBuilder(
+                local_size_wo_overlap,
+                local_size_with_overlap,
+                default_approximation.block_diagonal_hmatrix,
+                local_neumann_matrix,
+                geneo_threshold=100,
+            )
+        elif geneo_type == "custom_dense_geneo_nu":
+            geneo_space_operator_builder = CustomDenseGeneoCoarseSpaceDenseBuilder(
                 local_size_wo_overlap,
                 local_size_with_overlap,
                 default_approximation.block_diagonal_hmatrix.to_dense(),
@@ -392,8 +580,8 @@ def test_ddm_solver(
                 UPLO,
                 geneo_nu=2,
             )
-        elif geneo_type == "custom_geneo_threshold":
-            geneo_space_operator_builder = CustomGeneoCoarseSpaceDenseBuilder(
+        elif geneo_type == "custom_dense_geneo_threshold":
+            geneo_space_operator_builder = CustomDenseGeneoCoarseSpaceDenseBuilder(
                 local_size_wo_overlap,
                 local_size_with_overlap,
                 default_approximation.block_diagonal_hmatrix.to_dense(),