Notebooks · lib/power/splu
Sparse LU with a minimum-degree ordering and threshold pivoting
splu(n, colptr, rowind, vals, q, tol) factors a square sparse matrix A,
given in compressed sparse column form, as P A Q = L U: L unit lower
triangular, U upper triangular, Q the column order q and P the row order
chosen while factoring. lu_solve(f, b) then solves A x = b.
The factorization is left-looking, column by column (Gilbert and Peierls, "Sparse partial pivoting in time proportional to arithmetic operations", SIAM Journal on Scientific and Statistical Computing 9 (1988) 862-874). Column k of L and U is the solution of a triangular system with the L found so far; its nonzeros are the rows reachable from the nonzeros of A's column through the graph of L, found by a depth-first search, so the work is proportional to the arithmetic, never to n per column.
The pivot of column k is the entry of largest magnitude among the rows not
yet pivotal (partial pivoting), except that the diagonal entry, row q[k],
is kept when its magnitude is at least tol times that largest one
(threshold pivoting, as UMFPACK and SuperLU do). With tol = 1 this is
partial pivoting; with tol = 0.001 the diagonal is kept almost always, so
the order chosen to limit fill survives.
mindeg_order(n, colptr, rowind) chooses q: the minimum-degree ordering of
the graph of A + A^T, by explicit elimination (each eliminated node's
neighbors become a clique). It is the exact minimum degree that AMD
approximates; on power-network matrices, whose fill is small, the explicit
graph stays small.
Look inside this program
- Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.
This program declares a record type (SpLU), and records run only in compiled C for now, so the Run buttons here are off. Open it in the workbench and use its C tab to compile it.
Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.