Add LDLt factorization for HMatrix and use it in the local DDM solvers

CHANGELOG.md

@@ -30,6 +30,10 @@ All notable changes to this project will be documented in this file.
 
 - JOSS badge and update `CITATION.cff`.
 - Opt-in `tidy` CMake target running clang-tidy (`bugprone-*`/`performance-*`) over `examples/` and `tests/`, not part of the default build. Configured via `.clang-tidy`: findings are errors except `branch-clone`, `exception-escape` and `implicit-widening-of-multiplication-result`; the macro-generated `wrapper_blas.hpp`/`wrapper_lapack.hpp` are excluded from findings (needs clang-tidy >= 19). Run in CI on the Ubuntu clang/OpenMPI job.
+- LDLt factorization (Bunch-Kaufman) for symmetric, complex-symmetric and Hermitian matrices, selected by a `symmetry` of `'S'` or `'H'`:
+  - dense: `ldlt_factorization`/`ldlt_solve` (`sytrf`/`sytrs`, `hetrf`/`hetrs`), `triangular_ldlt_matrix_matrix_solve` (unit L/U with `transa` in N/T/C) and `apply_ldlt_diagonal` (D^-1), plus `Blas::geru`/`Blas::gerc` (BLAS has no plain complex `ger`).
+  - HMatrix: `ldlt_factorization`/`sequential_ldlt_factorization` and `ldlt_solve`, also on a full-storage HMatrix (symmetry `'N'`) like `cholesky_factorization`, with `internal_triangular_ldlt_hmatrix_{hmatrix,matrix,lrmat}_solve` and `internal_apply_ldlt_diagonal`.
+- `is_positive_definite` option (default `false`) in `LocalHMatrixSolver`, `LocalHMatrixPlusOverlapSolver`, `make_DDM_solver_w_custom_local_solver` and `DDMSolverBuilder`: a symmetric/Hermitian local HMatrix is factorized with Cholesky only when it is set, and with LDLt otherwise.
 
 ### Changed
 
@@ -39,6 +43,8 @@ All notable changes to this project will be documented in this file.
 - Replace `std::endl`/`endl` with `'\n'` across `include/`, `tests/` and `examples/`.
 - Fix unnecessary-value-param findings from clang-tidy: move sink parameters, take pure reads by const reference.
 - Add explicit `default` cases to `DDM`'s `schwarz_method`/`krylov_method` switches (behavior unchanged).
+- Symmetric/Hermitian local HMatrix solvers (`DDMSolverBuilder`) now use LDLt by default instead of Cholesky, which silently gave a wrong factorization for indefinite matrices; pass `is_positive_definite=true` to keep Cholesky.
+- Checks of `lu_factorization`/`cholesky_factorization` on HMatrix moved to the public functions taking an execution policy, so they run once for the sequential and task-based implementations.
 
 ### Removed
 
@@ -56,6 +62,9 @@ All notable changes to this project will be documented in this file.
 - Fix `Matrix`'s int-overflow risk in allocation sizes (rows*cols computed in `int` then used as a `new[]`/read/write size), and a zero-length-array leak in `resize()`.
 - Fix real overflow/mistyping bugs surfaced by clang-tidy's widening/narrowing-conversion checks in `matrix/utils/output.hpp` and `hmatrix_output.hpp`.
 - Compute the partition-level power of the number of children as an exact integer instead of via `std::pow` in cluster tree building.
