diff --git a/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp b/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp
index 9b7fb8a5..e5d0b2bc 100644
--- a/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp
+++ b/tests/functional_tests/hmatrix/test_hmatrix_factorization.hpp
@@ -44,7 +44,42 @@ bool test_hmatrix_lu(char trans, int n1, int n2, htool::underlying_type<T> epsil
     return is_error;
 }
 
-template <typename T, typename GeneratorTestType>
+template <typename T, typename GeneratorTestType, std::enable_if_t<!is_complex_t<T>::value, bool> = true>
+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;
+    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>> HA = hmatrix_tree_builder_A.build(*test_case.operator_A);
+
+    // Matrix
+    int ni_A = test_case.root_cluster_A_input->get_size();
+    int no_A = test_case.root_cluster_A_output->get_size();
+    int ni_X = test_case.root_cluster_X_input->get_size();
+    int no_X = test_case.root_cluster_X_output->get_size();
+    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_symmetric_matrix_matrix_product('L', UPLO, T(1.), A_dense, X_dense, T(0.), B_dense);
+
+    // Cholesky factorization
+    matrix_test = B_dense;
+    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;
+    cout << "> is_error: " << is_error << "\n";
+
+    return is_error;
+}
+
+template <typename T, typename GeneratorTestType, std::enable_if_t<is_complex_t<T>::value, bool> = true>
 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;
@@ -65,11 +100,7 @@ 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);
-    if (is_complex<T>()) {
-        add_hermitian_matrix_matrix_product('L', UPLO, T(1.), A_dense, X_dense, T(0.), B_dense);
-    } else {
-        add_symmetric_matrix_matrix_product('L', UPLO, T(1.), A_dense, X_dense, T(0.), B_dense);
-    }
+    add_hermitian_matrix_matrix_product('L', UPLO, T(1.), A_dense, X_dense, T(0.), B_dense);
 
     // Cholesky factorization
     matrix_test = B_dense;