Skip to content

Commit beedf2f

Browse files
PierreMarchand20PierreMarchand20
committed
improve testing and fix virtual low rank generator with threading
1 parent bae8293 commit beedf2f

2 files changed

Lines changed: 49 additions & 25 deletions

File tree

‎src/htool/hmatrix/interfaces/virtual_low_rank_generator.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,11 +23,11 @@ class VirtualLowRankGeneratorPython : public VirtualLowRankGenerator<Coefficient
2323
VirtualLowRankGeneratorPython(bool allow_copy = true) : VirtualLowRankGenerator<CoefficientPrecision>(), m_allow_copy(allow_copy) {}
2424

2525
bool copy_low_rank_approximation(int M, int N, const int *const rows, const int *const cols, LowRankMatrix<CoefficientPrecision> &lrmat) const override {
26+
py::gil_scoped_acquire acquire;
2627
auto &U = lrmat.get_U();
2728
auto &V = lrmat.get_V();
2829
py::array_t<int> py_rows(std::array<long int, 1>{M}, rows, py::capsule(rows));
2930
py::array_t<int> py_cols(std::array<long int, 1>{N}, cols, py::capsule(cols));
30-
3131
bool success = build_low_rank_approximation(py_rows, py_cols, lrmat.get_epsilon());
3232
if (success) {
3333
if (m_allow_copy) {

‎tests/test_hmatrix.py‎

Lines changed: 48 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -5,23 +5,27 @@
55
import pytest
66

77
import Htool
8-
from example.advanced.define_custom_low_rank_generator import CustomSVD
8+
from example.advanced.define_custom_low_rank_generator import (
9+
ComplexCustomSVD,
10+
CustomSVD,
11+
)
912
from example.create_geometry import create_random_geometries
10-
from example.define_generators import CustomGenerator
13+
from example.define_generators import ComplexCustomGenerator, CustomGenerator
1114

1215

1316
@pytest.mark.parametrize(
14-
"loglevel,symmetry",
17+
"loglevel,symmetry,is_complex",
1518
[
16-
(logging.INFO, "N"),
17-
(logging.DEBUG, "N"),
18-
(logging.WARNING, "N"),
19-
(logging.ERROR, "N"),
20-
(logging.CRITICAL, "N"),
21-
(logging.INFO, "S"),
19+
(logging.INFO, "N", False),
20+
(logging.DEBUG, "N", False),
21+
(logging.WARNING, "N", False),
22+
(logging.ERROR, "N", False),
23+
(logging.CRITICAL, "N", False),
24+
(logging.INFO, "S", False),
25+
(logging.INFO, "S", True),
2226
],
2327
)
24-
def test_hmatrix(loglevel, symmetry):
28+
def test_hmatrix(loglevel, symmetry, is_complex):
2529
logging.basicConfig(level=loglevel)
2630

2731
# Random geometry
@@ -52,17 +56,30 @@ def test_hmatrix(loglevel, symmetry):
5256
source_cluster = target_cluster
5357

5458
# Build generator
55-
if symmetry == "N":
56-
generator = CustomGenerator(target_points, source_points)
59+
if is_complex is False:
60+
if symmetry == "N":
61+
generator = CustomGenerator(target_points, source_points)
62+
else:
63+
generator = CustomGenerator(target_points, target_points)
5764
else:
58-
generator = CustomGenerator(target_points, target_points)
59-
65+
if symmetry == "N":
66+
generator = ComplexCustomGenerator(target_points, source_points)
67+
else:
68+
generator = ComplexCustomGenerator(target_points, target_points)
6069
# Custom low rank generator
61-
low_rank_generator = CustomSVD(generator, False)
70+
if is_complex is False:
71+
low_rank_generator = CustomSVD(generator, False)
72+
else:
73+
low_rank_generator = ComplexCustomSVD(generator, False)
6274

6375
# Build HMatrix
64-
hmatrix_builder = Htool.HMatrixTreeBuilder(epsilon, eta, "N", "N")
76+
if is_complex is False:
77+
hmatrix_builder = Htool.HMatrixTreeBuilder(epsilon, eta, "N", "N")
78+
else:
79+
hmatrix_builder = Htool.ComplexHMatrixTreeBuilder(epsilon, eta, "N", "N")
6580
hmatrix_builder.set_low_rank_generator(low_rank_generator)
81+
if symmetry == "S":
82+
hmatrix_builder.set_block_tree_consistency(True)
6683
hmatrix = hmatrix_builder.build(generator, target_cluster, source_cluster)
6784
assert hmatrix.shape == (nb_rows, nb_cols)
6885

@@ -74,39 +91,46 @@ def test_hmatrix(loglevel, symmetry):
7491
dense_in_user_numbering = hmatrix.to_dense_in_user_numbering()
7592

7693
# HMatrix vector product
94+
dtype = np.float64 if is_complex is False else np.complex128
7795
np.random.seed(0)
78-
x = np.random.rand(nb_cols)
96+
if is_complex is False:
97+
x = np.random.rand(nb_cols)
98+
else:
99+
x = np.random.rand(nb_cols) + 1j * np.random.rand(nb_cols)
79100
y = hmatrix * x
80101
y_exact = generator.mat_vec(x)
81102
y_dense = dense_in_user_numbering.dot(x)
82103
y_copy = copy_hmatrix * x
83104
assert np.linalg.norm(y - y_exact) / np.linalg.norm(y_exact) < epsilon
84105
assert np.linalg.norm(y - y_dense) / np.linalg.norm(y_dense) < 1e-10
85-
assert np.linalg.norm(y - y_copy) < 1e-10
106+
assert np.linalg.norm(y - y_copy) / np.linalg.norm(y) < 1e-10
86107

87108
# HMatrix matrix product
88109
np.random.seed(0)
89-
x = np.random.rand(nb_cols, 2)
110+
if is_complex is False:
111+
x = np.random.rand(nb_cols, 2)
112+
else:
113+
x = np.random.rand(nb_cols, 2) + 1j * np.random.rand(nb_cols, 2)
90114
y = hmatrix @ x
91115
y_exact = generator.mat_mat(x)
92116
y_dense = dense_in_user_numbering @ x
93117
y_copy = copy_hmatrix @ x
94118
assert np.linalg.norm(y - y_exact) / np.linalg.norm(y_exact) < epsilon
95119
assert np.linalg.norm(y - y_dense) / np.linalg.norm(y_dense) < 1e-10
96-
assert np.linalg.norm(y - y_copy) < 1e-10
120+
assert np.linalg.norm(y - y_copy) / np.linalg.norm(y) < 1e-10
97121

98122
if symmetry != "N":
99123
# HLU factorization
100124
copy_hmatrix.lu_factorization()
101125

102126
# HLU solve vec
103-
x_ref = np.ones(nb_cols)
127+
x_ref = np.ones(nb_cols, dtype=dtype)
104128
y = hmatrix * x_ref
105129
x = copy_hmatrix.lu_solve("N", y)
106130
assert np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref) < epsilon
107131

108132
# HLU solve mat
109-
x_ref = np.ones((nb_cols, 2))
133+
x_ref = np.ones((nb_cols, 2), dtype=dtype)
110134
y = hmatrix @ x_ref
111135
x = copy_hmatrix.lu_solve("N", y)
112136
assert np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref) < epsilon
@@ -116,13 +140,13 @@ def test_hmatrix(loglevel, symmetry):
116140
copy_hmatrix.cholesky_factorization("L")
117141

118142
# Cholesky solve vec
119-
x_ref = np.ones(nb_cols)
143+
x_ref = np.ones(nb_cols, dtype=dtype)
120144
y = hmatrix * x_ref
121145
x = copy_hmatrix.cholesky_solve("L", y)
122146
assert np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref) < epsilon
123147

124148
# Cholesky solve mat
125-
x_ref = np.ones((nb_cols, 2))
149+
x_ref = np.ones((nb_cols, 2), dtype=dtype)
126150
y = hmatrix @ x_ref
127151
x = copy_hmatrix.cholesky_solve("L", y)
128152
assert np.linalg.norm(x - x_ref) / np.linalg.norm(x_ref) < epsilon

0 commit comments

Comments
 (0)