+- Fix HMatrix `cholesky_factorization` symmetry check, which compared `UPLO` instead of the symmetry and could never trigger: it now rejects a complex symmetric (non-Hermitian) HMatrix.
+- Fix complex-symmetric (non-Hermitian) local HMatrix solvers, which were factorized as Hermitian with Cholesky.
+- Report a failed dense `cholesky_factorization` (`potrf`'s `info`), for a matrix that is not positive definite.
 
 ## [1.0.2] - 2026-02-14
 

include/htool/hmatrix/linalg/apply_ldlt_diagonal.hpp

@@ -0,0 +1,73 @@
+#ifndef HTOOL_HMATRIX_LINALG_APPLY_LDLT_DIAGONAL_HPP
+#define HTOOL_HMATRIX_LINALG_APPLY_LDLT_DIAGONAL_HPP
+
+#include "../../matrix/linalg/factorization.hpp"
+#include "../../matrix/matrix.hpp"
+#include "../../misc/logger.hpp"
+#include "../hmatrix.hpp"
+#include <algorithm>
+#include <memory>
+#include <vector>
+
+namespace htool {
+// B := D^-1 B (side='L') or B D^-1 (side='R'), D being the block-diagonal factor of an LDLt
+// factorization held at A's dense diagonal leaves (see apply_ldlt_diagonal). B's rows (side='L')
+// or columns (side='R') are numbered from offset.
+template <typename CoefficientPrecision, typename CoordinatePrecision>
+void internal_apply_ldlt_diagonal(char symmetry, char side, char UPLO, const HMatrix<CoefficientPrecision, CoordinatePrecision> &A, Matrix<CoefficientPrecision> &B, int offset) {
+    if (A.is_hierarchical()) {
+        for (auto &child : A.get_children()) {
+            if (child->get_target_cluster() == child->get_source_cluster()) {
+                internal_apply_ldlt_diagonal(symmetry, side, UPLO, *child, B, offset);
+            }
+        }
+    } else {
+        int local_offset = A.get_target_cluster().get_offset() - offset;
+        int size         = A.get_target_cluster().get_size();
+        if (side == 'R') {
+            Matrix<CoefficientPrecision> B_columns;
+            B_columns.assign(B.nb_rows(), size, B.data() + B.nb_rows() * local_offset, false);
+            apply_ldlt_diagonal(symmetry, 'R', UPLO, *A.get_dense_data(), B_columns);
+        } else {
+            Matrix<CoefficientPrecision> B_rows(size, B.nb_cols());
+            for (int j = 0; j < B.nb_cols(); j++) {
+                std::copy_n(B.data() + local_offset + j * B.nb_rows(), size, B_rows.data() + j * size);
+            }
+            apply_ldlt_diagonal(symmetry, 'L', UPLO, *A.get_dense_data(), B_rows);
+            for (int j = 0; j < B.nb_cols(); j++) {
+                std::copy_n(B_rows.data() + j * size, size, B.data() + local_offset + j * B.nb_rows());
+            }
+        }
+    }
+}
+
+// Same with an HMatrix B, following B's own block structure so that its blocks keep their compression.
+template <typename CoefficientPrecision, typename CoordinatePrecision>
+void internal_apply_ldlt_diagonal(char symmetry, char side, char UPLO, const HMatrix<CoefficientPrecision, CoordinatePrecision> &A, HMatrix<CoefficientPrecision, CoordinatePrecision> &B) {
+    int offset = side == 'L' ? B.get_target_cluster().get_offset() : B.get_source_cluster().get_offset();
+    if (B.is_dense()) {
+        internal_apply_ldlt_diagonal(symmetry, side, UPLO, A, *B.get_dense_data(), offset);
+    } else if (B.is_low_rank()) {
+        internal_apply_ldlt_diagonal(symmetry, side, UPLO, A, side == 'L' ? B.get_low_rank_data()->get_U() : B.get_low_rank_data()->get_V(), offset);
+    } else {
+        std::vector<const HMatrix<CoefficientPrecision, CoordinatePrecision> *> A_children;
+        for (auto &child : B.get_children()) {
+            const Cluster<CoordinatePrecision> &cluster = side == 'L' ? child->get_target_cluster() : child->get_source_cluster();
+            A_children.push_back(A.get_sub_hmatrix(cluster, cluster));
+        }
+        if (std::find(A_children.begin(), A_children.end(), nullptr) == A_children.end()) {
+            for (std::size_t c = 0; c < A_children.size(); c++) {
+                internal_apply_ldlt_diagonal(symmetry, side, UPLO, *A_children[c], *B.get_children()[c]);
+            }
+            return;
+        }
+        // B is refined below A's dense leaves, where a 2x2 pivot block of D may straddle B's children.
+        htool::Logger::get_instance().log(LogLevel::WARNING, "In internal_apply_ldlt_diagonal, B is refined below the dense diagonal leaves holding D and is densified, losing its compression.");
+        auto B_dense = std::make_unique<Matrix<CoefficientPrecision>>(B.nb_rows(), B.nb_cols());
+        copy_to_dense(B, B_dense->data());
+        internal_apply_ldlt_diagonal(symmetry, side, UPLO, A, *B_dense, offset);
+        B.set_dense_data(std::move(B_dense));
+    }
+}
+} // namespace htool
+#endif