linear_solve.h

Go to the documentation of this file.

#ifndef FUSED_ALGORITHMS_LINEAR_SOLVE_H
#define FUSED_ALGORITHMS_LINEAR_SOLVE_H
 
#include "config.h"
#include "utilities/expected.h"
#include "matrix_traits.h"
#include "algorithms/decomposition/lu.h"
#include "algorithms/solvers/triangular_solve.h"
#include "algorithms/solvers/cholesky_solve.h"
 
namespace matrix_algorithms
{
 
    using matrix_traits::MatrixStatus;
 
    template <typename T, my_size_t N>
    Expected<FusedVector<T, N>, MatrixStatus> lu_solve(
        const FusedMatrix<T, N, N> &A,
        const FusedVector<T, N> &b)
    {
        static_assert(is_floating_point_v<T>,
                      "lu_solve requires a floating-point scalar type");
 
        // 1. Decompose P·A = L·U
        auto lu_result = lu(A);
 
        if (!lu_result.has_value())
        {
            return Unexpected{lu_result.error()};
        }
 
        auto &decomp = lu_result.value();
 
        // 2. Permute b
        FusedVector<T, N> b_perm(T(0));
 
        for (my_size_t i = 0; i < N; ++i)
        {
            b_perm(i) = b(decomp.perm(i));
        }
 
        // 3. Extract L and U
        auto L = decomp.L();
        auto U = decomp.U();
 
        // 4. Solve L·y = b_perm (forward substitution, unit diagonal)
        auto fwd_result = forward_substitute<true>(L, b_perm);
 
        if (!fwd_result.has_value())
        {
            return Unexpected{fwd_result.error()};
        }
 
        auto &y = fwd_result.value();
 
        // 5. Solve U·x = y (back substitution)
        return back_substitute(U, y);
    }
 
    template <typename T, my_size_t N>
    Expected<FusedVector<T, N>, MatrixStatus> solve(
        const FusedMatrix<T, N, N> &A,
        const FusedVector<T, N> &b)
    {
        static_assert(is_floating_point_v<T>,
                      "solve requires a floating-point scalar type");
 
        // Try Cholesky path first (O(N³/3) — cheaper if SPD)
        auto chol_result = cholesky_solve(A, b);
 
        if (chol_result.has_value())
        {
            return chol_result;
        }
 
        // Fall back to LU path (O(2N³/3) — works for any non-singular A)
        return lu_solve(A, b);
    }
 
} // namespace matrix_algorithms
 
#endif // FUSED_ALGORITHMS_LINEAR_SOLVE_H