diff --git a/example/define_custom_generators.py b/example/define_custom_generators.py
index 2519eab..0703093 100644
--- a/example/define_custom_generators.py
+++ b/example/define_custom_generators.py
@@ -2,7 +2,7 @@
 import numpy as np
 
 
-class CustomGenerator(Htool.VirtualGeneratorInUserNumbering):
+class CustomGenerator(Htool.VirtualGenerator):
     def __init__(self, target_points, source_points):
         super().__init__()
         self.target_points = target_points
diff --git a/example/define_custom_local_operator.py b/example/define_custom_local_operator.py
index 526a38b..e2bdb28 100644
--- a/example/define_custom_local_operator.py
+++ b/example/define_custom_local_operator.py
@@ -5,7 +5,7 @@
 class CustomLocalOperator(Htool.LocalOperator):
     def __init__(
         self,
-        generator: Htool.VirtualGeneratorInUserNumbering,
+        generator: Htool.VirtualGenerator,
         target_cluster: Htool.Cluster,
         source_cluster: Htool.Cluster,
         symmetry: str = "N",
diff --git a/example/define_custom_low_rank_generator.py b/example/define_custom_low_rank_generator.py
index dab4378..f4f9123 100644
--- a/example/define_custom_low_rank_generator.py
+++ b/example/define_custom_low_rank_generator.py
@@ -5,12 +5,13 @@
 
 
 class CustomSVD(Htool.VirtualLowRankGenerator):
-    def __init__(self, generator: Htool.VirtualGeneratorInUserNumbering):
-        super().__init__(generator)
+    def __init__(self, generator: Htool.VirtualGenerator):
+        super().__init__()
+        self.generator = generator
 
     def build_low_rank_approximation(self, rows, cols, epsilon):
         submat = np.zeros((len(rows), len(cols)), order="F")
-        self.build_submatrix(rows, cols, submat)
+        self.generator.build_submatrix(rows, cols, submat)
         u, s, vh = np.linalg.svd(submat, full_matrices=False)
 
         norm = np.linalg.norm(submat)
@@ -22,4 +23,3 @@ def build_low_rank_approximation(self, rows, cols, epsilon):
         count = min(count + 1, min(len(rows), len(cols)))
         self.set_U(u[:, 0:count] * s[0:count])
         self.set_V(vh[0:count, :])
-        self.set_rank(count)
diff --git a/example/use_block_jacobi_solver.py b/example/use_block_jacobi_solver.py
index d3301a9..eb9c277 100644
--- a/example/use_block_jacobi_solver.py
+++ b/example/use_block_jacobi_solver.py
@@ -44,7 +44,6 @@
 
 # Solver with block Jacobi preconditionner
 block_diagonal_hmatrix = copy.deepcopy(default_approximation.block_diagonal_hmatrix)
-test = default_approximation.block_diagonal_hmatrix
 default_solver_builder = Htool.DDMSolverBuilder(
     default_approximation.distributed_operator, block_diagonal_hmatrix
 )
diff --git a/lib/htool b/lib/htool
index 29ce0df..1a3b198 160000
--- a/lib/htool
+++ b/lib/htool
@@ -1 +1 @@
-Subproject commit 29ce0dfe8125b6ebbd6efc1dfff6229dbaf92f38
+Subproject commit 1a3b198ffc6f73cd62059094ca7b606d151da976
diff --git a/src/htool/distributed_operator/utility.hpp b/src/htool/distributed_operator/utility.hpp
index 8b5f292..3eceecb 100644
--- a/src/htool/distributed_operator/utility.hpp
+++ b/src/htool/distributed_operator/utility.hpp
@@ -24,7 +24,7 @@ void declare_distributed_operator_utility(py::module &m, std::string prefix = ""
         "distributed_operator", [](const CustomApproximation &self) { return &self.distributed_operator; }, py::return_value_policy::reference_internal);
 
     py::class_<DefaultApproximation> default_approximation_class(m, default_approximation_name.c_str());
-    default_approximation_class.def(py::init<const VirtualGeneratorInUserNumbering<CoefficientPrecision> &, const Cluster<CoordinatePrecision> &, const Cluster<CoordinatePrecision> &, htool::underlying_type<CoefficientPrecision>, htool::underlying_type<CoefficientPrecision>, char, char, MPI_Comm_wrapper>());
+    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(
@@ -33,7 +33,7 @@ void declare_distributed_operator_utility(py::module &m, std::string prefix = ""
         "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 VirtualGeneratorInUserNumbering<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(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(
@@ -42,7 +42,7 @@ void declare_distributed_operator_utility(py::module &m, std::string prefix = ""
         "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 VirtualGeneratorInUserNumbering<CoefficientPrecision> &, const Cluster<CoordinatePrecision> &, const Cluster<CoordinatePrecision> &, const HMatrixTreeBuilder<CoefficientPrecision, CoordinatePrecision> &, MPI_Comm_wrapper>());
+    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(
diff --git a/src/htool/hmatrix/hmatrix_builder.hpp b/src/htool/hmatrix/hmatrix_builder.hpp
index 2a4578d..613c25c 100644
--- a/src/htool/hmatrix/hmatrix_builder.hpp
+++ b/src/htool/hmatrix/hmatrix_builder.hpp
@@ -16,7 +16,7 @@ void declare_hmatrix_builder(py::module &m, const std::string &className) {
 
     // Build
     // py_class.def("build", [](const Class &self, const VirtualGenerator<CoefficientPrecision> &generator) { return self.build(generator); });
-    py_class.def("build", [](const Class &self, const VirtualGeneratorInUserNumbering<CoefficientPrecision> &generator) { return self.build(GeneratorWithPermutation<CoefficientPrecision>(generator, self.get_target_cluster().get_permutation().data(), self.get_source_cluster().get_permutation().data())); });
+    py_class.def("build", [](const Class &self, const VirtualGenerator<CoefficientPrecision> &generator) { return self.build(InternalGeneratorWithPermutation<CoefficientPrecision>(generator, self.get_target_cluster().get_permutation().data(), self.get_source_cluster().get_permutation().data())); });
 
     // Setters
     py_class.def("set_minimal_source_depth", &Class::set_minimal_source_depth);
diff --git a/src/htool/hmatrix/interfaces/virtual_generator.hpp b/src/htool/hmatrix/interfaces/virtual_generator.hpp
index 8e4f816..1c66764 100644
--- a/src/htool/hmatrix/interfaces/virtual_generator.hpp
+++ b/src/htool/hmatrix/interfaces/virtual_generator.hpp
@@ -7,11 +7,11 @@
 using namespace htool;
 
 template <typename CoefficientPrecision>
-class VirtualGeneratorInUserNumberingPython : public htool::VirtualGeneratorInUserNumbering<CoefficientPrecision> {
+class VirtualGeneratorPython : public htool::VirtualGenerator<CoefficientPrecision> {
   public:
-    using VirtualGeneratorInUserNumbering<CoefficientPrecision>::VirtualGeneratorInUserNumbering;
+    using VirtualGenerator<CoefficientPrecision>::VirtualGenerator;
 
-    VirtualGeneratorInUserNumberingPython(const py::array_t<int> &target_permutation, const py::array_t<int> &source_permutation) : VirtualGeneratorInUserNumbering<CoefficientPrecision>() {}
+    VirtualGeneratorPython(const py::array_t<int> &target_permutation, const py::array_t<int> &source_permutation) : VirtualGenerator<CoefficientPrecision>() {}
 
     void copy_submatrix(int M, int N, const int *const rows, const int *const cols, CoefficientPrecision *ptr) const override {
         if (M * N > 0) {
@@ -29,15 +29,15 @@ class VirtualGeneratorInUserNumberingPython : public htool::VirtualGeneratorInUs
 };
 
 template <typename CoefficientPrecision>
-class PyVirtualGeneratorInUserNumbering : public VirtualGeneratorInUserNumberingPython<CoefficientPrecision> {
+class PyVirtualGenerator : public VirtualGeneratorPython<CoefficientPrecision> {
   public:
-    using VirtualGeneratorInUserNumberingPython<CoefficientPrecision>::VirtualGeneratorInUserNumberingPython;
+    using VirtualGeneratorPython<CoefficientPrecision>::VirtualGeneratorPython;
 
     /* Trampoline (need one for each virtual function) */
     virtual void build_submatrix(const py::array_t<int> &J, const py::array_t<int> &K, py::array_t<CoefficientPrecision, py::array::f_style> &mat) const override {
         PYBIND11_OVERRIDE_PURE(
             void,                                                        /* Return type */
-            VirtualGeneratorInUserNumberingPython<CoefficientPrecision>, /* Parent class */
+            VirtualGeneratorPython<CoefficientPrecision>, /* Parent class */
             build_submatrix,                                             /* Name of function in C++ (must match Python name) */
             J,
             K,
@@ -48,11 +48,11 @@ class PyVirtualGeneratorInUserNumbering : public VirtualGeneratorInUserNumbering
 
 template <typename CoefficientPrecision>
 void declare_virtual_generator(py::module &m, const std::string &className, const std::string &base_class_name) {
-    using BaseClass = VirtualGeneratorInUserNumbering<CoefficientPrecision>;
+    using BaseClass = VirtualGenerator<CoefficientPrecision>;
     py::class_<BaseClass>(m, base_class_name.c_str());
 
-    using Class = VirtualGeneratorInUserNumberingPython<CoefficientPrecision>;
-    py::class_<Class, BaseClass, PyVirtualGeneratorInUserNumbering<CoefficientPrecision>> py_class(m, className.c_str());
+    using Class = VirtualGeneratorPython<CoefficientPrecision>;
+    py::class_<Class, BaseClass, PyVirtualGenerator<CoefficientPrecision>> py_class(m, className.c_str());
     py_class.def(py::init<>());
     py_class.def("build_submatrix", &Class::build_submatrix);
 }
diff --git a/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp b/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp
index dcfdc8c..732a7ed 100644
--- a/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp
+++ b/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp
@@ -13,22 +13,19 @@ using namespace htool;
 
 template <typename CoefficientPrecision, typename CoordinatePrecision = CoefficientPrecision>
 class VirtualLowRankGeneratorPython : public VirtualLowRankGenerator<CoefficientPrecision, CoordinatePrecision> {
-    int m_rank;
     mutable std::vector<py::array_t<CoefficientPrecision, py::array::f_style>> m_mats_U; // owned by Python
     mutable std::vector<py::array_t<CoefficientPrecision, py::array::f_style>> m_mats_V; // owned by Python
-    const VirtualGeneratorInUserNumbering<CoefficientPrecision> &m_generator_in_user_numbering;
 
   public:
     using VirtualLowRankGenerator<CoefficientPrecision, CoordinatePrecision>::VirtualLowRankGenerator;
 
-    VirtualLowRankGeneratorPython(const VirtualGeneratorInUserNumbering<CoefficientPrecision> &generator_in_user_numbering) : m_generator_in_user_numbering(generator_in_user_numbering) {}
+    VirtualLowRankGeneratorPython() {}
 
-    void copy_low_rank_approximation(const VirtualGenerator<CoefficientPrecision> &A, const Cluster<CoordinatePrecision> &target_cluster, const Cluster<CoordinatePrecision> &source_cluster, underlying_type<CoefficientPrecision> epsilon, int &rank, Matrix<CoefficientPrecision> &U, Matrix<CoefficientPrecision> &V) const override {
-        py::array_t<int, py::array::f_style> rows(target_cluster.get_size(), target_cluster.get_permutation().data() + target_cluster.get_offset(), py::capsule(target_cluster.get_permutation().data()));
-        py::array_t<int, py::array::f_style> cols(source_cluster.get_size(), source_cluster.get_permutation().data() + source_cluster.get_offset(), py::capsule(source_cluster.get_permutation().data()));
+    void copy_low_rank_approximation(const VirtualInternalGenerator<CoefficientPrecision> &A, const Cluster<CoordinatePrecision> &target_cluster, const Cluster<CoordinatePrecision> &source_cluster, underlying_type<CoefficientPrecision> epsilon, int &rank, Matrix<CoefficientPrecision> &U, Matrix<CoefficientPrecision> &V) const override {
+        py::array_t<int, py::array::f_style> rows(std::array<int,1>{target_cluster.get_size()}, target_cluster.get_permutation().data() + target_cluster.get_offset(), py::capsule(target_cluster.get_permutation().data()));
+        py::array_t<int, py::array::f_style> cols(std::array<int,1>{source_cluster.get_size()}, source_cluster.get_permutation().data() + source_cluster.get_offset(), py::capsule(source_cluster.get_permutation().data()));
 
         build_low_rank_approximation(rows, cols, epsilon);
-        rank = m_rank;
         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);
     }
@@ -42,11 +39,7 @@ class VirtualLowRankGeneratorPython : public VirtualLowRankGenerator<Coefficient
         m_mats_U.push_back(U0); // no copy here
     }
     void set_V(py::array_t<CoefficientPrecision, py::array::f_style> V0) { m_mats_V.push_back(V0); }
-    void set_rank(int rank) { m_rank = rank; }
 
-    void build_submatrix(const py::array_t<int> rows, const py::array_t<int> cols, py::array_t<CoefficientPrecision, py::array::f_style> &mat) const {
-        m_generator_in_user_numbering.copy_submatrix(rows.size(), cols.size(), rows.data(), cols.data(), mat.mutable_data());
-    }
 };
 
 template <typename CoefficientPrecision, typename CoordinatePrecision>
@@ -70,12 +63,10 @@ template <typename CoefficientPrecision, typename CoordinatePrecision>
 void declare_custom_VirtualLowRankGenerator(py::module &m, const std::string &className) {
     using Class = VirtualLowRankGeneratorPython<CoefficientPrecision, CoordinatePrecision>;
     py::class_<Class, std::shared_ptr<Class>, PyVirtualLowRankGenerator<CoefficientPrecision, CoordinatePrecision>> py_class(m, className.c_str());
-    py_class.def(py::init<const VirtualGeneratorInUserNumbering<CoefficientPrecision> &>());
+    py_class.def(py::init<>());
     py_class.def("build_low_rank_approximation", &Class::build_low_rank_approximation);
     py_class.def("set_U", &Class::set_U);
     py_class.def("set_V", &Class::set_V);
-    py_class.def("set_rank", &Class::set_rank);
-    py_class.def("build_submatrix", &Class::build_submatrix);
 }
 
 #endif
diff --git a/src/htool/main.cpp b/src/htool/main.cpp
index a847324..8b6e606 100644
--- a/src/htool/main.cpp
+++ b/src/htool/main.cpp
@@ -56,7 +56,7 @@ PYBIND11_MODULE(Htool, m) {
 
     declare_hmatrix_builder<double, double>(m, "HMatrixBuilder");
     declare_HMatrix<double, double>(m, "HMatrix");
-    declare_virtual_generator<double>(m, "VirtualGeneratorInUserNumbering", "IGeneratorInUserNumbering");
+    declare_virtual_generator<double>(m, "VirtualGenerator", "IGenerator");
     declare_custom_VirtualLowRankGenerator<double, double>(m, "VirtualLowRankGenerator");
     declare_custom_VirtualDenseBlocksGenerator<double>(m, "VirtualDenseBlocksGenerator");
 
@@ -77,7 +77,7 @@ PYBIND11_MODULE(Htool, m) {
 
     declare_hmatrix_builder<std::complex<double>, double>(m, "ComplexHMatrixBuilder");
     declare_HMatrix<std::complex<double>, double>(m, "ComplexHMatrix");
-    declare_virtual_generator<std::complex<double>>(m, "ComplexVirtualGeneratorInUserNumbering", "IComplexGeneratorInUserNumbering");
+    declare_virtual_generator<std::complex<double>>(m, "ComplexVirtualGenerator", "IComplexGenerator");
 
     declare_distributed_operator<std::complex<double>>(m, "ComplexDistributedOperator");
     declare_distributed_operator_utility<std::complex<double>, double>(m, "Complex");
diff --git a/src/htool/solver/utility.hpp b/src/htool/solver/utility.hpp
index 47636fa..5a61b59 100644
--- a/src/htool/solver/utility.hpp
+++ b/src/htool/solver/utility.hpp
@@ -13,12 +13,13 @@ void declare_solver_utility(py::module &m, std::string prefix = "") {
     py::class_<DDMSolverBuilder> ddm_solver_class(m, ddm_solver_name.c_str());
     ddm_solver_class.def(py::init<DistributedOperator<CoefficientPrecision> &, HMatrix<CoefficientPrecision, CoordinatePrecision> &>());
 
-    ddm_solver_class.def(py::init<DistributedOperator<CoefficientPrecision> &, HMatrix<CoefficientPrecision, CoordinatePrecision> &, const VirtualGeneratorInUserNumbering<CoefficientPrecision> &, const std::vector<int> &, const std::vector<int> &, const std::vector<int> &, const std::vector<std::vector<int>> &>());
-    ddm_solver_class.def(py::init([](DistributedOperator<CoefficientPrecision> &distributed_operator, const std::vector<int> &ovr_subdomain_to_global, const std::vector<int> &cluster_to_ovr_subdomain, const std::vector<int> &neighbors, const std::vector<std::vector<int>> &intersections, const VirtualGeneratorInUserNumbering<CoefficientPrecision> &generator, const py::array_t<CoordinatePrecision, py::array::f_style | py::array::forcecast> coordinates, underlying_type<CoefficientPrecision> epsilon, CoordinatePrecision eta) {
+    ddm_solver_class.def(py::init<DistributedOperator<CoefficientPrecision> &, HMatrix<CoefficientPrecision, CoordinatePrecision> &, const VirtualGenerator<CoefficientPrecision> &, const std::vector<int> &, const std::vector<int> &, const std::vector<int> &, const std::vector<std::vector<int>> &>(),py::keep_alive<1, 2>(),py::keep_alive<1, 3>());
+    ddm_solver_class.def(py::init([](DistributedOperator<CoefficientPrecision> &distributed_operator, const std::vector<int> &ovr_subdomain_to_global, const std::vector<int> &cluster_to_ovr_subdomain, const std::vector<int> &neighbors, const std::vector<std::vector<int>> &intersections, const VirtualGenerator<CoefficientPrecision> &generator, const py::array_t<CoordinatePrecision, py::array::f_style | py::array::forcecast> coordinates, underlying_type<CoefficientPrecision> epsilon, CoordinatePrecision eta) {
         return new DDMSolverBuilder(distributed_operator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections, generator, coordinates.shape()[0], coordinates.data(), epsilon, eta);
-    }));
+    }),py::keep_alive<1, 2>());
 
     ddm_solver_class.def_property_readonly(
         "solver", [](const DDMSolverBuilder &self) { return &self.solver; }, py::return_value_policy::reference_internal);
+    ddm_solver_class.def_property_readonly("local_to_global_numbering",[](const DDMSolverBuilder &self) { return &self.local_to_global_numbering; }, py::return_value_policy::reference_internal);
 }
 #endif
diff --git a/tests/conftest.py b/tests/conftest.py
index db9e62b..1f3f626 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -13,7 +13,7 @@
 from example.define_custom_low_rank_generator import CustomSVD
 
 
-class GeneratorFromMatrix(Htool.ComplexVirtualGeneratorInUserNumbering):
+class GeneratorFromMatrix(Htool.ComplexVirtualGenerator):
     def __init__(self, matrix):
         super().__init__()
         self.matrix = matrix
@@ -27,7 +27,7 @@ def build_submatrix(self, J, K, mat):
                 mat[j, k] = self.get_coef(J[j], K[k])
 
 
-class LocalGeneratorFromMatrix(Htool.ComplexVirtualGeneratorInUserNumbering):
+class LocalGeneratorFromMatrix(Htool.ComplexVirtualGenerator):
     def __init__(
         self,
         permutation,