Computes the LU decomposition without pivoting.
[L, U] = linalg_factorlu ( A ) [L, U] = linalg_factorlu ( A , verbose )
a n-by-n matrix of doubles
a 1-by-1 matrix of boolean (default verbose=%f), set to %t to display intermediate messages
a n-by-n matrix of doubles, lower triangular.
a n-by-n matrix of doubles, upper triangular.
Decomposes A as A = L*U with L lower triangular and unit diagonal, U upper triangular and non-unit diagonal. This algorithm might be sensitive to inaccuracies if A near singular. There is no solution is A is singular. The algorithm might be unstable, even if A is not ill-conditionned. Indeed, it might produce large entries in L and U, even if entries in A are of moderate size. Uses an effective algorithm with vectorization.
// Basic test A=[ 1 1 1 1 1 1 2 2 2 2 1 2 3 3 3 1 2 3 4 4 1 2 3 4 5 ]; [L,U]=linalg_factorlu(A) Lexpected=[ 1 0 0 0 0 1 1 0 0 0 1 1 1 0 0 1 1 1 1 0 1 1 1 1 1 ]; Uexpected=[ 1 1 1 1 1 0 1 1 1 1 0 0 1 1 1 0 0 0 1 1 0 0 0 0 1 ]; // See what happens [L,U]=linalg_factorlu(A,%t) // Combine factorlu and solvelu A=[ 1 1 1 1 1 1 2 2 2 2 1 2 3 3 3 1 2 3 4 4 1 2 3 4 5 ]; b=[5;9;12;14;15]; [L,U]=linalg_factorlu(A) x=linalg_solvelu(L,U,b,%f) xexpected=[1;1;1;1;1] // A difficult case: pivoting is necessary A=[ 1.d-8 1 1 1 ]; [L,U]=linalg_factorlu(A) Lexpected=[ 1. 0 1.d8 1.] Uexpected=[ 1.d-8 1 0. -1.d8+1 ] cond(A) // 2.6180340238745070102766 cond(L) // 9999999949752408. cond(U) // 9999999900000000. // See the algorithm edit linalg_factorlu | ![]() | ![]() |