Skip to content

Commit

Permalink
fix cholesky factorization
Browse files Browse the repository at this point in the history
  • Loading branch information
@PierreMarchand20
PierreMarchand20 committed Jul 18, 2024
1 parent 794be2d commit ee8dbd3
Show file tree
Hide file tree
Showing 3 changed files with 12 additions and 24 deletions.
8 changes: 4 additions & 4 deletions include/htool/hmatrix/linalg/factorization.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,10 +104,10 @@ void cholesky_factorization(char UPLO, HMatrix<CoefficientPrecision, CoordinateP
if (other_cluster_child->get_offset() > cluster_child->get_offset()) {
if (UPLO == 'L') {
HMatrix<CoefficientPrecision, CoordinatePrecision> *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<CoefficientPrecision>() ? 'C' : 'T', 'N', CoefficientPrecision(1), *pivot, *L);
} else {
HMatrix<CoefficientPrecision, CoordinatePrecision> *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<CoefficientPrecision>() ? 'C' : 'T', 'N', CoefficientPrecision(1), *pivot, *U);
}
}
}
Expand All @@ -122,15 +122,15 @@ void cholesky_factorization(char UPLO, HMatrix<CoefficientPrecision, CoordinateP
// if (*output_cluster_child == *input_cluster_child) {
// symmetric_rank_k_update(UPLO, 'N', CoefficientPrecision(-1), *L, CoefficientPrecision(1), *A_child);
// } else {
add_hmatrix_hmatrix_product('N', 'C', CoefficientPrecision(-1), *L, *L, CoefficientPrecision(1), *A_child);
add_hmatrix_hmatrix_product('N', is_complex<CoefficientPrecision>() ? '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<CoefficientPrecision, CoordinatePrecision> *A_child = hmatrix.get_sub_hmatrix(*output_cluster_child, *input_cluster_child);
const HMatrix<CoefficientPrecision, CoordinatePrecision> *U = hmatrix.get_sub_hmatrix(*cluster_child, *input_cluster_child);
// if (*output_cluster_child == *input_cluster_child) {
// symmetric_rank_k_update(UPLO, 'C', CoefficientPrecision(-1), *U, CoefficientPrecision(1), *A_child);
// } else {
add_hmatrix_hmatrix_product('C', 'N', CoefficientPrecision(-1), *U, *U, CoefficientPrecision(1), *A_child);
add_hmatrix_hmatrix_product(is_complex<CoefficientPrecision>() ? 'C' : 'T', 'N', CoefficientPrecision(-1), *U, *U, CoefficientPrecision(1), *A_child);
// }
}
}
Expand Down
6 changes: 3 additions & 3 deletions tests/functional_tests/hmatrix/hmatrix_factorization/test_hmatrix_factorization_double.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,9 @@ int main(int, char *[]) {
for (auto trans : {'N', 'T'}) {
is_error = is_error || test_hmatrix_lu<double, GeneratorTestDoubleSymmetric>(trans, n1, n2, epsilon, margin);
}
// for (auto UPLO : {'L', 'U'}) {
// is_error = is_error || test_hmatrix_cholesky<double, GeneratorTestDoubleSymmetric>(UPLO, n1, n2, epsilon, margin);
// }
for (auto UPLO : {'L', 'U'}) {
is_error = is_error || test_hmatrix_cholesky<double, GeneratorTestDoubleSymmetric>(UPLO, n1, n2, epsilon, margin);
}
}

if (is_error) {
Expand Down
22 changes: 5 additions & 17 deletions tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,17 +47,15 @@ bool test_hmatrix_lu(char trans, int n1, int n2, htool::underlying_type<T> epsil
template <typename T, typename GeneratorTestType>
bool test_hmatrix_cholesky(char UPLO, int n1, int n2, htool::underlying_type<T> epsilon, htool::underlying_type<T> margin) {
bool is_error = false;
double eta = -100;
double eta = 100;
htool::underlying_type<T> error;

// Setup test case
htool::TestCaseSolve<T, GeneratorTestType> test_case('L', 'N', n1, n2, 1, -1);

// HMatrix
HMatrixTreeBuilder<T, htool::underlying_type<T>> hmatrix_tree_builder_A(*test_case.root_cluster_A_output, *test_case.root_cluster_A_input, epsilon, eta, is_complex<T>() ? 'H' : 'S', UPLO, -1, -1, -1);
HMatrix<T, htool::underlying_type<T>> A = hmatrix_tree_builder_A.build(*test_case.operator_A);

save_leaves_with_rank(A, "A");
HMatrix<T, htool::underlying_type<T>> HA = hmatrix_tree_builder_A.build(*test_case.operator_A);

// Matrix
int ni_A = test_case.root_cluster_A_input->get_size();
Expand All @@ -67,22 +65,12 @@ bool test_hmatrix_cholesky(char UPLO, int n1, int n2, htool::underlying_type<T>
Matrix<T> 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<T> 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;
Expand Down

0 comments on commit ee8dbd3

Please sign in to comment.