LU decomposition: Difference between revisions

</pre>

=={{header|zkl}}==

{{trans|Common Lisp}}

{{trans|D}} to get the conversion from column major to row major correct.

A matrix is a list of lists, ie list of rows in row major order.

<lang zkl>fcn make_array(n,m,v){ (m).pump(List.createLong(m).write,v)*n }

fcn eye(n){ // Creates a nxn identity matrix.

I:=make_array(n,n,0.0);

foreach j in (n){ I[j][j]=1.0 }

I

}

// Creates the pivoting matrix for A.

fcn pivotize(A){

n:=A.len(); // rows

P:=eye(n);

foreach i in (n){

max,row:=A[i][i],i;

foreach j in ([i..n-1]){

if(A[j][i]>max) max,row=A[j][i],j;

}

if(i!=row) P.swap(i,row);

}

// Return P.

P

}

// Decomposes a square matrix A by PA=LU and returns L, U and P.

fcn lu(A){

n:=A.len();

L:=make_array(n,n,0.0);

U:=make_array(n,n,0.0);

P:=pivotize(A);

A=matMult(P,A);

foreach j in (n){

L[j][j]=1.0;

foreach i in (j+1){

U[i][j]=A[i][j] -

(i).reduce('wrap(s,k){ s + U[k][j]*L[i][k] },0.0);

}

foreach i in ([j..n-1]){

L[i][j]=( A[i][j] -

(j).reduce('wrap(s,k){ s + U[k][j]*L[i][k] },0.0) ) /

U[j][j];

}

}

// Return L, U and P.

return(L,U,P);

}

fcn matMult(a,b){

n,m,p:=a[0].len(),a.len(),b[0].len();

ans:=make_array(n,m,0.0);

foreach i,j,k in (m,p,n){ ans[i][j]+=a[i][k]*b[k][j]; }

ans

}</lang>

Example 1

<lang zkl>g:=L(L(1.0,3.0,5.0),L(2.0,4.0,7.0),L(1.0,1.0,0.0));

lu(g).apply2("println");</lang>

{{out}}

<pre>

L(L(1,0,0),L(0.5,1,0),L(0.5,-1,1))

L(L(2,4,7),L(0,1,1.5),L(0,0,-2))

L(L(0,1,0),L(1,0,0),L(0,0,1))

</pre>

Example 2

<lang zkl>lu(L( L(11.0, 9.0, 24.0, 2.0),

L( 1.0, 5.0, 2.0, 6.0),

L( 3.0, 17.0, 18.0, 1.0),

L( 2.0, 5.0, 7.0, 1.0) )).apply2(T(printM,Console.writeln.fpM("-")));

fcn printM(m) { m.pump(Console.println,rowFmt) }

fcn rowFmt(row){ ("%9.5f "*row.len()).fmt(row.xplode()) }</lang>

The list apply2 method is side effects only, it doesn't aggregate results. When given a list of actions, it applies the action and passes the result to the next action. The fpM method is partial application with a mask, "-" truncates the parameters at that point (in this case, no parameters, ie just print a blank line, not the result of printM).

{{out}}

<pre>

1.00000 0.00000 0.00000 0.00000

0.27273 1.00000 0.00000 0.00000

0.09091 0.28750 1.00000 0.00000

0.18182 0.23125 0.00360 1.00000

11.00000 9.00000 24.00000 2.00000

0.00000 14.54545 11.45455 0.45455

0.00000 0.00000 -3.47500 5.68750

0.00000 0.00000 0.00000 0.51079

1.00000 0.00000 0.00000 0.00000

0.00000 0.00000 1.00000 0.00000

0.00000 1.00000 0.00000 0.00000

0.00000 0.00000 0.00000 1.00000