From 2f528bb11e9206f6050051e41da2950da794a6a0 Mon Sep 17 00:00:00 2001 From: Pierre Marchand Date: Sun, 10 Sep 2023 15:16:18 +0200 Subject: [PATCH] update --- ...block_jacobi_solver_with_default_build.py} | 3 + ...ock_jacobi_solver_with_no_default_build.py | 125 ++++++++++++++++++ lib/htool | 2 +- pyproject.toml | 2 + src/htool/hmatrix/hmatrix.hpp | 6 + src/htool/solver/solver.hpp | 4 +- 6 files changed, 140 insertions(+), 2 deletions(-) rename example/{use_block_jacobi_solver.py => use_block_jacobi_solver_with_default_build.py} (95%) create mode 100644 example/use_block_jacobi_solver_with_no_default_build.py diff --git a/example/use_block_jacobi_solver.py b/example/use_block_jacobi_solver_with_default_build.py similarity index 95% rename from example/use_block_jacobi_solver.py rename to example/use_block_jacobi_solver_with_default_build.py index 1b3a913..f594480 100644 --- a/example/use_block_jacobi_solver.py +++ b/example/use_block_jacobi_solver_with_default_build.py @@ -69,6 +69,8 @@ local_block_hmatrix = default_approximation.block_diagonal_hmatrix print(hmatrix.get_tree_parameters()) print(hmatrix.get_information()) + print(solver.get_information()) + print(solver.get_information("Nb_it")) fig = plt.figure() ax1 = None @@ -89,6 +91,7 @@ ax1.set_title("cluster at depth 1") ax2.set_title("cluster at depth 2") ax4.set_title("Hmatrix on rank 0") + ax4.set_title("Block diagonal Hmatrix on rank 0") Htool.plot(ax1, cluster, points, 1) Htool.plot(ax2, cluster, points, 2) Htool.plot(ax3, hmatrix) diff --git a/example/use_block_jacobi_solver_with_no_default_build.py b/example/use_block_jacobi_solver_with_no_default_build.py new file mode 100644 index 0000000..d1dc9b0 --- /dev/null +++ b/example/use_block_jacobi_solver_with_no_default_build.py @@ -0,0 +1,125 @@ +import matplotlib.pyplot as plt +import mpi4py +import numpy as np +from create_geometry import create_partitionned_geometries +from define_custom_generators import CustomGenerator + +import Htool + +# Random geometry +size = 500 +dimension = 3 +[points, _, partition] = create_partitionned_geometries( + dimension, size, size, mpi4py.MPI.COMM_WORLD.size +) + + +# Htool parameters +eta = 10 +epsilon = 1e-3 +minclustersize = 10 +number_of_children = 2 + +# Build clusters +cluster_builder = Htool.ClusterBuilder() +cluster_builder.set_minclustersize(minclustersize) +cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + points, number_of_children, mpi4py.MPI.COMM_WORLD.size, partition +) + + +# Build generator +generator = CustomGenerator(cluster, points, cluster, points) + +# Build HMatrix +hmatrix_builder = Htool.HMatrixBuilder( + cluster, + cluster, + epsilon, + eta, + "S", + "L", + -1, + mpi4py.MPI.COMM_WORLD.rank, +) + +hmatrix: Htool.HMatrix = hmatrix_builder.build(generator) + + +# Build local operator +local_operator = Htool.LocalHMatrix( + hmatrix, + cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank), + cluster, + "N", + "N", + False, + False, +) + +# Build distributed operator +partition_from_cluster = Htool.PartitionFromCluster(cluster) +distributed_operator = Htool.DistributedOperator( + partition_from_cluster, + partition_from_cluster, + "S", + "L", + mpi4py.MPI.COMM_WORLD, +) + +distributed_operator.add_local_operator(local_operator) + + +# Solver with block Jacobi preconditionner +default_solver_builder = Htool.DefaultSolverBuilder( + distributed_operator, + hmatrix.get_block_diagonal_hmatrix(), +) +solver = default_solver_builder.solver + +# Solver with block Jacobi +x_ref = np.random.random(size) +b = distributed_operator * x_ref +x = np.zeros(size) + +hpddm_args = "-hpddm_compute_residual l2 " +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + hpddm_args += "-hpddm_verbosity 10" +solver.set_hpddm_args(hpddm_args) +solver.facto_one_level() +solver.solve(x, b) + + +# Outputs +if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: + print(np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref)) + print(hmatrix.get_tree_parameters()) + print(hmatrix.get_information()) + print(solver.get_information()) + print(solver.get_information("Nb_it")) + + fig = plt.figure() + ax1 = None + ax2 = None + ax3 = None + ax4 = None + 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) + ax4 = fig.add_subplot(2, 2, 4) + + ax1.set_title("cluster at depth 1") + ax2.set_title("cluster at depth 2") + ax3.set_title("Hmatrix on rank 0") + ax4.set_title("Block diagonal Hmatrix on rank 0") + Htool.plot(ax1, cluster, points, 1) + Htool.plot(ax2, cluster, points, 2) + Htool.plot(ax3, hmatrix) + Htool.plot(ax4, hmatrix.get_block_diagonal_hmatrix()) + plt.show() diff --git a/lib/htool b/lib/htool index f10d16a..44e3fea 160000 --- a/lib/htool +++ b/lib/htool @@ -1 +1 @@ -Subproject commit f10d16a7b84eb99a8fb4028f2a8cb9712c381157 +Subproject commit 44e3feaa62ec80e0d08650988478983b37a3fea2 diff --git a/pyproject.toml b/pyproject.toml index d03184f..b04bc03 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,6 +34,8 @@ line-length = 88 [tool.isort] profile = "black" +known_local_folder = 'Htool,example' +skip = "lib" [build-system] diff --git a/src/htool/hmatrix/hmatrix.hpp b/src/htool/hmatrix/hmatrix.hpp index 16bc6ff..a332082 100644 --- a/src/htool/hmatrix/hmatrix.hpp +++ b/src/htool/hmatrix/hmatrix.hpp @@ -85,6 +85,12 @@ void declare_HMatrix(py::module &m, const std::string &className) { // py_class.def("__str__", [](const Class &self) { // return "HMatrix: (shape: " + htool::NbrToStr(self.nb_cols()) + "x" + htool::NbrToStr(self.nb_rows()) + ", nb_low_rank_blocks: " + htool::NbrToStr(self.get_nlrmat()) + ", nb_dense_blocks: " + htool::NbrToStr(self.get_ndmat()) + ")"; // }); + + py_class.def( + "get_block_diagonal_hmatrix", [](const HMatrix &hmatrix) { + return &*hmatrix.get_diagonal_hmatrix(); + }, + py::return_value_policy::reference_internal); py_class.def("get_tree_parameters", [](const HMatrix &hmatrix) { std::stringstream ss; htool::print_tree_parameters(hmatrix, ss); diff --git a/src/htool/solver/solver.hpp b/src/htool/solver/solver.hpp index 3cbf1b2..2435f40 100644 --- a/src/htool/solver/solver.hpp +++ b/src/htool/solver/solver.hpp @@ -2,6 +2,7 @@ #define HTOOL_DDM_SOLVER_CPP #include "htool/solvers/ddm.hpp" +#include #include namespace py = pybind11; @@ -80,7 +81,8 @@ void declare_DDM(py::module &m, const std::string &className) { opt.parse(hpddm_args); }); // py_class.def("print_infos", &Class::print_infos); - // py_class.def("get_infos", &Class::get_infos); + py_class.def("get_information", py::overload_cast<>(&Class::get_information, py::const_)); + py_class.def("get_information", py::overload_cast(&Class::get_information, py::const_)); } #endif