blob: 6717733046d8ef87ff64c9c6efc50d0d040e7ff1 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
|
#pragma once
#include <cstddef>
#include <vector>
#include "matrix.hpp"
#include "vector.hpp"
namespace linalg {
// Result of LU factorization with partial pivoting.
//
// The factorization satisfies PA = LU, where:
// P is the permutation matrix encoded by `perm`
// L is unit lower triangular (L[i][i] == 1)
// U is upper triangular
//
// `perm[i]` = index of the original row that ended up at position i.
// Applying P to a vector b means: (Pb)[i] = b[perm[i]].
//
// `sign` is the sign of the permutation: +1 if an even number of row
// swaps were made, -1 if odd. Useful for computing det(A) = sign * prod(diag(U)).
struct LUResult {
Matrix L;
Matrix U;
std::vector<std::size_t> perm;
int sign;
};
// Compute the LU factorization of A with partial pivoting.
//
// Throws DimensionMismatchError if A is not square.
// Throws SingularMatrixError if A is (numerically) singular, i.e. any
// pivot is smaller in magnitude than `singular_tolerance`.
LUResult lu_factor(const Matrix& A, double singular_tolerance = 1e-12);
// Solve Ax = b given a precomputed LU factorization.
//
// Applies the stored permutation, then forward / backward substitution.
// Throws DimensionMismatchError if b.size() != lu.L.rows().
Vector lu_solve(const LUResult& lu, const Vector& b);
} // namespace linalg
|