lu.h File Reference

LU decomposition with partial pivoting for square matrices. More...

#include "config.h"
#include "utilities/expected.h"
#include "matrix_traits.h"
#include "fused/fused_matrix.h"
#include "fused/fused_vector.h"
#include "math/math_utils.h"

Include dependency graph for lu.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	matrix_algorithms::LUResult< T, N >
	Result of LU decomposition with partial pivoting. More...

Namespaces
namespace	matrix_algorithms

Functions
template<typename T , my_size_t N>
Expected< LUResult< T, N >, MatrixStatus >	matrix_algorithms::lu (const FusedMatrix< T, N, N > &A, T tol=T(PRECISION_TOLERANCE))
	Compute the LU decomposition of a square matrix with partial pivoting.
template<typename T , my_size_t N>
LUResult< T, N >	matrix_algorithms::lu_or_die (const FusedMatrix< T, N, N > &A)
	LU decomposition — abort on failure.

Detailed Description

LU decomposition with partial pivoting for square matrices.

Decomposes A into P·A = L·U where:

• L is lower-triangular with unit diagonal (implicit, stored below diagonal)
• U is upper-triangular (stored on and above diagonal)
• P is a row permutation represented as an index vector

The compact representation stores L and U in a single N×N matrix (LAPACK-style). Accessor methods L() and U() extract separate matrices when needed.

---

ALGORITHM

For each column j = 0 … N−1:
  1. Find pivot: row p ≥ j with max |A(p,j)|
  2. Swap rows j and p in LU and record in perm; flip sign
  3. Check for singularity: |U(j,j)| ≤ PRECISION_TOLERANCE
  4. Eliminate: for each row i > j:
       factor = LU(i,j) / LU(j,j)
       LU(i,j) = factor                        (stored in L part)
       LU(i,k) -= factor * LU(j,k)  for k > j  (update U part)

Complexity: O(2N³/3) multiply-adds, O(N) comparisons for pivoting.

---

FAILURE MODES

• MatrixStatus::Singular — pivot magnitude ≤ tol (default PRECISION_TOLERANCE)

---