add hmat/mat product, h-lu and h-cholesky + ruff clean

.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'

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__(

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):

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__), "..")))
 

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__), "..")))
 

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

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

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

example/advanced/use_local_hmatrix_compression.py

@@ -1,13 +1,14 @@
 import logging
 
-import Htool
 import matplotlib.pyplot as plt
 import mpi4py
 import numpy as np
 from create_geometry import create_random_geometries
 from define_custom_local_operator import CustomGlobalToLocalOperator
 from define_generators import CustomGenerator
 
+import Htool
+
 logging.basicConfig(level=logging.INFO)
 
 # Random geometry

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):

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

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

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

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)
 
 

pyproject.toml

@@ -47,6 +47,8 @@ ignore = []
 fixable = ["ALL"]
 unfixable = []
 
+[tool.ruff.lint.isort]
+known-first-party = ["Htool"]
 
 [build-system]
 requires = [

src/htool/hmatrix/hmatrix.hpp

@@ -10,7 +10,9 @@
 #include <htool/hmatrix/hmatrix.hpp>
 #include <htool/hmatrix/hmatrix_output.hpp>
 #include <htool/hmatrix/linalg/add_hmatrix_vector_product.hpp>
+#include <htool/hmatrix/linalg/factorization.hpp>
 #include <htool/hmatrix/utils/recompression.hpp>
+#include <htool/matrix/matrix_view.hpp>
 
 #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<CoefficientPrecision, CoordinatePrecision>::get_target_cluster, py::return_value_policy::reference_internal);
     py_class.def("get_source_cluster", &HMatrix<CoefficientPrecision, CoordinatePrecision>::get_source_cluster, py::return_value_policy::reference_internal);
 
+    py_class.def("lu_factorization", [](HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix) {
+        htool::lu_factorization(hmatrix);
+    });
+    py_class.def("cholesky_factorization", [](HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix, char UPLO) {
+        htool::cholesky_factorization(UPLO, hmatrix);
+    });
+    py_class.def("lu_solve", [](const Class &self, char trans, const py::array_t<CoefficientPrecision, py::array::f_style> &input) {
+        std::vector<ssize_t> 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<CoefficientPrecision, py::array::f_style> result(shape);
+        std::copy_n(input.data(), input.size(), result.mutable_data());
+        htool::MatrixView<CoefficientPrecision> 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<CoefficientPrecision, py::array::f_style> &input) {
+        std::vector<ssize_t> 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<CoefficientPrecision, py::array::f_style> result(shape);
+        std::copy_n(input.data(), input.size(), result.mutable_data());
+        htool::MatrixView<CoefficientPrecision> 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<CoefficientPrecision, CoordinatePrecision, std::function<void(LowRankMatrix<CoefficientPrecision> &)>>);
     m.def("recompression", [](HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix) { recompression(hmatrix); });
     m.def("openmp_recompression", &htool::openmp_recompression<CoefficientPrecision, CoordinatePrecision, std::function<void(LowRankMatrix<CoefficientPrecision> &)>>);
@@ -69,14 +109,33 @@ void declare_HMatrix(py::module &m, const std::string &className) {
             py::array_t<CoefficientPrecision, py::array::f_style> result(self.get_target_cluster().get_size());
             std::fill_n(result.mutable_data(), self.get_target_cluster().get_size(), CoefficientPrecision(0));
 
-            htool::Matrix<CoefficientPrecision> 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<CoefficientPrecision, py::array::f_style> 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<CoefficientPrecision, py::array::f_style> result({input.shape()[0], input.shape()[1]});
+            std::fill_n(result.mutable_data(), input.shape()[0] * input.shape()[1], CoefficientPrecision(0));
+
+            htool::MatrixView<const CoefficientPrecision> input_view(input.shape()[0], input.shape()[1], input.data());
+            htool::MatrixView<CoefficientPrecision> 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

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<VirtualLowRankGeneratorPython<CoefficientPrecision>> low_rank_generator) { self.set_low_rank_generator(low_rank_generator); });
     py_class.def("set_dense_blocks_generator", [](Class &self, std::shared_ptr<VirtualDenseBlocksGeneratorPython<CoefficientPrecision>> 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

src/htool/main.cpp

@@ -84,6 +84,7 @@ PYBIND11_MODULE(Htool, m) {
 
     declare_matplotlib_cluster<double>(m);
     declare_matplotlib_hmatrix<double, double>(m);
+    declare_matplotlib_hmatrix<std::complex<double>, double>(m);
 
     declare_virtual_partitioning<std::complex<double>>(m, "Complex");
     declare_LowRankMatrix<std::complex<double>>(m, "ComplexLowRankMatrix");

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,
 )

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",

tests/test_ddm_solver.py

@@ -1,14 +1,15 @@
 import copy
 import logging
 
-import Htool
 import mpi4py
 import numpy as np
 import pytest
 from conftest import ComplexGeneratorFromMatrix, GeneratorFromMatrix
 from scipy.linalg import eig, eigh, ldl, solve_triangular
 from scipy.sparse.linalg import LinearOperator, eigsh
 
+import Htool
+
 
 class CustomDenseGeneoCoarseSpaceDenseBuilder(
     Htool.VirtualGeneoCoarseSpaceDenseBuilder

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])