LU decomposition: Difference between revisions

Every square matrix <math>A</math> can be decomposed into a product of a lower triangular matrix <math>L</math> and a upper triangular matrix <math>U</math>, as described in [[wp:LU decomposition|LU decomposition]].

It is a modified form of Gaussian elimination. While the [[Cholesky decomposition]] only works for symmetric, positive definite matrices, the more general LU decomposition works for any square matrix.

There are several algorithms for calculating L and U. To derive ''Crout's algorithm'' for a 3x3 example, we have to solve the following system:

We see in the second formula that to get the <math>l_{ij}</math> below the diagonal, we have to divide by the diagonal element (pivot) <math>u_{ijjj}</math>, so we get problems when <math>u_{ijjj}</math> is either 0 or very small, which leads to numerical instability.

The solution to this problem is ''pivoting'' <math>A</math>, which means rearranging the rows of <math>A</math>, prior to the <math>LU</math> decomposition, in a way that the largest element of each column gets onto the diagonal of <math>A</math>. Rearranging the columnsrows means to multiply <math>A</math> by a permutation matrix <math>P</math>:

''';Task description''':

The task is to implement a routine which will take a square nxn matrix <math>A</math> and return a lower triangular matrix <math>L</math>, a upper triangular matrix <math>U</math> and a permutation matrix <math>P</math>, so that the above equation is fullfilled. You should then test it on the following two examples and include your output.

;Example 2:

<langsyntaxhighlight Adalang="ada">with Ada.Numerics.Generic_Real_Arrays;

end Decomposition;</langsyntaxhighlight>

<langsyntaxhighlight Adalang="ada">package body Decomposition is

end Decomposition;</langsyntaxhighlight>

<langsyntaxhighlight Adalang="ada">with Ada.Numerics.Real_Arrays;

end Decompose_Example;</langsyntaxhighlight>

<langsyntaxhighlight lang="bbcbasic"> DIM A1(2,2)

= a$</langsyntaxhighlight>

Compiled with <code>gcc -std=gnu99 -Wall -lm -pedantic</code>. Demonstrating how to do LU decomposition, and how (not) to use macros. <langsyntaxhighlight Clang="c">#include <stdio.h>

}</langsyntaxhighlight>

<langsyntaxhighlight lang="lisp">;; Creates a nxn identity matrix.

(values L U P)))</langsyntaxhighlight>

<langsyntaxhighlight lang="lisp">(setf g (make-array '(3 3) :initial-contents '((1 3 5) (2 4 7)(1 1 0))))

#2A((0 1 0) (1 0 0) (0 0 1))</langsyntaxhighlight>

<langsyntaxhighlight lang="lisp">(setf h (make-array '(4 4) :initial-contents '((11 9 24 2)(1 5 2 6)(3 17 18 1)(2 5 7 1))))

#2A((1 0 0 0) (0 0 1 0) (0 1 0 0) (0 0 0 1))</langsyntaxhighlight>

<langsyntaxhighlight lang="d">import std.stdio, std.algorithm, std.typecons, std.numeric,

std.array, std.conv, std.string, std.range, std.math;

