lu.cpp

Go to the documentation of this file.

#include "linalg/factorization/lu.hpp"
#include <cmath>
#include <algorithm>
 
namespace num {
 
// lu()   --  Doolittle algorithm with partial pivoting
//
// Gaussian elimination transforms A into U by subtracting multiples of pivot
// rows from rows below.  The multipliers l_{ik} = A(i,k)/A(k,k) are stored
// in the lower triangle of the working matrix (they are the entries of L).
//
// Partial pivoting: at step k, scan column k from row k downward and swap
// the largest-magnitude element into the pivot position.  This guarantees
// |L(i,j)| <= 1 for all i > j, bounding the growth factor to 2^{n-1} in
// the worst case (random matrices grow far less in practice).
//
// After n steps the working matrix holds L (below diagonal) and U (diagonal
// and above) packed together.  The diagonal of L is implicitly 1.
 
LUResult lu(const Matrix& A) {
    constexpr real singular_tol = 1e-14;
    const idx n = A.rows();
    LUResult f;
    f.LU  = A;           // working copy
    f.piv.resize(n);
    f.singular = false;
 
    Matrix& M = f.LU;
 
    for (idx k = 0; k < n; ++k) {
 
        idx  pivot_row = k;
        real pivot_val = std::abs(M(k, k));
        for (idx i = k + 1; i < n; ++i) {
            real v = std::abs(M(i, k));
            if (v > pivot_val) { pivot_val = v; pivot_row = i; }
        }
        f.piv[k] = pivot_row;
 
        if (pivot_row != k)
            for (idx j = 0; j < n; ++j)
                std::swap(M(k, j), M(pivot_row, j));
 
        if (std::abs(M(k, k)) < singular_tol) {
            f.singular = true;
            continue;
        }
 
        const real inv_ukk = real(1) / M(k, k);
        for (idx i = k + 1; i < n; ++i)
            M(i, k) *= inv_ukk;
 
        // M[i,j] -= L[i,k] * U[k,j]  for i,j > k
        for (idx i = k + 1; i < n; ++i) {
            const real lik = M(i, k);
            for (idx j = k + 1; j < n; ++j)
                M(i, j) -= lik * M(k, j);
        }
    }
 
    return f;
}
 
// lu_solve()   --  apply P, then forward/backward substitution
 
void lu_solve(const LUResult& f, const Vector& b, Vector& x) {
    const idx n    = f.LU.rows();
    const Matrix& M = f.LU;
 
    Vector y = b;   // working copy; will become x
 
    for (idx k = 0; k < n; ++k)
        if (f.piv[k] != k)
            std::swap(y[k], y[f.piv[k]]);
 
    for (idx i = 1; i < n; ++i)
        for (idx j = 0; j < i; ++j)
            y[i] -= M(i, j) * y[j];
 
    for (idx i = n; i-- > 0; ) {
        for (idx j = i + 1; j < n; ++j)
            y[i] -= M(i, j) * y[j];
        y[i] /= M(i, i);
    }
 
    x = std::move(y);
}
 
void lu_solve(const LUResult& f, const Matrix& B, Matrix& X) {
    const idx nrhs = B.cols();
    const idx n    = B.rows();
    Vector col(n), xcol(n);
    for (idx j = 0; j < nrhs; ++j) {
        for (idx i = 0; i < n; ++i) col[i] = B(i, j);
        lu_solve(f, col, xcol);
        for (idx i = 0; i < n; ++i) X(i, j) = xcol[i];
    }
}
 
// lu_det()   --  determinant from the diagonal of U and pivot parity
//
// det(A) = det(P^{-1}) * det(L) * det(U)
//        = (-1)^{swaps} * 1 * prod(U[i,i])
//
// Each transposition recorded in piv contributes a factor of -1.
// Count only entries where piv[k] != k (an identity "swap" contributes +1).
 
real lu_det(const LUResult& f) {
    const idx n = f.LU.rows();
    real det = real(1);
    for (idx i = 0; i < n; ++i)
        det *= f.LU(i, i);
    idx swaps = 0;
    for (idx k = 0; k < n; ++k)
        if (f.piv[k] != k) ++swaps;
    return (swaps % 2 == 0) ? det : -det;
}
 
// lu_inv()   --  A^{-1} by solving A * X = I column by column
 
Matrix lu_inv(const LUResult& f) {
    const idx n = f.LU.rows();
    Matrix inv(n, n, real(0));
    Vector e(n, real(0)), col(n);
    for (idx j = 0; j < n; ++j) {
        e[j] = real(1);
        lu_solve(f, e, col);
        for (idx i = 0; i < n; ++i) inv(i, j) = col[i];
        e[j] = real(0);
    }
    return inv;
}
 
} // namespace num