From e37d075f324ccf9b0495fc9c0d1d90e475cf6e59 Mon Sep 17 00:00:00 2001 From: Pierre Marchand Date: Sat, 28 Mar 2026 20:04:02 +0100 Subject: [PATCH] add hmat/mat product, h-lu and h-cholesky + ruff clean --- .github/workflows/CI.yml | 6 +- .../advanced/define_custom_local_operator.py | 3 +- .../define_custom_low_rank_generator.py | 3 +- .../use_cluster_with_given_partition.py | 3 +- .../use_custom_dense_block_generator.py | 3 +- .../use_custom_global_to_local_operator.py | 3 +- .../use_custom_local_to_local_operator.py | 3 +- .../use_custom_low_rank_approximation.py | 3 +- .../advanced/use_local_hmatrix_compression.py | 3 +- example/define_generators.py | 3 +- example/use_cluster.py | 3 +- example/use_ddm_solver.py | 3 +- example/use_distributed_operator.py | 3 +- example/use_hmatrix.py | 15 +++- lib/htool | 2 +- pyproject.toml | 2 + src/htool/hmatrix/hmatrix.hpp | 63 ++++++++++++++- src/htool/hmatrix/hmatrix_tree_builder.hpp | 1 + src/htool/main.cpp | 1 + tests/conftest.py | 2 +- tests/test_cluster.py | 3 +- tests/test_ddm_solver.py | 3 +- tests/test_distributed_operator.py | 3 +- tests/test_hmatrix.py | 76 ++++++++++++++++--- 24 files changed, 179 insertions(+), 34 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 631319f..184e60d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -161,8 +161,10 @@ jobs: git diff --exit-code - name: Check with ruff - run: | - ruff check example/ tests/ + uses: astral-sh/ruff-action@v3 + with: + version: "0.15.8" + args: "check ./example ./tests" - name: Generate coverage reports if: matrix.CODE_COVERAGE == 'ON' diff --git a/example/advanced/define_custom_local_operator.py b/example/advanced/define_custom_local_operator.py index 73c8414..eb80a43 100644 --- a/example/advanced/define_custom_local_operator.py +++ b/example/advanced/define_custom_local_operator.py @@ -1,6 +1,7 @@ -import Htool import numpy as np +import Htool + class CustomRestrictedGlobalToLocalOperator(Htool.RestrictedGlobalToLocalOperator): def __init__( diff --git a/example/advanced/define_custom_low_rank_generator.py b/example/advanced/define_custom_low_rank_generator.py index 917c4f8..f9ed9e7 100644 --- a/example/advanced/define_custom_low_rank_generator.py +++ b/example/advanced/define_custom_low_rank_generator.py @@ -1,8 +1,9 @@ import math -import Htool import numpy as np +import Htool + class CustomSVD(Htool.VirtualLowRankGenerator): def __init__(self, generator: Htool.VirtualGenerator, allow_copy: bool = True): diff --git a/example/advanced/use_cluster_with_given_partition.py b/example/advanced/use_cluster_with_given_partition.py index 4914fc9..986f0fe 100644 --- a/example/advanced/use_cluster_with_given_partition.py +++ b/example/advanced/use_cluster_with_given_partition.py @@ -1,11 +1,12 @@ import os import sys -import Htool import matplotlib.pyplot as plt import mpi4py import numpy as np +import Htool + # Add the parent directory to sys.path sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) diff --git a/example/advanced/use_custom_dense_block_generator.py b/example/advanced/use_custom_dense_block_generator.py index 651af50..c886296 100644 --- a/example/advanced/use_custom_dense_block_generator.py +++ b/example/advanced/use_custom_dense_block_generator.py @@ -2,12 +2,13 @@ import os import sys -import Htool import mpi4py import numpy as np from define_custom_dense_blocks_generator import CustomDenseBlocksGenerator from matplotlib import pyplot as plt +import Htool + # Add the parent directory to sys.path sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) diff --git a/example/advanced/use_custom_global_to_local_operator.py b/example/advanced/use_custom_global_to_local_operator.py index 9300005..c64a1a1 100644 --- a/example/advanced/use_custom_global_to_local_operator.py +++ b/example/advanced/use_custom_global_to_local_operator.py @@ -1,11 +1,12 @@ import os import sys -import Htool import mpi4py import numpy as np from define_custom_local_operator import CustomRestrictedGlobalToLocalOperator +import Htool + # Add the parent directory to sys.path sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) from create_geometry import create_partitionned_geometries diff --git a/example/advanced/use_custom_local_to_local_operator.py b/example/advanced/use_custom_local_to_local_operator.py index d61f6b1..d651eea 100644 --- a/example/advanced/use_custom_local_to_local_operator.py +++ b/example/advanced/use_custom_local_to_local_operator.py @@ -1,11 +1,12 @@ import os import sys -import Htool import mpi4py import numpy as np from define_custom_local_operator import CustomRestrictedGlobalToLocalOperator +import Htool + # Add the parent directory to sys.path sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) from create_geometry import create_partitionned_geometries diff --git a/example/advanced/use_custom_low_rank_approximation.py b/example/advanced/use_custom_low_rank_approximation.py index 7e8f2c2..3d1a6b6 100644 --- a/example/advanced/use_custom_low_rank_approximation.py +++ b/example/advanced/use_custom_low_rank_approximation.py @@ -2,12 +2,13 @@ import os import sys -import Htool import mpi4py import numpy as np from define_custom_low_rank_generator import CustomSVD from matplotlib import pyplot as plt +import Htool + # Add the parent directory to sys.path sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))) from create_geometry import create_partitionned_geometries diff --git a/example/advanced/use_local_hmatrix_compression.py b/example/advanced/use_local_hmatrix_compression.py index aefbfca..5f487c9 100644 --- a/example/advanced/use_local_hmatrix_compression.py +++ b/example/advanced/use_local_hmatrix_compression.py @@ -1,6 +1,5 @@ import logging -import Htool import matplotlib.pyplot as plt import mpi4py import numpy as np @@ -8,6 +7,8 @@ from define_custom_local_operator import CustomGlobalToLocalOperator from define_generators import CustomGenerator +import Htool + logging.basicConfig(level=logging.INFO) # Random geometry diff --git a/example/define_generators.py b/example/define_generators.py index 0703093..9349cce 100644 --- a/example/define_generators.py +++ b/example/define_generators.py @@ -1,6 +1,7 @@ -import Htool import numpy as np +import Htool + class CustomGenerator(Htool.VirtualGenerator): def __init__(self, target_points, source_points): diff --git a/example/use_cluster.py b/example/use_cluster.py index e227065..6336f73 100644 --- a/example/use_cluster.py +++ b/example/use_cluster.py @@ -1,7 +1,8 @@ -import Htool import matplotlib.pyplot as plt from create_geometry import create_random_geometries +import Htool + # Random geometry nb_rows = 500 nb_cols = 500 diff --git a/example/use_ddm_solver.py b/example/use_ddm_solver.py index 68f8c41..194606a 100644 --- a/example/use_ddm_solver.py +++ b/example/use_ddm_solver.py @@ -1,13 +1,14 @@ import copy import logging -import Htool import matplotlib.pyplot as plt import mpi4py import numpy as np from create_geometry import create_random_geometries from define_generators import CustomGenerator +import Htool + logging.basicConfig(level=logging.INFO) # Random geometry diff --git a/example/use_distributed_operator.py b/example/use_distributed_operator.py index f8f957f..9e5a05e 100644 --- a/example/use_distributed_operator.py +++ b/example/use_distributed_operator.py @@ -1,12 +1,13 @@ import logging -import Htool import matplotlib.pyplot as plt import mpi4py import numpy as np from create_geometry import create_partitionned_geometries from define_generators import CustomGenerator +import Htool + logging.basicConfig(level=logging.INFO) # Random geometry diff --git a/example/use_hmatrix.py b/example/use_hmatrix.py index 46bf32c..95231fd 100644 --- a/example/use_hmatrix.py +++ b/example/use_hmatrix.py @@ -1,11 +1,12 @@ import logging -import Htool import matplotlib.pyplot as plt import numpy as np from create_geometry import create_random_points_in_disk, create_random_points_in_sphere from define_generators import CustomGenerator +import Htool + logging.basicConfig(level=logging.INFO) # Random geometry @@ -44,11 +45,17 @@ ) # HMatrix vector product -dense_in_user_numbering = hmatrix.to_dense_in_user_numbering() np.random.seed(0) x = np.random.rand(size) -y = generator.mat_vec(x) -y_dense = dense_in_user_numbering.dot(x) +y_dense = generator.mat_vec(x) +y = hmatrix * x +print(np.linalg.norm(y - y_dense) / np.linalg.norm(y_dense), epsilon) + +# HMatrix matrix product +np.random.seed(0) +x = np.random.rand(size, 2) +y_dense = generator.mat_mat(x) +y = hmatrix @ x print(np.linalg.norm(y - y_dense) / np.linalg.norm(y_dense), epsilon) diff --git a/lib/htool b/lib/htool index d7c0fa8..f0b1541 160000 --- a/lib/htool +++ b/lib/htool @@ -1 +1 @@ -Subproject commit d7c0fa8b42c461446b92a4891dea78532eeea6b5 +Subproject commit f0b154130a0f542b19ad6d3b4e6986539fb702b7 diff --git a/pyproject.toml b/pyproject.toml index 63721da..26fc408 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -47,6 +47,8 @@ ignore = [] fixable = ["ALL"] unfixable = [] +[tool.ruff.lint.isort] +known-first-party = ["Htool"] [build-system] requires = [ diff --git a/src/htool/hmatrix/hmatrix.hpp b/src/htool/hmatrix/hmatrix.hpp index 2d626b3..27d7752 100644 --- a/src/htool/hmatrix/hmatrix.hpp +++ b/src/htool/hmatrix/hmatrix.hpp @@ -10,7 +10,9 @@ #include #include #include +#include #include +#include #ifdef HAVE_MPI # include "../misc/wrapper_mpi.hpp" @@ -53,6 +55,44 @@ void declare_HMatrix(py::module &m, const std::string &className) { py_class.def("get_target_cluster", &HMatrix::get_target_cluster, py::return_value_policy::reference_internal); py_class.def("get_source_cluster", &HMatrix::get_source_cluster, py::return_value_policy::reference_internal); + py_class.def("lu_factorization", [](HMatrix &hmatrix) { + htool::lu_factorization(hmatrix); + }); + py_class.def("cholesky_factorization", [](HMatrix &hmatrix, char UPLO) { + htool::cholesky_factorization(UPLO, hmatrix); + }); + py_class.def("lu_solve", [](const Class &self, char trans, const py::array_t &input) { + std::vector shape; + if (input.ndim() == 1) { + shape = {input.shape()[0]}; + } else if (input.ndim() == 2) { + shape = {input.shape()[0], input.shape()[1]}; + } else { + throw std::runtime_error("Wrong dimension for HMatrix-LU input"); // LCOV_EXCL_LINE + } + py::array_t result(shape); + std::copy_n(input.data(), input.size(), result.mutable_data()); + htool::MatrixView output_view(result.shape()[0], input.ndim() == 1 ? 1 : result.shape()[1], result.mutable_data()); + htool::lu_solve(trans, self, output_view); + return result; + }); + + py_class.def("cholesky_solve", [](const Class &self, char UPLO, const py::array_t &input) { + std::vector shape; + if (input.ndim() == 1) { + shape = {input.shape()[0]}; + } else if (input.ndim() == 2) { + shape = {input.shape()[0], input.shape()[1]}; + } else { + throw std::runtime_error("Wrong dimension for HMatrix-Cholesky input"); // LCOV_EXCL_LINE + } + py::array_t result(shape); + std::copy_n(input.data(), input.size(), result.mutable_data()); + htool::MatrixView output_view(result.shape()[0], input.ndim() == 1 ? 1 : result.shape()[1], result.mutable_data()); + htool::cholesky_solve(UPLO, self, output_view); + return result; + }); + m.def("recompression", &htool::recompression &)>>); m.def("recompression", [](HMatrix &hmatrix) { recompression(hmatrix); }); m.def("openmp_recompression", &htool::openmp_recompression &)>>); @@ -69,14 +109,33 @@ void declare_HMatrix(py::module &m, const std::string &className) { py::array_t result(self.get_target_cluster().get_size()); std::fill_n(result.mutable_data(), self.get_target_cluster().get_size(), CoefficientPrecision(0)); - htool::Matrix dense_mat(self.get_target_cluster().get_size(), self.get_source_cluster().get_size()); - copy_to_dense_in_user_numbering(self, dense_mat.data()); char trans = 'N'; htool::add_hmatrix_vector_product(trans, CoefficientPrecision(1), self, input.data(), CoefficientPrecision(0), result.mutable_data()); return result; }, "in"_a); + + py_class.def( + "__matmul__", [](const Class &self, const py::array_t input) { + if (input.ndim() != 2) { + throw std::runtime_error("Wrong dimension for HMatrix-matrix product"); // LCOV_EXCL_LINE + } + if (input.shape()[0] != self.get_source_cluster().get_size()) { + throw std::runtime_error("Wrong size for HMatrix-matrix product"); // LCOV_EXCL_LINE + } + py::array_t result({input.shape()[0], input.shape()[1]}); + std::fill_n(result.mutable_data(), input.shape()[0] * input.shape()[1], CoefficientPrecision(0)); + + htool::MatrixView input_view(input.shape()[0], input.shape()[1], input.data()); + htool::MatrixView output_view(input.shape()[0], input.shape()[1], result.mutable_data()); + char transa = 'N'; + char transb = 'N'; + htool::add_hmatrix_matrix_product(transa, transb, CoefficientPrecision(1), self, input_view, CoefficientPrecision(0), output_view); + + return result; + }, + "in"_a); } #endif diff --git a/src/htool/hmatrix/hmatrix_tree_builder.hpp b/src/htool/hmatrix/hmatrix_tree_builder.hpp index a77567a..217fe28 100644 --- a/src/htool/hmatrix/hmatrix_tree_builder.hpp +++ b/src/htool/hmatrix/hmatrix_tree_builder.hpp @@ -40,5 +40,6 @@ void declare_hmatrix_builder(py::module &m, const std::string &className) { py_class.def("set_minimal_target_depth", &Class::set_minimal_target_depth); py_class.def("set_low_rank_generator", [](Class &self, std::shared_ptr> low_rank_generator) { self.set_low_rank_generator(low_rank_generator); }); py_class.def("set_dense_blocks_generator", [](Class &self, std::shared_ptr> dense_blocks_generator) { self.set_dense_blocks_generator(dense_blocks_generator); }); + py_class.def("set_block_tree_consistency", [](Class &self, bool consistency) { self.set_block_tree_consistency(consistency); }); } #endif diff --git a/src/htool/main.cpp b/src/htool/main.cpp index 89aee6e..548445d 100644 --- a/src/htool/main.cpp +++ b/src/htool/main.cpp @@ -84,6 +84,7 @@ PYBIND11_MODULE(Htool, m) { declare_matplotlib_cluster(m); declare_matplotlib_hmatrix(m); + declare_matplotlib_hmatrix, double>(m); declare_virtual_partitioning>(m, "Complex"); declare_LowRankMatrix>(m, "ComplexLowRankMatrix"); diff --git a/tests/conftest.py b/tests/conftest.py index ed640fb..23db6d4 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -2,11 +2,11 @@ import pathlib import struct -import Htool import mpi4py import numpy as np import pytest +import Htool from example.advanced.define_custom_dense_blocks_generator import ( CustomDenseBlocksGenerator, ) diff --git a/tests/test_cluster.py b/tests/test_cluster.py index 6161878..b701180 100644 --- a/tests/test_cluster.py +++ b/tests/test_cluster.py @@ -1,8 +1,9 @@ -import Htool import matplotlib.pyplot as plt import mpi4py import pytest +import Htool + @pytest.mark.parametrize( "dimension,nb_rows,nb_cols,symmetry,partition_type,number_of_children", diff --git a/tests/test_ddm_solver.py b/tests/test_ddm_solver.py index 3b9a93c..09938ce 100644 --- a/tests/test_ddm_solver.py +++ b/tests/test_ddm_solver.py @@ -1,7 +1,6 @@ import copy import logging -import Htool import mpi4py import numpy as np import pytest @@ -9,6 +8,8 @@ from scipy.linalg import eig, eigh, ldl, solve_triangular from scipy.sparse.linalg import LinearOperator, eigsh +import Htool + class CustomDenseGeneoCoarseSpaceDenseBuilder( Htool.VirtualGeneoCoarseSpaceDenseBuilder diff --git a/tests/test_distributed_operator.py b/tests/test_distributed_operator.py index 2fc4b05..30de09d 100644 --- a/tests/test_distributed_operator.py +++ b/tests/test_distributed_operator.py @@ -1,9 +1,10 @@ -import Htool import matplotlib.pyplot as plt import mpi4py import numpy as np import pytest +import Htool + @pytest.mark.parametrize("epsilon", [1e-3, 1e-6]) @pytest.mark.parametrize("eta", [10]) diff --git a/tests/test_hmatrix.py b/tests/test_hmatrix.py index b495bcf..7eed45a 100644 --- a/tests/test_hmatrix.py +++ b/tests/test_hmatrix.py @@ -1,20 +1,27 @@ import copy import logging -import Htool import numpy as np import pytest +import Htool from example.advanced.define_custom_low_rank_generator import CustomSVD from example.create_geometry import create_random_geometries from example.define_generators import CustomGenerator @pytest.mark.parametrize( - "loglevel", - [logging.INFO, logging.DEBUG, logging.WARNING, logging.ERROR, logging.CRITICAL], + "loglevel,symmetry", + [ + (logging.INFO, "N"), + (logging.DEBUG, "N"), + (logging.WARNING, "N"), + (logging.ERROR, "N"), + (logging.CRITICAL, "N"), + (logging.INFO, "S"), + ], ) -def test_hmatrix(loglevel): +def test_hmatrix(loglevel, symmetry): logging.basicConfig(level=loglevel) # Random geometry @@ -37,12 +44,18 @@ def test_hmatrix(loglevel): target_points, number_of_children, ) - source_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( - source_points, number_of_children - ) + if symmetry == "N": + source_cluster: Htool.Cluster = cluster_builder.create_cluster_tree( + source_points, number_of_children + ) + else: + source_cluster = target_cluster # Build generator - generator = CustomGenerator(target_points, source_points) + if symmetry == "N": + generator = CustomGenerator(target_points, source_points) + else: + generator = CustomGenerator(target_points, target_points) # Custom low rank generator low_rank_generator = CustomSVD(generator, False) @@ -64,11 +77,56 @@ def test_hmatrix(loglevel): np.random.seed(0) x = np.random.rand(nb_cols) y = hmatrix * x + y_exact = generator.mat_vec(x) y_dense = dense_in_user_numbering.dot(x) y_copy = copy_hmatrix * x - assert np.linalg.norm(y - y_dense) / np.linalg.norm(y_dense) < epsilon + assert np.linalg.norm(y - y_exact) / np.linalg.norm(y_exact) < epsilon + assert np.linalg.norm(y - y_dense) / np.linalg.norm(y_dense) < 1e-10 + assert np.linalg.norm(y - y_copy) < 1e-10 + + # HMatrix matrix product + np.random.seed(0) + x = np.random.rand(nb_cols, 2) + y = hmatrix @ x + y_exact = generator.mat_mat(x) + y_dense = dense_in_user_numbering @ x + y_copy = copy_hmatrix @ x + assert np.linalg.norm(y - y_exact) / np.linalg.norm(y_exact) < epsilon + assert np.linalg.norm(y - y_dense) / np.linalg.norm(y_dense) < 1e-10 assert np.linalg.norm(y - y_copy) < 1e-10 + if symmetry != "N": + # HLU factorization + copy_hmatrix.lu_factorization() + + # HLU solve vec + x_ref = np.ones(nb_cols) + y = hmatrix * x_ref + x = copy_hmatrix.lu_solve("N", y) + assert np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref) < epsilon + + # HLU solve mat + x_ref = np.ones((nb_cols, 2)) + y = hmatrix @ x_ref + x = copy_hmatrix.lu_solve("N", y) + assert np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref) < epsilon + + # Cholesky factorization + copy_hmatrix = copy.deepcopy(hmatrix) + copy_hmatrix.cholesky_factorization("L") + + # Cholesky solve vec + x_ref = np.ones(nb_cols) + y = hmatrix * x_ref + x = copy_hmatrix.cholesky_solve("L", y) + assert np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref) < epsilon + + # Cholesky solve mat + x_ref = np.ones((nb_cols, 2)) + y = hmatrix @ x_ref + x = copy_hmatrix.cholesky_solve("L", y) + assert np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref) < epsilon + # Output hmatrix_tree_parameter = hmatrix.get_tree_parameters() hmatrix_local_information = hmatrix.get_local_information()