From ee8dbd39991b4638502c3909a4d91f467e2c82e6 Mon Sep 17 00:00:00 2001 From: Pierre Marchand Date: Thu, 18 Jul 2024 19:25:04 +0200 Subject: [PATCH] fix cholesky factorization --- .../htool/hmatrix/linalg/factorization.hpp | 8 +++---- .../test_hmatrix_factorization_double.cpp | 6 ++--- .../hmatrix/test_hmatrix_factorization.hpp | 22 +++++-------------- 3 files changed, 12 insertions(+), 24 deletions(-) diff --git a/include/htool/hmatrix/linalg/factorization.hpp b/include/htool/hmatrix/linalg/factorization.hpp index 5feb76cd..807ef145 100644 --- a/include/htool/hmatrix/linalg/factorization.hpp +++ b/include/htool/hmatrix/linalg/factorization.hpp @@ -104,10 +104,10 @@ void cholesky_factorization(char UPLO, HMatrixget_offset() > cluster_child->get_offset()) { if (UPLO == 'L') { HMatrix *L = hmatrix.get_sub_hmatrix(*other_cluster_child, *cluster_child); - triangular_hmatrix_hmatrix_solve('R', UPLO, 'C', 'N', CoefficientPrecision(1), *pivot, *L); + triangular_hmatrix_hmatrix_solve('R', UPLO, is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(1), *pivot, *L); } else { HMatrix *U = hmatrix.get_sub_hmatrix(*cluster_child, *other_cluster_child); - triangular_hmatrix_hmatrix_solve('L', UPLO, 'C', 'N', CoefficientPrecision(1), *pivot, *U); + triangular_hmatrix_hmatrix_solve('L', UPLO, is_complex() ? 'C' : 'T', 'N', CoefficientPrecision(1), *pivot, *U); } } } @@ -122,7 +122,7 @@ void cholesky_factorization(char UPLO, HMatrix() ? 'C' : 'T', CoefficientPrecision(-1), *L, *L, CoefficientPrecision(1), *A_child); // } } else if (UPLO == 'U' && output_cluster_child->get_offset() > cluster_child->get_offset() && input_cluster_child->get_offset() > cluster_child->get_offset() && input_cluster_child->get_offset() >= output_cluster_child->get_offset()) { HMatrix *A_child = hmatrix.get_sub_hmatrix(*output_cluster_child, *input_cluster_child); @@ -130,7 +130,7 @@ void cholesky_factorization(char UPLO, HMatrix() ? 'C' : 'T', 'N', CoefficientPrecision(-1), *U, *U, CoefficientPrecision(1), *A_child); // } } } diff --git a/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_double.cpp b/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_double.cpp index f68037f4..76b6505a 100644 --- a/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_double.cpp +++ b/tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_double.cpp @@ -18,9 +18,9 @@ int main(int, char *[]) { for (auto trans : {'N', 'T'}) { is_error = is_error || test_hmatrix_lu(trans, n1, n2, epsilon, margin); } - // for (auto UPLO : {'L', 'U'}) { - // is_error = is_error || test_hmatrix_cholesky(UPLO, n1, n2, epsilon, margin); - // } + for (auto UPLO : {'L', 'U'}) { + is_error = is_error || test_hmatrix_cholesky(UPLO, n1, n2, epsilon, margin); + } } if (is_error) { diff --git a/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp b/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp index 771b9089..9380dd55 100644 --- a/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp +++ b/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp @@ -47,7 +47,7 @@ bool test_hmatrix_lu(char trans, int n1, int n2, htool::underlying_type epsil template bool test_hmatrix_cholesky(char UPLO, int n1, int n2, htool::underlying_type epsilon, htool::underlying_type margin) { bool is_error = false; - double eta = -100; + double eta = 100; htool::underlying_type error; // Setup test case @@ -55,9 +55,7 @@ bool test_hmatrix_cholesky(char UPLO, int n1, int n2, htool::underlying_type // HMatrix HMatrixTreeBuilder> hmatrix_tree_builder_A(*test_case.root_cluster_A_output, *test_case.root_cluster_A_input, epsilon, eta, is_complex() ? 'H' : 'S', UPLO, -1, -1, -1); - HMatrix> A = hmatrix_tree_builder_A.build(*test_case.operator_A); - - save_leaves_with_rank(A, "A"); + HMatrix> HA = hmatrix_tree_builder_A.build(*test_case.operator_A); // Matrix int ni_A = test_case.root_cluster_A_input->get_size(); @@ -67,22 +65,12 @@ bool test_hmatrix_cholesky(char UPLO, int n1, int n2, htool::underlying_type Matrix A_dense(no_A, ni_A), X_dense(no_X, ni_X), B_dense(X_dense), densified_hmatrix_test(B_dense), matrix_test; test_case.operator_A->copy_submatrix(no_A, ni_A, test_case.root_cluster_A_output->get_offset(), test_case.root_cluster_A_input->get_offset(), A_dense.data()); generate_random_matrix(X_dense); - // add_matrix_matrix_product('N', 'N', T(1.), A_dense, X_dense, T(0.), B_dense); add_symmetric_matrix_matrix_product('L', UPLO, T(1.), A_dense, X_dense, T(0.), B_dense); - // LU factorization + // Cholesky factorization matrix_test = B_dense; - cholesky_factorization(UPLO, A); - save_leaves_with_rank(A, "A_chol"); - Matrix A_chol_dense(A.nb_rows(), A.nb_cols()); - copy_to_dense(A, A_chol_dense.data()); - std::ofstream A_chol_file("A_chol_dense"); - A_chol_dense.print(A_chol_file, ","); - cholesky_solve(UPLO, A, matrix_test); - std::ofstream matrix_test_file("matrix_test"); - matrix_test.print(matrix_test_file, ","); - std::ofstream X_dense_file("X_dense"); - X_dense.print(X_dense_file, ","); + cholesky_factorization(UPLO, HA); + cholesky_solve(UPLO, HA, matrix_test); error = normFrob(X_dense - matrix_test) / normFrob(X_dense); is_error = is_error || !(error < epsilon * margin); cout << "> Errors on hmatrix cholesky solve: " << error << endl;