From 2207e6c4b89a1ff68b7e24f0ed7b297680c042ca Mon Sep 17 00:00:00 2001 From: Pierre Marchand Date: Wed, 7 Aug 2024 19:52:36 +0200 Subject: [PATCH] update htool 0.9.0 --- CMakeLists.txt | 2 +- example/ddm_solver.py | 97 -------- example/define_custom_generators.py | 61 +---- example/define_custom_local_operator.py | 12 +- example/define_custom_low_rank_generator.py | 13 +- ...lt_build.py => use_block_jacobi_solver.py} | 15 +- ...ock_jacobi_solver_with_no_default_build.py | 134 ---------- example/use_cluster.py | 2 - example/use_cluster_with_given_partition.py | 2 - example/use_custom_dense_block_generator.py | 4 +- example/use_custom_local_operator.py | 20 +- example/use_custom_low_rank_approximation.py | 35 +-- ...lt_build.py => use_hmatrix_compression.py} | 9 +- ...ld.py => use_local_hmatrix_compression.py} | 28 +-- example/use_no_default_build.py | 128 ---------- lib/hpddm | 2 +- lib/htool | 2 +- pyproject.toml | 2 +- src/htool/clustering/cluster_node.hpp | 15 -- .../distributed_operator.hpp | 39 ++- .../implementation/partition_from_cluster.hpp | 16 -- .../interface/partition.hpp | 32 --- src/htool/distributed_operator/utility.hpp | 32 ++- src/htool/hmatrix/hmatrix.hpp | 122 ++------- src/htool/hmatrix/hmatrix_builder.hpp | 6 +- .../hmatrix/interfaces/virtual_generator.hpp | 61 +---- .../interfaces/virtual_low_rank_generator.hpp | 45 ++-- src/htool/local_operator/local_operator.hpp | 8 +- src/htool/main.cpp | 21 +- src/htool/misc/utility.hpp | 2 +- .../solver/geneo/coarse_space_builder.hpp | 35 ++- .../geneo/coarse_space_dense_builder.hpp | 133 ++++++++++ src/htool/solver/solver.hpp | 19 +- src/htool/solver/utility.hpp | 51 +--- tests/conftest.py | 59 ++--- tests/test_cluster.py | 9 +- tests/test_ddm_solver.py | 232 ++++++++---------- 37 files changed, 455 insertions(+), 1050 deletions(-) delete mode 100644 example/ddm_solver.py rename example/{use_block_jacobi_solver_with_default_build.py => use_block_jacobi_solver.py} (89%) delete mode 100644 example/use_block_jacobi_solver_with_no_default_build.py rename example/{use_default_build.py => use_hmatrix_compression.py} (94%) rename example/{use_local_build.py => use_local_hmatrix_compression.py} (91%) delete mode 100644 example/use_no_default_build.py delete mode 100644 src/htool/distributed_operator/implementation/partition_from_cluster.hpp delete mode 100644 src/htool/distributed_operator/interface/partition.hpp create mode 100644 src/htool/solver/geneo/coarse_space_dense_builder.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index e41773d..e8419d0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -87,7 +87,7 @@ if("${BLA_VENDOR}" STREQUAL "Intel10_32" target_compile_definitions(htool PRIVATE "-DHPDDM_MKL -DHTOOL_MKL") endif() -target_compile_definitions(Htool PRIVATE "-DPYTHON_INTERFACE" "-DWITH_HPDDM") +target_compile_definitions(Htool PRIVATE "-DHTOOL_WITH_PYTHON_INTERFACE" "-DHTOOL_WITH_HPDDM") if(CODE_COVERAGE AND (CMAKE_C_COMPILER_ID MATCHES "GNU" OR CMAKE_CXX_COMPILER_ID MATCHES "GNU")) target_compile_options(Htool PRIVATE -fprofile-arcs -ftest-coverage) diff --git a/example/ddm_solver.py b/example/ddm_solver.py deleted file mode 100644 index 8dd208d..0000000 --- a/example/ddm_solver.py +++ /dev/null @@ -1,97 +0,0 @@ -import Htool -import matplotlib.pyplot as plt -import mpi4py -import numpy as np -from create_geometry import create_random_geometries -from define_custom_generators import CustomGeneratorWithPermutation - -# Random geometry -size = 500 -dimension = 3 - -[points, _] = create_random_geometries(dimension, size, 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 -) - -# Build generator -generator = CustomGeneratorWithPermutation( - cluster.get_permutation(), points, cluster.get_permutation(), points -) - -# Build distributed operator -default_approximation = Htool.DefaultApproximationBuilder( - generator, - cluster, - cluster, - epsilon, - eta, - "S", - "L", - mpi4py.MPI.COMM_WORLD, -) - -# Solver with block Jacobi preconditionner -default_solver_builder = Htool.DefaultSolverBuilder( - default_approximation.distributed_operator, - default_approximation.block_diagonal_hmatrix, -) -solver = default_solver_builder.solver - - -# Solver with block Jacobi -x_ref = np.random.random(size) -b = default_approximation.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) - - -# Several ways to display information -if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: - print(np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref)) - hmatrix = default_approximation.hmatrix - local_block_hmatrix = default_approximation.block_diagonal_hmatrix - print(hmatrix.get_tree_parameters()) - print(hmatrix.get_information()) - - 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") - ax4.set_title("Hmatrix on rank 0") - Htool.plot(ax1, cluster, points, 1) - Htool.plot(ax2, cluster, points, 2) - Htool.plot(ax3, hmatrix) - Htool.plot(ax4, local_block_hmatrix) - plt.show() diff --git a/example/define_custom_generators.py b/example/define_custom_generators.py index 334e7a6..2519eab 100644 --- a/example/define_custom_generators.py +++ b/example/define_custom_generators.py @@ -2,66 +2,13 @@ import numpy as np -class CustomGenerator(Htool.VirtualGenerator): - def __init__(self, target_cluster, target_points, source_cluster, source_points): +class CustomGenerator(Htool.VirtualGeneratorInUserNumbering): + def __init__(self, target_points, source_points): super().__init__() self.target_points = target_points - self.target_permutation = target_cluster.get_permutation() self.source_points = source_points - self.source_permutation = source_cluster.get_permutation() - self.nb_rows = len(self.target_permutation) - self.nb_cols = len(self.source_permutation) - - def get_coef(self, i, j): - return 1.0 / ( - 1e-1 - + np.linalg.norm( - self.target_points[:, self.target_permutation[i]] - - self.source_points[:, self.source_permutation[j]] - ) - ) - - def build_submatrix(self, row_offset, col_offset, mat): - for j in range(0, mat.shape[0]): - for k in range(0, mat.shape[1]): - mat[j, k] = 1.0 / ( - 1.0e-1 - + np.linalg.norm( - self.target_points[:, self.target_permutation[j + row_offset]] - - self.source_points[:, self.source_permutation[k + col_offset]] - ) - ) - - def mat_vec(self, x): - y = np.zeros(self.nb_rows) - for i in range(0, self.nb_rows): - for j in range(0, self.nb_cols): - y[self.target_permutation[i]] += ( - self.get_coef(i, j) * x[self.source_permutation[j]] - ) - return y - - def mat_mat(self, X): - Y = np.zeros((self.nb_rows, X.shape[1])) - - for i in range(0, self.nb_rows): - for j in range(0, X.shape[1]): - for k in range(0, self.nb_cols): - Y[self.target_permutation[i], j] += ( - self.get_coef(i, k) * X[self.source_permutation[k], j] - ) - return Y - - -class CustomGeneratorWithPermutation(Htool.VirtualGeneratorWithPermutation): - def __init__( - self, target_permutation, target_points, source_permutation, source_points - ): - super().__init__(target_permutation, source_permutation) - self.target_points = target_points - self.source_points = source_points - self.nb_rows = len(target_permutation) - self.nb_cols = len(source_permutation) + self.nb_rows = target_points.shape[1] + self.nb_cols = source_points.shape[1] def get_coef(self, i, j): return 1.0 / ( diff --git a/example/define_custom_local_operator.py b/example/define_custom_local_operator.py index c28c4ce..526a38b 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.VirtualGenerator, + generator: Htool.VirtualGeneratorInUserNumbering, target_cluster: Htool.Cluster, source_cluster: Htool.Cluster, symmetry: str = "N", @@ -23,7 +23,15 @@ def __init__( ) self.data = np.zeros((target_cluster.get_size(), source_cluster.get_size())) generator.build_submatrix( - target_cluster.get_offset(), source_cluster.get_offset(), self.data + target_cluster.get_permutation()[ + target_cluster.get_offset() : target_cluster.get_offset() + + target_cluster.get_size() + ], + source_cluster.get_permutation()[ + source_cluster.get_offset() : source_cluster.get_offset() + + source_cluster.get_size() + ], + self.data, ) def add_vector_product( diff --git a/example/define_custom_low_rank_generator.py b/example/define_custom_low_rank_generator.py index 6b02f5b..dab4378 100644 --- a/example/define_custom_low_rank_generator.py +++ b/example/define_custom_low_rank_generator.py @@ -5,11 +5,12 @@ class CustomSVD(Htool.VirtualLowRankGenerator): - def build_low_rank_approximation( - self, generator, target_size, source_size, target_offset, source_offset, epsilon - ): - submat = np.zeros((target_size, source_size), order="F") - generator.build_submatrix(target_offset, source_offset, submat) + def __init__(self, generator: Htool.VirtualGeneratorInUserNumbering): + super().__init__(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) u, s, vh = np.linalg.svd(submat, full_matrices=False) norm = np.linalg.norm(submat) @@ -18,7 +19,7 @@ def build_low_rank_approximation( while count > 0 and math.sqrt(svd_norm) / norm < epsilon: svd_norm += s[count] ** 2 count -= 1 - count = min(count + 1, min(target_size, source_size)) + 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_with_default_build.py b/example/use_block_jacobi_solver.py similarity index 89% rename from example/use_block_jacobi_solver_with_default_build.py rename to example/use_block_jacobi_solver.py index d7352ed..d3301a9 100644 --- a/example/use_block_jacobi_solver_with_default_build.py +++ b/example/use_block_jacobi_solver.py @@ -1,3 +1,5 @@ +import copy + import Htool import matplotlib.pyplot as plt import mpi4py @@ -26,7 +28,7 @@ ) # Build generator -generator = CustomGenerator(cluster, points, cluster, points) +generator = CustomGenerator(points, points) # Build distributed operator default_approximation = Htool.DefaultApproximationBuilder( @@ -41,9 +43,10 @@ ) # Solver with block Jacobi preconditionner -default_solver_builder = Htool.DefaultSolverBuilder( - default_approximation.distributed_operator, - default_approximation.block_diagonal_hmatrix, +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 ) solver = default_solver_builder.solver @@ -79,10 +82,6 @@ print(solver_information) 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) diff --git a/example/use_block_jacobi_solver_with_no_default_build.py b/example/use_block_jacobi_solver_with_no_default_build.py deleted file mode 100644 index e5ccc11..0000000 --- a/example/use_block_jacobi_solver_with_no_default_build.py +++ /dev/null @@ -1,134 +0,0 @@ -import Htool -import matplotlib.pyplot as plt -import mpi4py -import numpy as np -from create_geometry import create_partitionned_geometries -from define_custom_generators import CustomGenerator - -# 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 -block_diagonal_hmatrix = hmatrix.get_sub_hmatrix( - cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank), - cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank), -) -default_solver_builder = Htool.DefaultSolverBuilder( - distributed_operator, - 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 -hmatrix_distributed_information = hmatrix.get_distributed_information( - mpi4py.MPI.COMM_WORLD -) -hmatrix_tree_parameter = hmatrix.get_tree_parameters() -hmatrix_local_information = hmatrix.get_local_information() -solver_information = solver.get_information() -if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: - print(np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref)) - print(hmatrix_distributed_information) - print(hmatrix_local_information) - print(hmatrix_tree_parameter) - print(solver_information) - - 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, block_diagonal_hmatrix) - plt.show() diff --git a/example/use_cluster.py b/example/use_cluster.py index e975e89..1f3589a 100644 --- a/example/use_cluster.py +++ b/example/use_cluster.py @@ -41,8 +41,6 @@ if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: fig = plt.figure() - ax1 = None - ax2 = None if dimension == 2: ax1 = fig.add_subplot(1, 2, 1) diff --git a/example/use_cluster_with_given_partition.py b/example/use_cluster_with_given_partition.py index 2306eae..8499605 100644 --- a/example/use_cluster_with_given_partition.py +++ b/example/use_cluster_with_given_partition.py @@ -31,8 +31,6 @@ if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: fig = plt.figure() - ax1 = None - ax2 = None if dimension == 2: ax1 = fig.add_subplot(1, 2, 1) diff --git a/example/use_custom_dense_block_generator.py b/example/use_custom_dense_block_generator.py index 5514172..05639ba 100644 --- a/example/use_custom_dense_block_generator.py +++ b/example/use_custom_dense_block_generator.py @@ -32,9 +32,7 @@ # Build generator -generator = CustomGenerator( - target_cluster, target_points, source_cluster, source_points -) +generator = CustomGenerator(target_points, source_points) # custom_dense_blocks_generator = CustomDenseBlocksGenerator(generator) diff --git a/example/use_custom_local_operator.py b/example/use_custom_local_operator.py index 9d952d3..d9814fe 100644 --- a/example/use_custom_local_operator.py +++ b/example/use_custom_local_operator.py @@ -32,9 +32,7 @@ # Build generator -generator = CustomGenerator( - target_cluster, target_points, source_cluster, source_points -) +generator = CustomGenerator(target_points, source_points) # Build local operator local_operator = CustomLocalOperator( @@ -44,22 +42,14 @@ "N", "N", False, - False, + True, ) # Build distributed operator -target_partition_from_cluster = Htool.PartitionFromCluster(target_cluster) -source_partition_from_cluster = Htool.PartitionFromCluster(source_cluster) -distributed_operator = Htool.DistributedOperator( - target_partition_from_cluster, - source_partition_from_cluster, - "N", - "N", - mpi4py.MPI.COMM_WORLD, +custom_local_approximation = Htool.CustomApproximationBuilder( + target_cluster, source_cluster, "N", "N", mpi4py.MPI.COMM_WORLD, local_operator ) - -distributed_operator.add_local_operator(local_operator) - +distributed_operator = custom_local_approximation.distributed_operator # Test matrix vector product np.random.seed(0) diff --git a/example/use_custom_low_rank_approximation.py b/example/use_custom_low_rank_approximation.py index 1ed76b5..ac05a9c 100644 --- a/example/use_custom_low_rank_approximation.py +++ b/example/use_custom_low_rank_approximation.py @@ -32,13 +32,11 @@ # Build generator -generator = CustomGenerator( - target_cluster, target_points, source_cluster, source_points -) +generator = CustomGenerator(target_points, source_points) # Low rank generator -low_rank_generator = CustomSVD() +low_rank_generator = CustomSVD(generator) # Build HMatrix hmatrix_builder = Htool.HMatrixBuilder( @@ -50,37 +48,16 @@ "N", -1, mpi4py.MPI.COMM_WORLD.rank, + mpi4py.MPI.COMM_WORLD.rank, ) hmatrix_builder.set_low_rank_generator(low_rank_generator) -hmatrix: Htool.HMatrix = hmatrix_builder.build(generator) -del hmatrix_builder -low_rank_generator.test() -del low_rank_generator - - -# Build local operator -local_operator = Htool.LocalHMatrix( - hmatrix, - target_cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank), - source_cluster, - "N", - "N", - False, - False, -) # Build distributed operator -target_partition_from_cluster = Htool.PartitionFromCluster(target_cluster) -source_partition_from_cluster = Htool.PartitionFromCluster(source_cluster) -distributed_operator = Htool.DistributedOperator( - target_partition_from_cluster, - source_partition_from_cluster, - "N", - "N", - mpi4py.MPI.COMM_WORLD, +distributed_operator_from_hmatrix = Htool.DistributedOperatorFromHMatrix( + generator, target_cluster, source_cluster, hmatrix_builder, mpi4py.MPI.COMM_WORLD ) -distributed_operator.add_local_operator(local_operator) +distributed_operator = distributed_operator_from_hmatrix.distributed_operator # Test matrix vector product diff --git a/example/use_default_build.py b/example/use_hmatrix_compression.py similarity index 94% rename from example/use_default_build.py rename to example/use_hmatrix_compression.py index e8278af..932a118 100644 --- a/example/use_default_build.py +++ b/example/use_hmatrix_compression.py @@ -32,9 +32,7 @@ # Build generator -generator = CustomGenerator( - target_cluster, target_points, source_cluster, source_points -) +generator = CustomGenerator(target_points, source_points) # Build distributed operator default_approximation = Htool.DefaultApproximationBuilder( @@ -77,10 +75,7 @@ print(hmatrix_tree_parameter) 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) diff --git a/example/use_local_build.py b/example/use_local_hmatrix_compression.py similarity index 91% rename from example/use_local_build.py rename to example/use_local_hmatrix_compression.py index 2607705..60d9e23 100644 --- a/example/use_local_build.py +++ b/example/use_local_hmatrix_compression.py @@ -53,9 +53,7 @@ permuted_source_points[:, i] = source_points[:, source_permutation[i]] # Build generator -generator = CustomGenerator( - target_cluster, target_points, source_cluster, source_points -) +generator = CustomGenerator(target_points, source_points) # Build distributed operator @@ -90,12 +88,7 @@ permuted_source_points, number_of_children, 2, off_diagonal_partition ) -off_diagonal_generator = CustomGenerator( - target_cluster, - target_points, - off_diagonal_cluster, - permuted_source_points, -) +off_diagonal_generator = CustomGenerator(target_points, permuted_source_points) local_operator_1 = None if off_diagonal_nc_1 > 0: @@ -154,10 +147,6 @@ print(hmatrix_tree_parameter) 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) @@ -175,11 +164,12 @@ ax4.set_title("Hmatrix on rank 0") Htool.plot(ax1, source_cluster, source_points, 1) Htool.plot(ax2, source_cluster, source_points, 2) - Htool.plot( - ax3, - off_diagonal_cluster.get_cluster_on_partition(1), - permuted_source_points, - 2, - ) + if mpi4py.MPI.COMM_WORLD.Get_size() > 1: + Htool.plot( + ax3, + off_diagonal_cluster.get_cluster_on_partition(1), + permuted_source_points, + 2, + ) Htool.plot(ax4, hmatrix) plt.show() diff --git a/example/use_no_default_build.py b/example/use_no_default_build.py deleted file mode 100644 index 6d1b1c8..0000000 --- a/example/use_no_default_build.py +++ /dev/null @@ -1,128 +0,0 @@ -import Htool -import matplotlib.pyplot as plt -import mpi4py -import numpy as np -from create_geometry import create_partitionned_geometries -from define_custom_generators import CustomGenerator - -# Random geometry -nb_rows = 500 -nb_cols = 500 -dimension = 3 -[target_points, source_points, target_partition] = create_partitionned_geometries( - dimension, nb_rows, nb_cols, 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) -target_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( - target_points, number_of_children, mpi4py.MPI.COMM_WORLD.size, target_partition -) -source_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( - source_points, number_of_children, mpi4py.MPI.COMM_WORLD.size -) - - -# Build generator -generator = CustomGenerator( - target_cluster, target_points, source_cluster, source_points -) - -# Build HMatrix -hmatrix_builder = Htool.HMatrixBuilder( - target_cluster, - source_cluster, - epsilon, - eta, - "N", - "N", - -1, - mpi4py.MPI.COMM_WORLD.rank, -) - -hmatrix: Htool.HMatrix = hmatrix_builder.build(generator) - - -# Build local operator -local_operator = Htool.LocalHMatrix( - hmatrix, - target_cluster.get_cluster_on_partition(mpi4py.MPI.COMM_WORLD.rank), - source_cluster, - "N", - "N", - False, - False, -) - -# Build distributed operator -target_partition_from_cluster = Htool.PartitionFromCluster(target_cluster) -source_partition_from_cluster = Htool.PartitionFromCluster(source_cluster) -distributed_operator = Htool.DistributedOperator( - target_partition_from_cluster, - source_partition_from_cluster, - "N", - "N", - mpi4py.MPI.COMM_WORLD, -) - -distributed_operator.add_local_operator(local_operator) - - -# Test matrix vector product -np.random.seed(0) -x = np.random.rand(nb_cols) -y_1 = distributed_operator * x -y_2 = generator.mat_vec(x) -print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(y_1 - y_2) / np.linalg.norm(y_2)) - - -# Test matrix matrix product -X = np.asfortranarray(np.random.rand(nb_cols, 2)) -Y_1 = distributed_operator @ X -Y_2 = generator.mat_mat(X) -print(mpi4py.MPI.COMM_WORLD.rank, np.linalg.norm(Y_1 - Y_2) / np.linalg.norm(Y_2)) - - -# Outputs -hmatrix_distributed_information = hmatrix.get_distributed_information( - mpi4py.MPI.COMM_WORLD -) -hmatrix_tree_parameter = hmatrix.get_tree_parameters() -hmatrix_local_information = hmatrix.get_local_information() -if mpi4py.MPI.COMM_WORLD.Get_rank() == 0: - print(hmatrix_distributed_information) - print(hmatrix_local_information) - print(hmatrix_tree_parameter) - 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, projection="3d") - ax4 = fig.add_subplot(2, 2, 4) - - ax1.set_title("target cluster at depth 1") - ax2.set_title("target cluster at depth 2") - ax3.set_title("source cluster at depth 1") - ax4.set_title("Hmatrix on rank 0") - Htool.plot(ax1, target_cluster, target_points, 1) - Htool.plot(ax2, target_cluster, target_points, 2) - Htool.plot(ax3, source_cluster, source_points, 1) - Htool.plot(ax4, hmatrix) - plt.show() diff --git a/lib/hpddm b/lib/hpddm index 55c0af0..5890d5a 160000 --- a/lib/hpddm +++ b/lib/hpddm @@ -1 +1 @@ -Subproject commit 55c0af02e5b78d2cb1e466b684e06ce2f111d2ac +Subproject commit 5890d5addf3962d539dc25c441ec3ff4af93b3ab diff --git a/lib/htool b/lib/htool index dd96716..29ce0df 160000 --- a/lib/htool +++ b/lib/htool @@ -1 +1 @@ -Subproject commit dd96716e2d522b2ca50949ddc1c15581c0d0a8f8 +Subproject commit 29ce0dfe8125b6ebbd6efc1dfff6229dbaf92f38 diff --git a/pyproject.toml b/pyproject.toml index 3ae84e8..5537558 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -21,7 +21,7 @@ dependencies = ["mpi4py", "numpy"] [project.optional-dependencies] -dev = ["ruff", "matplotlib>=3.0.0", "pytest"] +dev = ["ruff", "matplotlib>=3.0.0", "pytest", "scipy"] [project.urls] diff --git a/src/htool/clustering/cluster_node.hpp b/src/htool/clustering/cluster_node.hpp index abfdf92..fdd5861 100644 --- a/src/htool/clustering/cluster_node.hpp +++ b/src/htool/clustering/cluster_node.hpp @@ -20,21 +20,6 @@ void declare_cluster_node(py::module &m, const std::string &className) { py_class.def("get_offset", &Class::get_offset); py_class.def("get_permutation", [](const Class &self) { auto &permutation = self.get_permutation(); - // int rankWorld; - // MPI_Comm_rank(MPI_COMM_WORLD, &rankWorld); - // if (rankWorld == 0) { - // for (auto &elt : permutation) { - // std::cout << elt << " "; - // } - // std::cout << "\n"; - // } - // MPI_Barrier(MPI_COMM_WORLD); - // if (rankWorld == 1) { - // for (auto &elt : permutation) { - // std::cout << elt << " "; - // } - // std::cout << "\n"; - // } return py::array_t(std::array{permutation.size()}, permutation.data(), py::capsule(permutation.data())); ; }); diff --git a/src/htool/distributed_operator/distributed_operator.hpp b/src/htool/distributed_operator/distributed_operator.hpp index ceccf29..08b6215 100644 --- a/src/htool/distributed_operator/distributed_operator.hpp +++ b/src/htool/distributed_operator/distributed_operator.hpp @@ -2,7 +2,6 @@ #define HTOOL_PYBIND11_DISTRIBUTED_OPERATOR_HPP #include "../misc/utility.hpp" #include "../misc/wrapper_mpi.hpp" -#include "./interface/partition.hpp" #include #include #include @@ -17,21 +16,24 @@ void declare_distributed_operator(py::module &m, const std::string &class_name) // Linear algebra py_class.def( - "__mul__", [](const Class &self, std::vector in) { - std::vector out(self.get_target_partition().get_global_size()); - self.vector_product_global_to_global(in.data(), out.data()); - return as_pyarray(std::move(out)); - ; + "__mul__", [](const Class &self, const py::array_t input) { + if (input.ndim() != 1) { + throw std::runtime_error("Wrong dimension for HMatrix-vector product"); // LCOV_EXCL_LINE + } + if (input.shape()[0] != self.get_source_partition().get_global_size()) { + throw std::runtime_error("Wrong size for HMatrix-vector product"); // LCOV_EXCL_LINE + } + py::array_t result(std::array{self.get_target_partition().get_global_size()}); + self.vector_product_global_to_global(input.data(), result.mutable_data()); + + return result; }, "in"_a); py_class.def( "__matmul__", [](const Class &self, py::array_t input) { int mu; - - if (input.ndim() == 1) { - mu = 1; - } else if (input.ndim() == 2) { + if (input.ndim() == 2) { mu = input.shape()[1]; } else { throw std::runtime_error("Wrong dimension for HMatrix-matrix product"); // LCOV_EXCL_LINE @@ -40,21 +42,12 @@ void declare_distributed_operator(py::module &m, const std::string &class_name) throw std::runtime_error("Wrong size for HMatrix-matrix product"); // LCOV_EXCL_LINE } - std::vector result(self.get_target_partition().get_global_size() * mu, 0); - - const CoefficientPrecision *in = input.data(); - int nc = input.shape()[0]; - const CoefficientPrecision *test = input.data(); + std::array shape{self.get_target_partition().get_global_size(), mu}; + py::array_t result(shape); - self.matrix_product_global_to_global(input.data(), result.data(), mu); + self.matrix_product_global_to_global(input.data(), result.mutable_data(), mu); - if (input.ndim() == 1) { - std::array shape{self.get_target_partition().get_global_size()}; - return py::array_t(shape, result.data()); - } else { - std::array shape{self.get_target_partition().get_global_size(), mu}; - return py::array_t(shape, result.data()); - } + return result; }, py::arg("input").noconvert(true)); } diff --git a/src/htool/distributed_operator/implementation/partition_from_cluster.hpp b/src/htool/distributed_operator/implementation/partition_from_cluster.hpp deleted file mode 100644 index 497790f..0000000 --- a/src/htool/distributed_operator/implementation/partition_from_cluster.hpp +++ /dev/null @@ -1,16 +0,0 @@ -#ifndef HTOOL_DISTRIBUTED_OPERATOR_PARTITION_FROM_CLUSTER_CPP -#define HTOOL_DISTRIBUTED_OPERATOR_PARTITION_FROM_CLUSTER_CPP - -#include -#include -#include - -template -void declare_partition_from_cluster(py::module &m, const std::string &class_name) { - using Class = htool::PartitionFromCluster; - - py::class_> py_class(m, class_name.c_str()); - py_class.def(py::init &>(), py::keep_alive<1, 2>()); -} - -#endif diff --git a/src/htool/distributed_operator/interface/partition.hpp b/src/htool/distributed_operator/interface/partition.hpp deleted file mode 100644 index ed42d39..0000000 --- a/src/htool/distributed_operator/interface/partition.hpp +++ /dev/null @@ -1,32 +0,0 @@ -#ifndef HTOOL_PYBIND11_DISTRIBUTED_OPERATOR_PARTITION_HPP -#define HTOOL_PYBIND11_DISTRIBUTED_OPERATOR_PARTITION_HPP - -#include -#include - -// class PyIPartition : public htool::IPartition { -// public: -// /* Inherit the constructors */ -// using htool::IPartition::IPartition; - -// /* Trampoline (need one for each virtual function) */ -// std::string go(int n_times) override { -// PYBIND11_OVERRIDE_PURE( -// std::string, /* Return type */ -// IPartition, /* Parent class */ -// go, /* Name of function in C++ (must match Python name) */ -// n_times /* Argument(s) */ -// ); -// } -// }; - -template -void declare_interface_partition(py::module &m, const std::string &class_name) { - using Class = htool::IPartition; - py::class_ py_class(m, class_name.c_str()); - py_class.def("get_size_of_partition", &Class::get_size_of_partition); - py_class.def("get_offset_of_partition", &Class::get_offset_of_partition); - py_class.def("get_global_size", &Class::get_global_size); -} - -#endif diff --git a/src/htool/distributed_operator/utility.hpp b/src/htool/distributed_operator/utility.hpp index 660c009..8b5f292 100644 --- a/src/htool/distributed_operator/utility.hpp +++ b/src/htool/distributed_operator/utility.hpp @@ -8,12 +8,23 @@ template void declare_distributed_operator_utility(py::module &m, std::string prefix = "") { - using DefaultApproximation = DefaultApproximationBuilder; - using LocalDefaultApproximation = DefaultLocalApproximationBuilder; - std::string default_approximation_name = prefix + "DefaultApproximationBuilder"; - std::string default_local_approximation_name = prefix + "DefaultLocalApproximationBuilder"; + using CustomApproximation = CustomApproximationBuilder; + using DefaultApproximation = DefaultApproximationBuilder; + using LocalDefaultApproximation = DefaultLocalApproximationBuilder; + using DistributedOperatorFromHMatrix = DistributedOperatorFromHMatrix; + + std::string custom_approximation_name = prefix + "CustomApproximationBuilder"; + std::string default_approximation_name = prefix + "DefaultApproximationBuilder"; + std::string default_local_approximation_name = prefix + "DefaultLocalApproximationBuilder"; + std::string distributed_operator_from_hmatrix_name = prefix + "DistributedOperatorFromHMatrix"; + + py::class_ custom_approximation_class(m, custom_approximation_name.c_str()); + custom_approximation_class.def(py::init &, const Cluster &, char, char, MPI_Comm_wrapper, const VirtualLocalOperator &>()); + custom_approximation_class.def_property_readonly( + "distributed_operator", [](const CustomApproximation &self) { return &self.distributed_operator; }, py::return_value_policy::reference_internal); + py::class_ default_approximation_class(m, default_approximation_name.c_str()); - default_approximation_class.def(py::init &, const Cluster &, const Cluster &, htool::underlying_type, htool::underlying_type, char, char, MPI_Comm_wrapper>()); + default_approximation_class.def(py::init &, const Cluster &, const Cluster &, htool::underlying_type, htool::underlying_type, 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( @@ -22,12 +33,21 @@ 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_ default_local_approximation_class(m, default_local_approximation_name.c_str()); - default_local_approximation_class.def(py::init &, const Cluster &, const Cluster &, htool::underlying_type, htool::underlying_type, char, char, MPI_Comm_wrapper>()); + default_local_approximation_class.def(py::init &, const Cluster &, const Cluster &, htool::underlying_type, htool::underlying_type, 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( "hmatrix", [](const LocalDefaultApproximation &self) { return &self.hmatrix; }, py::return_value_policy::reference_internal); default_local_approximation_class.def_property_readonly( "block_diagonal_hmatrix", [](const LocalDefaultApproximation &self) { return &*self.block_diagonal_hmatrix; }, py::return_value_policy::reference_internal); + + py::class_ distributed_operator_from_hmatrix_class(m, distributed_operator_from_hmatrix_name.c_str()); + distributed_operator_from_hmatrix_class.def(py::init &, const Cluster &, const Cluster &, const HMatrixTreeBuilder &, 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( + "hmatrix", [](const DistributedOperatorFromHMatrix &self) { return &self.hmatrix; }, py::return_value_policy::reference_internal); + distributed_operator_from_hmatrix_class.def_property_readonly( + "block_diagonal_hmatrix", [](const DistributedOperatorFromHMatrix &self) { return &*self.block_diagonal_hmatrix; }, py::return_value_policy::reference_internal); } #endif diff --git a/src/htool/hmatrix/hmatrix.hpp b/src/htool/hmatrix/hmatrix.hpp index 4dc9c3b..15c124f 100644 --- a/src/htool/hmatrix/hmatrix.hpp +++ b/src/htool/hmatrix/hmatrix.hpp @@ -22,22 +22,23 @@ void declare_HMatrix(py::module &m, const std::string &className) { using Class = HMatrix; py::class_ py_class(m, className.c_str()); - // py_class.def("build", [](Class &self, VirtualGeneratorCpp &mat, const py::array_t &x) { - // self.build(mat, x.data()); - // }); + py_class.def("to_dense", [](const Class &self) { + std::array shape{self.get_target_cluster().get_size(), self.get_source_cluster().get_size()}; + py::array_t dense(shape); + std::fill_n(dense.mutable_data(), dense.size(), CoefficientPrecision(0)); + copy_to_dense(self, dense.mutable_data()); + return dense; + }); - // py_class.def("build_dense_blocks", [](Class &self, VirtualDenseBlocksGeneratorCpp &dense_block_generator) { - // self.build_dense_blocks(dense_block_generator); - // }); + py_class.def("to_dense_in_user_numbering", [](const Class &self) { + std::array shape{self.get_target_cluster().get_size(), self.get_source_cluster().get_size()}; + py::array_t dense(shape); + std::fill_n(dense.mutable_data(), dense.size(), CoefficientPrecision(0)); + copy_to_dense_in_user_numbering(self, dense.mutable_data()); + return dense; + }); - // // Setters - // py_class.def("set_maxblocksize", &Class::set_maxblocksize); - // py_class.def("set_minsourcedepth", &Class::set_minsourcedepth); - // py_class.def("set_mintargetdepth", &Class::set_mintargetdepth); - // py_class.def("set_delay_dense_computation", &Class::set_delay_dense_computation); - // py_class.def("set_compression", [](Class &self, std::shared_ptr> mat) { - // self.set_compression(mat); - // }); + py_class.def("__deepcopy__", [](const Class &self, py::dict) { return Class(self); }, "memo"_a); // // Getters // py_class.def_property_readonly("shape", [](const Class &self) { @@ -82,13 +83,6 @@ void declare_HMatrix(py::module &m, const std::string &className) { // } // }); - // // Print information - // py_class.def("print_infos", &Class::print_infos); - // py_class.def("get_infos", overload_cast_()(&Class::get_infos, py::const_)); - // 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_sub_hmatrix", [](const HMatrix &hmatrix, const Cluster &target_cluster, const Cluster &source_cluster) { return &*hmatrix.get_sub_hmatrix(target_cluster, source_cluster); @@ -106,92 +100,6 @@ void declare_HMatrix(py::module &m, const std::string &className) { auto information = htool::get_distributed_hmatrix_information(hmatrix, comm); return information; }); - - // Plot pattern - py_class.def( - "display", [](const HMatrix &hmatrix, bool show = true) { - std::vector buf; - int nb_leaves = 0; - htool::preorder_tree_traversal( - hmatrix, - [&buf, &nb_leaves, &hmatrix](const HMatrix ¤t_hmatrix) { - if (current_hmatrix.is_leaf()) { - nb_leaves += 1; - buf.push_back(current_hmatrix.get_target_cluster().get_offset() - hmatrix.get_target_cluster().get_offset()); - buf.push_back(current_hmatrix.get_target_cluster().get_size()); - buf.push_back(current_hmatrix.get_source_cluster().get_offset() - hmatrix.get_source_cluster().get_offset()); - buf.push_back(current_hmatrix.get_source_cluster().get_size()); - buf.push_back(current_hmatrix.get_rank()); - } - }); - - // Import - py::object plt = py::module::import("matplotlib.pyplot"); - py::object patches = py::module::import("matplotlib.patches"); - py::object colors = py::module::import("matplotlib.colors"); - py::object numpy = py::module::import("numpy"); - - // First Data - int nr = hmatrix.get_target_cluster().get_size(); - int nc = hmatrix.get_source_cluster().get_size(); - py::array_t matrix({nr, nc}); - py::array_t mask_matrix({nr, nc}); - mask_matrix.attr("fill")(false); - - // Figure - py::tuple sublots_output = plt.attr("subplots")(1, 1); - py::object fig = sublots_output[0]; - py::object axes = sublots_output[1]; - // axes.attr() - - // Issue: there a shift of one pixel along the y-axis... - // int shift = axes.transData.transform([(0,0), (1,1)]) - // shift = shift[1,1] - shift[0,1] # 1 unit in display coords - int shift = 0; - - int max_rank = 0; - for (int p = 0; p < nb_leaves; p++) { - int i_row = buf[5 * p]; - int nb_row = buf[5 * p + 1]; - int i_col = buf[5 * p + 2]; - int nb_col = buf[5 * p + 3]; - int rank = buf[5 * p + 4]; - - if (rank > max_rank) { - max_rank = rank; - } - for (int i = 0; i < nb_row; i++) { - for (int j = 0; j < nb_col; j++) { - matrix.mutable_at(i_row + i, i_col + j) = rank; - if (rank == -1) { - mask_matrix.mutable_at(i_row + i, i_col + j) = true; - } - } - } - - py::object rect = patches.attr("Rectangle")(py::make_tuple(i_col - 0.5, i_row - 0.5 + shift), nb_col, nb_row, "linewidth"_a = 0.75, "edgecolor"_a = 'k', "facecolor"_a = "none"); - axes.attr("add_patch")(rect); - - if (rank >= 0 && nb_col / double(nc) > 0.05 && nb_row / double(nc) > 0.05) { - axes.attr("annotate")(rank, py::make_tuple(i_col + nb_col / 2., i_row + nb_row / 2.), "color"_a = "white", "size"_a = 10, "va"_a = "center", "ha"_a = "center"); - } - } - - // Colormap - py::object cmap = plt.attr("get_cmap")("YlGn"); - py::object new_cmap = colors.attr("LinearSegmentedColormap").attr("from_list")("trunc(YlGn,0.4,1)", cmap(numpy.attr("linspace")(0.4, 1, 100))); - - // Plot - py::object masked_matrix = numpy.attr("ma").attr("array")(matrix, "mask"_a = mask_matrix); - new_cmap.attr("set_bad")("color"_a = "red"); - - plt.attr("imshow")(masked_matrix, "cmap"_a = new_cmap, "vmin"_a = 0, "vmax"_a = 10); - plt.attr("draw")(); - if (show) { - plt.attr("show")(); // LCOV_EXCL_LINE - } - }, - py::arg("show") = true); } #endif diff --git a/src/htool/hmatrix/hmatrix_builder.hpp b/src/htool/hmatrix/hmatrix_builder.hpp index 75a1f6f..2a4578d 100644 --- a/src/htool/hmatrix/hmatrix_builder.hpp +++ b/src/htool/hmatrix/hmatrix_builder.hpp @@ -12,13 +12,13 @@ void declare_hmatrix_builder(py::module &m, const std::string &className) { py::class_ py_class(m, className.c_str()); // Constructor - py_class.def(py::init &, const htool::Cluster &, htool::underlying_type, CoordinatePrecision, char, char, int, int>()); + py_class.def(py::init &, const htool::Cluster &, htool::underlying_type, CoordinatePrecision, char, char, int, int, int>()); // Build - py_class.def("build", [](const Class &self, const VirtualGenerator &generator) { return self.build(generator); }); + // py_class.def("build", [](const Class &self, const VirtualGenerator &generator) { return self.build(generator); }); + py_class.def("build", [](const Class &self, const VirtualGeneratorInUserNumbering &generator) { return self.build(GeneratorWithPermutation(generator, self.get_target_cluster().get_permutation().data(), self.get_source_cluster().get_permutation().data())); }); // Setters - py_class.def("set_maximal_block_size", &Class::set_maximal_block_size); py_class.def("set_minimal_source_depth", &Class::set_minimal_source_depth); py_class.def("set_minimal_target_depth", &Class::set_minimal_target_depth); py_class.def("set_low_rank_generator", [](Class &self, const std::shared_ptr> &low_rank_generator) { self.set_low_rank_generator(low_rank_generator); }); diff --git a/src/htool/hmatrix/interfaces/virtual_generator.hpp b/src/htool/hmatrix/interfaces/virtual_generator.hpp index f96c5ee..8e4f816 100644 --- a/src/htool/hmatrix/interfaces/virtual_generator.hpp +++ b/src/htool/hmatrix/interfaces/virtual_generator.hpp @@ -7,46 +7,13 @@ using namespace htool; template -class VirtualGeneratorPython : public htool::VirtualGenerator { +class VirtualGeneratorInUserNumberingPython : public htool::VirtualGeneratorInUserNumbering { public: - using VirtualGenerator::VirtualGenerator; + using VirtualGeneratorInUserNumbering::VirtualGeneratorInUserNumbering; - void copy_submatrix(int M, int N, int row_offset, int col_offset, CoefficientPrecision *ptr) const override { + VirtualGeneratorInUserNumberingPython(const py::array_t &target_permutation, const py::array_t &source_permutation) : VirtualGeneratorInUserNumbering() {} - py::array_t mat(std::array{M, N}, ptr, py::capsule(ptr)); - build_submatrix(row_offset, col_offset, mat); - } - - // lcov does not see it because of trampoline I assume - virtual void build_submatrix(int row_offset, int col_offset, py::array_t &mat) const = 0; // LCOV_EXCL_LINE -}; - -template -class PyVirtualGenerator : public VirtualGeneratorPython { - public: - using VirtualGeneratorPython::VirtualGeneratorPython; - - /* Trampoline (need one for each virtual function) */ - virtual void build_submatrix(int row_offset, int col_offset, py::array_t &mat) const override { - PYBIND11_OVERRIDE_PURE( - void, /* Return type */ - VirtualGeneratorPython, /* Parent class */ - build_submatrix, /* Name of function in C++ (must match Python name) */ - row_offset, - col_offset, - mat /* Argument(s) */ - ); - } -}; - -template -class VirtualGeneratorWithPermutationPython : public htool::VirtualGeneratorWithPermutation { - public: - using VirtualGeneratorWithPermutation::VirtualGeneratorWithPermutation; - - VirtualGeneratorWithPermutationPython(const py::array_t &target_permutation, const py::array_t &source_permutation) : VirtualGeneratorWithPermutation(target_permutation.data(), source_permutation.data()) {} - - void copy_submatrix_from_user_numbering(int M, int N, const int *const rows, const int *const cols, CoefficientPrecision *ptr) const override { + void copy_submatrix(int M, int N, const int *const rows, const int *const cols, CoefficientPrecision *ptr) const override { if (M * N > 0) { py::array_t mat(std::array{M, N}, ptr, py::capsule(ptr)); @@ -62,15 +29,15 @@ class VirtualGeneratorWithPermutationPython : public htool::VirtualGeneratorWith }; template -class PyVirtualGeneratorWithPermutation : public VirtualGeneratorWithPermutationPython { +class PyVirtualGeneratorInUserNumbering : public VirtualGeneratorInUserNumberingPython { public: - using VirtualGeneratorWithPermutationPython::VirtualGeneratorWithPermutationPython; + using VirtualGeneratorInUserNumberingPython::VirtualGeneratorInUserNumberingPython; /* Trampoline (need one for each virtual function) */ virtual void build_submatrix(const py::array_t &J, const py::array_t &K, py::array_t &mat) const override { PYBIND11_OVERRIDE_PURE( void, /* Return type */ - VirtualGeneratorWithPermutationPython, /* Parent class */ + VirtualGeneratorInUserNumberingPython, /* Parent class */ build_submatrix, /* Name of function in C++ (must match Python name) */ J, K, @@ -81,21 +48,13 @@ class PyVirtualGeneratorWithPermutation : public VirtualGeneratorWithPermutation template void declare_virtual_generator(py::module &m, const std::string &className, const std::string &base_class_name) { - using BaseClass = VirtualGenerator; + using BaseClass = VirtualGeneratorInUserNumbering; py::class_(m, base_class_name.c_str()); - using Class = VirtualGeneratorPython; - py::class_> py_class(m, className.c_str()); + using Class = VirtualGeneratorInUserNumberingPython; + py::class_> py_class(m, className.c_str()); py_class.def(py::init<>()); py_class.def("build_submatrix", &Class::build_submatrix); - - using BaseClassWithPermutation = VirtualGeneratorWithPermutation; - py::class_(m, (base_class_name + "WithPermutation").c_str()); - - using ClassWithPermutation = VirtualGeneratorWithPermutationPython; - py::class_> py_class_with_permutation(m, (className + "WithPermutation").c_str()); - py_class_with_permutation.def(py::init &, const py::array_t &>()); - py_class_with_permutation.def("build_submatrix", &ClassWithPermutation::build_submatrix); } #endif diff --git a/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp b/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp index 18fa23b..dcfdc8c 100644 --- a/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp +++ b/src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp @@ -16,14 +16,18 @@ class VirtualLowRankGeneratorPython : public VirtualLowRankGenerator> m_mats_U; // owned by Python mutable std::vector> m_mats_V; // owned by Python - - // py::array_t m_mat_U, m_mat_V; + const VirtualGeneratorInUserNumbering &m_generator_in_user_numbering; public: using VirtualLowRankGenerator::VirtualLowRankGenerator; + VirtualLowRankGeneratorPython(const VirtualGeneratorInUserNumbering &generator_in_user_numbering) : m_generator_in_user_numbering(generator_in_user_numbering) {} + void copy_low_rank_approximation(const VirtualGenerator &A, const Cluster &target_cluster, const Cluster &source_cluster, underlying_type epsilon, int &rank, Matrix &U, Matrix &V) const override { - build_low_rank_approximation(A, target_cluster.get_size(), source_cluster.get_size(), target_cluster.get_offset(), source_cluster.get_offset(), epsilon); + py::array_t rows(target_cluster.get_size(), target_cluster.get_permutation().data() + target_cluster.get_offset(), py::capsule(target_cluster.get_permutation().data())); + py::array_t cols(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); @@ -32,13 +36,17 @@ class VirtualLowRankGeneratorPython : public VirtualLowRankGenerator &A, int target_size, int source_size, int target_offset, int source_offset, underlying_type epsilon) const = 0; // LCOV_EXCL_LINE + virtual void build_low_rank_approximation(const py::array_t &rows, const py::array_t &cols, underlying_type epsilon) const = 0; // LCOV_EXCL_LINE void set_U(py::array_t U0) { m_mats_U.push_back(U0); // no copy here } void set_V(py::array_t V0) { m_mats_V.push_back(V0); } void set_rank(int rank) { m_rank = rank; } + + void build_submatrix(const py::array_t rows, const py::array_t cols, py::array_t &mat) const { + m_generator_in_user_numbering.copy_submatrix(rows.size(), cols.size(), rows.data(), cols.data(), mat.mutable_data()); + } }; template @@ -47,46 +55,27 @@ class PyVirtualLowRankGenerator : public VirtualLowRankGeneratorPython::VirtualLowRankGeneratorPython; /* Trampoline (need one for each virtual function) */ - virtual void build_low_rank_approximation(const VirtualGenerator &A, int target_size, int source_size, int target_offset, int source_offset, underlying_type epsilon) const override { + virtual void build_low_rank_approximation(const py::array_t &rows, const py::array_t &cols, underlying_type epsilon) const override { PYBIND11_OVERRIDE_PURE( void, /* Return type */ VirtualLowRankGeneratorPython, /* Parent class */ build_low_rank_approximation, /* Name of function in C++ (must match Python name) */ - A, - target_size, - source_size, - target_offset, - source_offset, + rows, + cols, epsilon); } }; template void declare_custom_VirtualLowRankGenerator(py::module &m, const std::string &className) { - // using BaseClass = VirtualLowRankGenerator; - // py::class_>(m, base_class_name.c_str()); - using Class = VirtualLowRankGeneratorPython; py::class_, PyVirtualLowRankGenerator> py_class(m, className.c_str()); - py_class.def(py::init<>()); + 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); } -// template -// void declare_VirtualLowRankGenerator(py::module &m, const std::string &className) { -// using Class = VirtualLowRankGenerator; -// py::class_ py_class(m, className.c_str()); -// } - -// template