update solvers

.github/workflows/CI.yml

@@ -133,8 +133,8 @@ jobs:
           fetch-depth: 0 # to get tags
       - name: Checkout hpddm
         run: |
-          git clone https://github.com/PierreMarchand20/hpddm.git hpddm
-          cd hpddm && git checkout 43a77d0afcd169e88a0eb062a517f41ec2b6ede9
+          git clone https://github.com/hpddm/hpddm.git hpddm
+          cd hpddm && git checkout 7e8c98deb030ec078bbb36bee36ea2a1910fc8e3
         # uses: actions/checkout@v3
         # with:
         #   path: "hpddm"

include/htool/distributed_operator/utility.hpp

@@ -7,10 +7,21 @@
 #include "implementations/partition_from_cluster.hpp"
 namespace htool {
 
+template <typename CoefficientPrecision, typename CoordinatePrecision = underlying_type<CoefficientPrecision>>
+class CustomApproximationBuilder {
+    const PartitionFromCluster<CoefficientPrecision, CoordinatePrecision> target_partition, source_partition;
+
+  public:
+    DistributedOperator<CoefficientPrecision> distributed_operator;
+
+    explicit CustomApproximationBuilder(const Cluster<CoordinatePrecision> &target_cluster, const Cluster<CoordinatePrecision> &source_cluster, char symmetry, char UPLO, MPI_Comm communicator, const VirtualLocalOperator<CoefficientPrecision> &local_operator) : target_partition(target_cluster), source_partition(source_cluster), distributed_operator(target_partition, source_partition, symmetry, UPLO, communicator) {
+        distributed_operator.add_local_operator(&local_operator);
+    }
+};
+
 template <typename CoefficientPrecision, typename CoordinatePrecision = underlying_type<CoefficientPrecision>>
 class DistributedOperatorFromHMatrix {
   private:
-    const PartitionFromCluster<CoefficientPrecision, CoordinatePrecision> target_partition, source_partition;
     std::function<int(MPI_Comm)> get_rankWorld = [](MPI_Comm comm) {
     int rankWorld;
     MPI_Comm_rank(comm, &rankWorld); 
@@ -21,13 +32,13 @@ class DistributedOperatorFromHMatrix {
 
   private:
     const LocalHMatrix<CoefficientPrecision, CoordinatePrecision> local_hmatrix;
+    CustomApproximationBuilder<CoefficientPrecision> distributed_operator_holder;
 
   public:
-    DistributedOperator<CoefficientPrecision> distributed_operator;
+    DistributedOperator<CoefficientPrecision> &distributed_operator;
     const HMatrix<CoefficientPrecision, CoordinatePrecision> *block_diagonal_hmatrix{nullptr};
 
-    DistributedOperatorFromHMatrix(const VirtualGenerator<CoefficientPrecision> &generator, const Cluster<CoordinatePrecision> &target_cluster, const Cluster<CoordinatePrecision> &source_cluster, const HMatrixTreeBuilder<CoefficientPrecision, CoordinatePrecision> &hmatrix_builder, MPI_Comm communicator) : target_partition(target_cluster), source_partition(source_cluster), hmatrix(hmatrix_builder.build(generator)), local_hmatrix(hmatrix, hmatrix_builder.get_target_cluster().get_cluster_on_partition(get_rankWorld(communicator)), hmatrix_builder.get_source_cluster(), hmatrix_builder.get_symmetry(), hmatrix_builder.get_UPLO(), false, false), distributed_operator(target_partition, source_partition, hmatrix_builder.get_symmetry(), hmatrix_builder.get_UPLO(), communicator) {
-        distributed_operator.add_local_operator(&local_hmatrix);
+    DistributedOperatorFromHMatrix(const VirtualGenerator<CoefficientPrecision> &generator, const Cluster<CoordinatePrecision> &target_cluster, const Cluster<CoordinatePrecision> &source_cluster, const HMatrixTreeBuilder<CoefficientPrecision, CoordinatePrecision> &hmatrix_builder, MPI_Comm communicator) : hmatrix(hmatrix_builder.build(generator)), local_hmatrix(hmatrix, hmatrix_builder.get_target_cluster().get_cluster_on_partition(get_rankWorld(communicator)), hmatrix_builder.get_source_cluster(), hmatrix_builder.get_symmetry(), hmatrix_builder.get_UPLO(), false, false), distributed_operator_holder(target_cluster, source_cluster, hmatrix_builder.get_symmetry(), hmatrix_builder.get_UPLO(), communicator, local_hmatrix), distributed_operator(distributed_operator_holder.distributed_operator) {
         block_diagonal_hmatrix = hmatrix.get_sub_hmatrix(hmatrix_builder.get_target_cluster().get_cluster_on_partition(get_rankWorld(communicator)), hmatrix_builder.get_source_cluster().get_cluster_on_partition(get_rankWorld(communicator)));
     }
 

include/htool/solvers/ddm.hpp

@@ -463,7 +463,7 @@ DDM<CoefficientPrecision, HPDDMCustomLocalSolver> make_DDM_solver_w_custom_local
     hpddm_op->initialize(n, sym, nullptr, neighbors, intersections); // we should not give a local dense matrix
 
     auto local_hmatrix_solver = std::make_unique<LocalHMatrixSolver<CoefficientPrecision, CoordinatePrecision>>(local_hmatrix, use_permutation);
-    hpddm_op->get_local_solver().set_local_solver(std::move(local_hmatrix_solver));
+    hpddm_op->getSolver().set_local_solver(std::move(local_hmatrix_solver));
 
     mytime = MPI_Wtime() - time;
 
@@ -494,7 +494,7 @@ DDM<CoefficientPrecision, HPDDMCustomLocalSolver> make_DDM_solver_w_custom_local
     hpddm_op->initialize(n, sym, nullptr, neighbors, intersections); // we should not give a local dense matrix
 
     auto local_hmatrix_solver = std::make_unique<LocalHMatrixPlusOverlapSolver<CoefficientPrecision, CoordinatePrecision>>(local_hmatrix, B, C, D);
-    hpddm_op->get_local_solver().set_local_solver(std::move(local_hmatrix_solver));
+    hpddm_op->getSolver().set_local_solver(std::move(local_hmatrix_solver));
 
     mytime = MPI_Wtime() - time;
 

include/htool/solvers/utility.hpp

@@ -184,7 +184,7 @@ class DDMSolverWithDenseLocalSolver {
 };
 
 template <typename CoefficientPrecision, typename CoordinatePrecision = underlying_type<CoefficientPrecision>>
-class DDMSolver {
+class DDMSolverBuilder {
   private:
     std::vector<int> m_neighbors;
     std::vector<std::vector<int>> m_intersections;
@@ -273,13 +273,13 @@ class DDMSolver {
     DDM<CoefficientPrecision, HPDDMCustomLocalSolver> solver;
 
     // Block Jacobi
-    DDMSolver(DistributedOperator<CoefficientPrecision> &distributed_operator, HMatrix<CoefficientPrecision, CoordinatePrecision> &block_diagonal_hmatrix) : local_to_global_numbering(block_diagonal_hmatrix.get_target_cluster().get_permutation()), solver(make_DDM_solver_w_custom_local_solver(distributed_operator, block_diagonal_hmatrix, m_neighbors, m_intersections, false)) {}
+    DDMSolverBuilder(DistributedOperator<CoefficientPrecision> &distributed_operator, HMatrix<CoefficientPrecision, CoordinatePrecision> &block_diagonal_hmatrix) : local_to_global_numbering(block_diagonal_hmatrix.get_target_cluster().get_permutation()), solver(make_DDM_solver_w_custom_local_solver(distributed_operator, block_diagonal_hmatrix, m_neighbors, m_intersections, false)) {}
 
     // DDM building local hmatrix adding overlap
-    DDMSolver(DistributedOperator<CoefficientPrecision> &distributed_operator, HMatrix<CoefficientPrecision, CoordinatePrecision> &block_diagonal_hmatrix, const VirtualGeneratorInUserNumbering<CoefficientPrecision> &generator, 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) : m_local_numbering(ovr_subdomain_to_global, cluster_to_ovr_subdomain, intersections), local_to_global_numbering(m_local_numbering.local_to_global_numbering), blocks_in_overlap(initialize_blocks_in_overlap(distributed_operator, block_diagonal_hmatrix, generator)), solver(make_DDM_solver_w_custom_local_solver(distributed_operator, block_diagonal_hmatrix, blocks_in_overlap[0], blocks_in_overlap[1], blocks_in_overlap[2], neighbors, m_local_numbering.intersections)) {}
+    DDMSolverBuilder(DistributedOperator<CoefficientPrecision> &distributed_operator, HMatrix<CoefficientPrecision, CoordinatePrecision> &block_diagonal_hmatrix, const VirtualGeneratorInUserNumbering<CoefficientPrecision> &generator, 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) : m_local_numbering(ovr_subdomain_to_global, cluster_to_ovr_subdomain, intersections), local_to_global_numbering(m_local_numbering.local_to_global_numbering), blocks_in_overlap(initialize_blocks_in_overlap(distributed_operator, block_diagonal_hmatrix, generator)), solver(make_DDM_solver_w_custom_local_solver(distributed_operator, block_diagonal_hmatrix, blocks_in_overlap[0], blocks_in_overlap[1], blocks_in_overlap[2], neighbors, m_local_numbering.intersections)) {}
 
     // DDM building local hmatrix with overlap
-    DDMSolver(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, int spatial_dimension, const CoordinatePrecision *global_geometry, underlying_type<CoefficientPrecision> epsilon, CoordinatePrecision eta) : m_local_numbering(ovr_subdomain_to_global, cluster_to_ovr_subdomain, intersections), local_to_global_numbering(m_local_numbering.local_to_global_numbering), local_cluster(std::make_unique<Cluster<CoordinatePrecision>>(initialize_local_cluster(spatial_dimension, global_geometry))), local_hmatrix(std::make_unique<HMatrix<CoefficientPrecision, CoordinatePrecision>>(initialize_local_hmatrix(generator, epsilon, eta, distributed_operator.get_symmetry_type()))), solver(make_DDM_solver_w_custom_local_solver(distributed_operator, *local_hmatrix, neighbors, m_local_numbering.intersections, true)) {}
+    DDMSolverBuilder(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, int spatial_dimension, const CoordinatePrecision *global_geometry, underlying_type<CoefficientPrecision> epsilon, CoordinatePrecision eta) : m_local_numbering(ovr_subdomain_to_global, cluster_to_ovr_subdomain, intersections), local_to_global_numbering(m_local_numbering.local_to_global_numbering), local_cluster(std::make_unique<Cluster<CoordinatePrecision>>(initialize_local_cluster(spatial_dimension, global_geometry))), local_hmatrix(std::make_unique<HMatrix<CoefficientPrecision, CoordinatePrecision>>(initialize_local_hmatrix(generator, epsilon, eta, distributed_operator.get_symmetry_type()))), solver(make_DDM_solver_w_custom_local_solver(distributed_operator, *local_hmatrix, neighbors, m_local_numbering.intersections, true)) {}
 };
 
 } // namespace htool

tests/functional_tests/solvers/test_solver.cpp

@@ -30,9 +30,9 @@ int main(int argc, char *argv[]) {
                 is_error = is_error || test_solver_wo_overlap<DDMSolverWithDenseLocalSolver<complex<double>, double>>(argc, argv, nb_rhs, symmetry, datapath_final);
                 is_error = is_error || test_solver_ddm_adding_overlap<DDMSolverWithDenseLocalSolver<std::complex<double>>>(argc, argv, nb_rhs, data_symmetry, symmetry, datapath_final);
                 is_error = is_error || test_solver_ddm<DDMSolverWithDenseLocalSolver<complex<double>>>(argc, argv, nb_rhs, data_symmetry, symmetry, datapath_final);
-                is_error = is_error || test_solver_wo_overlap<DDMSolver<complex<double>, double>>(argc, argv, nb_rhs, symmetry, datapath_final);
-                test_solver_ddm_adding_overlap<DDMSolver<std::complex<double>>>(argc, argv, nb_rhs, data_symmetry, symmetry, datapath_final);
-                is_error = is_error || test_solver_ddm<DDMSolver<complex<double>>>(argc, argv, nb_rhs, data_symmetry, symmetry, datapath_final);
+                is_error = is_error || test_solver_wo_overlap<DDMSolverBuilder<complex<double>, double>>(argc, argv, nb_rhs, symmetry, datapath_final);
+                test_solver_ddm_adding_overlap<DDMSolverBuilder<std::complex<double>>>(argc, argv, nb_rhs, data_symmetry, symmetry, datapath_final);
+                is_error = is_error || test_solver_ddm<DDMSolverBuilder<complex<double>>>(argc, argv, nb_rhs, data_symmetry, symmetry, datapath_final);
             }
         }
     }

tests/functional_tests/solvers/test_solver_ddm.hpp

@@ -232,7 +232,7 @@ int test_solver_ddm(int argc, char *argv[], int mu, char data_symmetry, char sym
 }
 
 template <>
-int test_solver_ddm<DDMSolver<complex<double>>>(int argc, char *argv[], int mu, char, char symmetric, std::string datapath) {
+int test_solver_ddm<DDMSolverBuilder<complex<double>>>(int argc, char *argv[], int mu, char, char symmetric, std::string datapath) {
 
     // Get the number of processes
     int size;
@@ -324,7 +324,7 @@ int test_solver_ddm<DDMSolver<complex<double>>>(int argc, char *argv[], int mu,
         std::cout << "Creating HMatrix" << std::endl;
     std::vector<double> geometry(n);
     bytes_to_vector(geometry, datapath + "/geometry.bin");
-    DDMSolver<complex<double>> default_ddm_solver(Operator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections, Generator, 3, geometry.data(), epsilon, eta);
+    DDMSolverBuilder<complex<double>> default_ddm_solver(Operator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections, Generator, 3, geometry.data(), epsilon, eta);
     auto &ddm_with_overlap = default_ddm_solver.solver;
 
     // No precond with overlap

tests/functional_tests/solvers/test_solver_ddm_adding_overlap.hpp

@@ -230,7 +230,7 @@ int test_solver_ddm_adding_overlap(int argc, char *argv[], int mu, char data_sym
 }
 
 template <>
-int test_solver_ddm_adding_overlap<DDMSolver<std::complex<double>>>(int argc, char *argv[], int mu, char, char symmetric, std::string datapath) {
+int test_solver_ddm_adding_overlap<DDMSolverBuilder<std::complex<double>>>(int argc, char *argv[], int mu, char, char symmetric, std::string datapath) {
 
     // Get the number of processes
     int size;
@@ -320,7 +320,7 @@ int test_solver_ddm_adding_overlap<DDMSolver<std::complex<double>>>(int argc, ch
     // Solve
     if (rank == 0)
         std::cout << "Creating HMatrix" << std::endl;
-    DDMSolver<std::complex<double>> default_ddm_solver(Operator, local_block_diagonal_hmatrix, Generator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections);
+    DDMSolverBuilder<std::complex<double>> default_ddm_solver(Operator, local_block_diagonal_hmatrix, Generator, ovr_subdomain_to_global, cluster_to_ovr_subdomain, neighbors, intersections);
     auto &ddm_with_overlap = default_ddm_solver.solver;
 
     // No precond with overlap