bool isRectangular(T)(in T[][] m) pure nothrow @nogc {

bool isSquare(T)(in T[][] m) pure nothrow @nogc {

T[][] pivotize(T)(immutable T[][] m) /*pure nothrow*/

.map!((in j) => n.iota.map!(i => cast(T)(i == j)).array)

T maxm = abs(m[i][i]);

if (abs(m[j][i]) > maxm) {

lu(T)(immutable T[][] A) /*pure nothrow*/

constimmutable P = A.pivotize!T;

constimmutable A2 = matrixMul!T(P, A);

//auto f = "[%([%(%.1f, %)],\n %)]]\n\n".replicate(3);

}</langsyntaxhighlight>

=={{header|FortranEchoLisp}}==

call check( reshape([real*(8)::1,2,1,3,4,1,5,7,0 :: aa(ubound(a,1), ubound(a ],2[3,3]) )

<langsyntaxhighlight lang="go">package main

}</langsyntaxhighlight>

<langsyntaxhighlight lang="go">package main

m3 = &matrix{m2.stride, make([]float64, (len(m1.ele)/m1.stride)*m2.stride), m2.stride}

return &matrix{cols, make([]float64, rows*cols), cols}

}

}

showLU(matrixFromRows([]&matrix{3, []float64{

{1, 3, 5},

{2, 4, 7},

{1, 1, 0}}))

showLU(matrixFromRows([]&matrix{4, []float64{

{11, 9, 24, 2},

{1, 5, 2, 6},

{3, 17, 18, 1},

{2, 5, 7, 1}}))

l, u, p := a.lu()

}</langsyntaxhighlight>

}</langsyntaxhighlight>

<langsyntaxhighlight lang="j">mp=: +/ .*

LUdecompose=: (permtomat&.>@{. , }.)@:LU</langsyntaxhighlight>

<langsyntaxhighlight lang="j"> A=:3 3$1 3 5 2 4 7 1 1 0

2 5 7 1</langsyntaxhighlight>

=={{header|Mathematica}}/{{header|Wolfram Language}}==

<langsyntaxhighlight Mathematicalang="mathematica">(*Ex1*)a = {{1, 3, 5}, {2, 4, 7}, {1, 1, 0}};

<langsyntaxhighlight Matlablang="matlab"> A = [

[L,U,P] = lu(A)</langsyntaxhighlight>

<langsyntaxhighlight Matlablang="matlab"> A = [

[L,U,P] = lu(A)</langsyntaxhighlight>

<langsyntaxhighlight lang="maxima">/* LU decomposition is built-in */

/* matrix([-204], [2100], [-4740], [2940]) */</langsyntaxhighlight>

<langsyntaxhighlight PLlang="pl/Ii">(subscriptrange, fofl, size): /* 2 Nov. 2013 */

<langsyntaxhighlight lang="python">from pprint import pprint

print</langsyntaxhighlight>

<langsyntaxhighlight lang="racket">

<langsyntaxhighlight lang="rexx">/*REXX pgmprogram makescreates a matrix from console input, performs/shows LU decomposition. */

#= 0; P.= 0; PA.= 0; L.= 0; U.= 0 /*initialize some variables to 0zero. */

parse arg x /*get theobtain matrix elements from CLthe C.L. */

call makeMat call bldAMat; call showMat 'A' /*makebuild and thedisplay A matrix from numbers.*/

call showMat 'A', N call bldPmat; call showMat 'P' /*display the " " " A P matrix. " */

call manPmat call multMat; call showMat 'PA' /*manufacture " " " PA " P (permutation).*/

call showMat 'P',do y=1 for N; call bldUmat; call bldLmat /*displaybuild theU and P L matrix. " */

call showMat 'PA', N call showMat 'L'; call showMat 'U' /*display the L and PAU matrix. " */

exit do y=1 for N; call manUmat y /*manufacturestick a fork Uin it, matrix,we're all done. parts*/

call showMat 'L', N end /*display the L matrix. N*/

call showMat 'U', N if N**2\==? then do; say /'***error***display thewrong # of Uelements entered:' matrix.?; exit */9

exit /*stick a fork in it, we're done.*/end

er: say; say '***error!***' do c=1 for N; say #= # + 1; say arg_= word(1x, #); say; exitA.r.c= 13_

makeMat: ?=words(x); do N=1 for ?; if N**2==? then leave; end /*c*/

bldLmat: do do cr=1 for N; #=#+1; _=word(x,#); A /*build lower matrix.r.c=_*/

if \datatype(_,' do c=1 for N'); if r==c then calldo; L.r.c= 1; er "elementiterate; isn't numeric: " _end

end /* if c*/\==y | r==c | c>r then iterate

end /* _= PA.r*/.c

manLmat: arg ? L.r.c= _ /*manufacture L (lower) matrixU.*/c.c

do r=1 forend N /*c*/

end do c=1 for N /*r*/; if r==c then do; L.r.c=1; iterate; endreturn

bldPmat: c= N; do r=N by -1 for N; P.r.c= 1; _ c=PA c + 1 /*build perm.r matrix.c*/

if c>N then c= N%2; doif kc=1=N for then c-1; _=_-U.k.c*L.r.k; end /*k*/1

end L. /*r.c=_*/U.c.c; return

bldUmat: do r=1 endfor N; if r\==y then iterate /*rbuild upper matrix.*/

manPmat: c=N; do rk=N1 by for r-1; for N _= _ /*manufacture P- (permutation) U.k.c */ L.r.k

P.r.c=1; c=c+1; if c>N then c=N%2; if c==N then c=1 end /*k*/

end /* U.r*.c= _ / 1

multMat: do do ri=1 for N; if r\==? then iterate /*multiply matrix P and A ──► PA */

do c j=1 for N; if c<r then iterate

_ do k=PA1 for N; pa.ri.cj= (pa.i.j + p.i.k * a.k.j) / 1

do k=1 for r-1; _=_-U.k.c*L.r.k; end /*k*/

U.r end /*j*/ /*÷ by one does normalization [↑].c=_ */1

end /*ci*/; return

multMat: do i =1 for N do /*multiplyc=1 matrix for cols; w= max(w, P &length( Avalue( mat'.'r"."c ──►) PA) */)

do j =1 for N end /*c*/

do k=1 for N; pa.i.jend = (pa.i.j + p.i.k /*r* a.k.j) / 1; end

say center(mat 'matrix', end cols * (w + 1) + 7, "─") /*jdisplay the header.*/

end do /*i*/ r=1 for rows; _=

showMat: parse arg mat,rows,cols; w=0; cols=word(cols rows,1); say _

do end /*r*/; =1 for rowsreturn</syntaxhighlight>

'''{{out|output'''|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 11 9 24 2 &nbsp; &nbsp; 1 5 2 6 &nbsp; &nbsp; 3 17 18 1 &nbsp; &nbsp; 2 5 7 1 </tt>}}