Solve Ax = b for general or symmetric positive-definite A. More...
#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"
Include dependency graph for linear_solve.h:
Go to the source code of this file.
Namespaces | |
| namespace | matrix_algorithms |
Functions | |
| template<typename T , my_size_t N> | |
| Expected< FusedVector< T, N >, MatrixStatus > | matrix_algorithms::lu_solve (const FusedMatrix< T, N, N > &A, const FusedVector< T, N > &b) |
| Solve Ax = b via LU decomposition with partial pivoting. | |
| template<typename T , my_size_t N> | |
| Expected< FusedVector< T, N >, MatrixStatus > | matrix_algorithms::solve (const FusedMatrix< T, N, N > &A, const FusedVector< T, N > &b) |
| Solve Ax = b with automatic algorithm selection. | |
Detailed Description
Solve Ax = b for general or symmetric positive-definite A.
Provides:
- lu_solve(A, b): always uses LU path — works for any non-singular A
- solve(A, b): auto-dispatches to Cholesky (if SPD) or LU
ALGORITHM — lu_solve
- Decompose P·A = L·U via lu(A)
- Permute b: b_perm(i) = b(perm(i))
- Solve L·y = b_perm via forward substitution (UnitDiag=true)
- Solve U·x = y via back substitution
Complexity: O(2N³/3) for LU + O(N²) for substitutions.
ALGORITHM — solve (dispatcher)
- Attempt cholesky_solve(A, b) — O(N³/3) if A is SPD
- If it fails (not symmetric or not positive definite), fall back to lu_solve(A, b) — O(2N³/3) for any non-singular A
The Cholesky attempt is not wasted work — if A is SPD, Cholesky is ~2× cheaper than LU. If A is not SPD, the failure is detected quickly (symmetry check is O(N²), positive-definiteness fails at the first non-positive diagonal during factorization).
FAILURE MODES
- MatrixStatus::Singular — from lu() if A is singular
- MatrixStatus::DimensionMismatch — if matrix/vector sizes don't match
