Skip to content

Commit a91ef8f

Browse files
PierreMarchand20PierreMarchand20
committed
add hmat/mat product, h-lu and h-cholesky + ruff clean
1 parent 18eb16c commit a91ef8f

14 files changed

Lines changed: 157 additions & 24 deletions

‎example/define_generators.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
1-
import Htool
21
import numpy as np
32

3+
import Htool
4+
45

56
class CustomGenerator(Htool.VirtualGenerator):
67
def __init__(self, target_points, source_points):

‎example/use_cluster.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,8 @@
1-
import Htool
21
import matplotlib.pyplot as plt
32
from create_geometry import create_random_geometries
43

4+
import Htool
5+
56
# Random geometry
67
nb_rows = 500
78
nb_cols = 500

‎example/use_ddm_solver.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,14 @@
11
import copy
22
import logging
33

4-
import Htool
54
import matplotlib.pyplot as plt
65
import mpi4py
76
import numpy as np
87
from create_geometry import create_random_geometries
98
from define_generators import CustomGenerator
109

10+
import Htool
11+
1112
logging.basicConfig(level=logging.INFO)
1213

1314
# Random geometry

‎example/use_distributed_operator.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,13 @@
11
import logging
22

3-
import Htool
43
import matplotlib.pyplot as plt
54
import mpi4py
65
import numpy as np
76
from create_geometry import create_partitionned_geometries
87
from define_generators import CustomGenerator
98

9+
import Htool
10+
1011
logging.basicConfig(level=logging.INFO)
1112

1213
# Random geometry

‎example/use_hmatrix.py‎

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,12 @@
11
import logging
22

3-
import Htool
43
import matplotlib.pyplot as plt
54
import numpy as np
65
from create_geometry import create_random_points_in_disk, create_random_points_in_sphere
76
from define_generators import CustomGenerator
87

8+
import Htool
9+
910
logging.basicConfig(level=logging.INFO)
1011

1112
# Random geometry
@@ -44,11 +45,17 @@
4445
)
4546

4647
# HMatrix vector product
47-
dense_in_user_numbering = hmatrix.to_dense_in_user_numbering()
4848
np.random.seed(0)
4949
x = np.random.rand(size)
50-
y = generator.mat_vec(x)
51-
y_dense = dense_in_user_numbering.dot(x)
50+
y_dense = generator.mat_vec(x)
51+
y = hmatrix * x
52+
print(np.linalg.norm(y - y_dense) / np.linalg.norm(y_dense), epsilon)
53+
54+
# HMatrix matrix product
55+
np.random.seed(0)
56+
x = np.random.rand(size, 2)
57+
y_dense = generator.mat_mat(x)
58+
y = hmatrix @ x
5259
print(np.linalg.norm(y - y_dense) / np.linalg.norm(y_dense), epsilon)
5360

5461

‎lib/htool‎

Submodule htool updated 35 files

‎src/htool/hmatrix/hmatrix.hpp‎

Lines changed: 61 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,9 @@
1010
#include <htool/hmatrix/hmatrix.hpp>
1111
#include <htool/hmatrix/hmatrix_output.hpp>
1212
#include <htool/hmatrix/linalg/add_hmatrix_vector_product.hpp>
13+
#include <htool/hmatrix/linalg/factorization.hpp>
1314
#include <htool/hmatrix/utils/recompression.hpp>
15+
#include <htool/matrix/matrix_view.hpp>
1416

