LU decomposition: Difference between revisions

=={{header|PARI/GP}}==

<lang parigp>matlup(M) =

{

my (P = matid(#M), L = P, U = M);

for (i = 1, #M-1, \\ pivoting

p = M[z=i,i];

for (k = i, #M, if (M[k,i] > p, p = M[z=k,i]));

if (i != z, \\ swap rows

k = U[i,]; U[i,] = U[z,]; U[z,] = k;

k = P[i,]; P[i,] = P[z,]; P[z,] = k;

);

);

for (i = 1, #M-1, \\ decompose

for (k = i+1, #M,

L[k,i] = U[k,i] / U[i,i];

for (j = i, #M, U[k,j] -= L[k,i] * U[i,j])

)

);

[L,U,P] \\ return L,U,P triple matrix

}</lang>

Output:

<pre>

gp > [L,U,P] = matlup([1,3,5;2,4,7;1,1,0]);

gp > L

[ 1 0 0]

[1/2 1 0]

[1/2 -1 1]

gp > U

[2 4 7]

[0 1 3/2]

[0 0 -2]

gp > P

[0 1 0]

[1 0 0]

[0 0 1]

gp > [L,U,P] = matlup([11,9,24,2;1,5,2,6;3,17,18,1;2,5,7,1]);

gp > L

[ 1 0 0 0]

[3/11 1 0 0]

[1/11 23/80 1 0]

[2/11 37/160 1/278 1]

gp > U

[11 9 24 2]

[ 0 160/11 126/11 5/11]

[ 0 0 -139/40 91/16]

[ 0 0 0 71/139]

gp > P

[1 0 0 0]

[0 0 1 0]

[0 1 0 0]

[0 0 0 1]

</pre>