cholesky_solve.h

Go to the documentation of this file.

#ifndef FUSED_ALGORITHMS_CHOLESKY_SOLVE_H
#define FUSED_ALGORITHMS_CHOLESKY_SOLVE_H
 
#include "config.h"
#include "utilities/expected.h"
#include "matrix_traits.h"
#include "simple_type_traits.h"
#include "algorithms/decomposition/cholesky.h"
#include "algorithms/solvers/triangular_solve.h"
 
namespace matrix_algorithms
{
 
    using matrix_traits::MatrixStatus;
 
    template <typename T, my_size_t N>
    Expected<FusedVector<T, N>, MatrixStatus> cholesky_solve(
        const FusedMatrix<T, N, N> &A,
        const FusedVector<T, N> &b)
    {
        static_assert(is_floating_point_v<T>,
                      "cholesky_solve requires a floating-point scalar type");
 
        // 1. Decompose A = LLᵀ
        auto chol_result = cholesky(A);
 
        if (!chol_result.has_value())
        {
            return Unexpected{chol_result.error()};
        }
 
        auto &L = chol_result.value();
 
        // 2. Solve Ly = b (forward substitution)
        //    L comes from cholesky — guaranteed non-singular, but forward_substitute
        //    checks anyway; no overhead for the unrolled paths.
        auto fwd_result = forward_substitute(L, b);
 
        if (!fwd_result.has_value())
        {
            return Unexpected{fwd_result.error()};
        }
 
        auto &y = fwd_result.value();
 
        // 3. Solve Lᵀx = y (back substitution on L transposed)
        //    Access L(k,i) instead of U(i,k) to avoid creating a transpose.
        //    L diagonal is guaranteed positive from cholesky — no singular check needed.
        FusedVector<T, N> x(T(0));
 
        for (my_size_t i = N; i-- > 0;)
        {
            T sum = y(i);
 
            for (my_size_t k = i + 1; k < N; ++k)
            {
                sum -= L(k, i) * x(k); // L(k,i) == Lᵀ(i,k)
            }
 
            x(i) = sum / L(i, i);
        }
 
        return move(x);
    }
 
} // namespace matrix_algorithms
 
#endif // FUSED_ALGORITHMS_CHOLESKY_SOLVE_H