1517
#ifdef HAVE_MPI
1618
# include "../misc/wrapper_mpi.hpp"
@@ -53,6 +55,44 @@ void declare_HMatrix(py::module &m, const std::string &className) {
5355
py_class.def("get_target_cluster", &HMatrix<CoefficientPrecision, CoordinatePrecision>::get_target_cluster, py::return_value_policy::reference_internal);
5456
py_class.def("get_source_cluster", &HMatrix<CoefficientPrecision, CoordinatePrecision>::get_source_cluster, py::return_value_policy::reference_internal);
5557

58+
py_class.def("lu_factorization", [](HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix) {
59+
htool::lu_factorization(hmatrix);
60+
});
61+
py_class.def("cholesky_factorization", [](HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix, char UPLO) {
62+
htool::cholesky_factorization(UPLO, hmatrix);
63+
});
64+
py_class.def("lu_solve", [](const Class &self, char trans, const py::array_t<CoefficientPrecision, py::array::f_style> &input) {
65+
std::vector<ssize_t> shape;
66+
if (input.ndim() == 1) {
67+
shape = {input.shape()[0]};
68+
} else if (input.ndim() == 2) {
69+
shape = {input.shape()[0], input.shape()[1]};
70+
} else {
71+
throw std::runtime_error("Wrong dimension for HMatrix-LU input"); // LCOV_EXCL_LINE
72+
}
73+
py::array_t<CoefficientPrecision, py::array::f_style> result(shape);
74+
std::copy_n(input.data(), input.size(), result.mutable_data());
75+
htool::MatrixView<CoefficientPrecision> output_view(result.shape()[0], input.ndim() == 1 ? 1 : result.shape()[1], result.mutable_data());
76+
htool::lu_solve(trans, self, output_view);
77+
return result;
78+
});
79+
80+
py_class.def("cholesky_solve", [](const Class &self, char UPLO, const py::array_t<CoefficientPrecision, py::array::f_style> &input) {
81+
std::vector<ssize_t> shape;
82+
if (input.ndim() == 1) {
83+
shape = {input.shape()[0]};
84+
} else if (input.ndim() == 2) {
85+
shape = {input.shape()[0], input.shape()[1]};
86+
} else {
87+
throw std::runtime_error("Wrong dimension for HMatrix-Cholesky input"); // LCOV_EXCL_LINE
88+
}
89+
py::array_t<CoefficientPrecision, py::array::f_style> result(shape);
90+
std::copy_n(input.data(), input.size(), result.mutable_data());
91+
htool::MatrixView<CoefficientPrecision> output_view(result.shape()[0], input.ndim() == 1 ? 1 : result.shape()[1], result.mutable_data());
92+
htool::cholesky_solve(UPLO, self, output_view);
93+
return result;
94+
});
95+
5696
m.def("recompression", &htool::recompression<CoefficientPrecision, CoordinatePrecision, std::function<void(LowRankMatrix<CoefficientPrecision> &)>>);
5797
m.def("recompression", [](HMatrix<CoefficientPrecision, CoordinatePrecision> &hmatrix) { recompression(hmatrix); });
5898
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) {
69109
py::array_t<CoefficientPrecision, py::array::f_style> result(self.get_target_cluster().get_size());
70110
std::fill_n(result.mutable_data(), self.get_target_cluster().get_size(), CoefficientPrecision(0));
71111

72-
htool::Matrix<CoefficientPrecision> dense_mat(self.get_target_cluster().get_size(), self.get_source_cluster().get_size());
73-
copy_to_dense_in_user_numbering(self, dense_mat.data());
74112
char trans = 'N';
75113
htool::add_hmatrix_vector_product(trans, CoefficientPrecision(1), self, input.data(), CoefficientPrecision(0), result.mutable_data());
76114

77115
return result;
78116
},
79117
"in"_a);
118+
119+
py_class.def(
120+
"__matmul__", [](const Class &self, const py::array_t<CoefficientPrecision, py::array::f_style> input) {
121+
if (input.ndim() != 2) {
122+
throw std::runtime_error("Wrong dimension for HMatrix-matrix product"); // LCOV_EXCL_LINE
123+
}
124+
if (input.shape()[0] != self.get_source_cluster().get_size()) {
125+
throw std::runtime_error("Wrong size for HMatrix-matrix product"); // LCOV_EXCL_LINE
126+
}
127+
py::array_t<CoefficientPrecision, py::array::f_style> result({input.shape()[0], input.shape()[1]});
128+
std::fill_n(result.mutable_data(), input.shape()[0] * input.shape()[1], CoefficientPrecision(0));
129+
130+
htool::MatrixView<const CoefficientPrecision> input_view(input.shape()[0], input.shape()[1], input.data());
131+
htool::MatrixView<CoefficientPrecision> output_view(input.shape()[0], input.shape()[1], result.mutable_data());
132+
char transa = 'N';
133+
char transb = 'N';
134+
htool::add_hmatrix_matrix_product(transa, transb, CoefficientPrecision(1), self, input_view, CoefficientPrecision(0), output_view);
135+
136+
return result;
137+
},
138+
"in"_a);
80139
}
81140

82141
#endif

‎src/htool/hmatrix/hmatrix_tree_builder.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,5 +40,6 @@ void declare_hmatrix_builder(py::module &m, const std::string &className) {
4040
py_class.def("set_minimal_target_depth", &Class::set_minimal_target_depth);
4141
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); });
4242
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); });
43+
py_class.def("set_block_tree_consistency", [](Class &self, bool consistency) { self.set_block_tree_consistency(consistency); });
4344
}
4445
#endif

‎src/htool/main.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -84,6 +84,7 @@ PYBIND11_MODULE(Htool, m) {
8484

8585
declare_matplotlib_cluster<double>(m);
8686
declare_matplotlib_hmatrix<double, double>(m);
87+
declare_matplotlib_hmatrix<std::complex<double>, double>(m);
8788

8889
declare_virtual_partitioning<std::complex<double>>(m, "Complex");
8990
declare_LowRankMatrix<std::complex<double>>(m, "ComplexLowRankMatrix");

‎tests/conftest.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,11 +2,11 @@
22
import pathlib
33
import struct
44

5-
import Htool
65
import mpi4py
76
import numpy as np
87
import pytest
98

9+
import Htool
1010
from example.advanced.define_custom_dense_blocks_generator import (
1111
CustomDenseBlocksGenerator,
1212
)

0 commit comments

Comments
 (0)