LU decomposition: Difference between revisions
Added Sidef |
|||
Line 2,164:
lu_backsub(lup, transpose([1, 1, -1, -1]));
/* matrix([-204], [2100], [-4740], [2940]) */</lang>
=={{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>
=={{header|Perl 6}}==
|
Revision as of 12:58, 29 June 2016
You are encouraged to solve this task according to the task description, using any language you may know.
Every square matrix can be decomposed into a product of a lower triangular matrix and a upper triangular matrix , as described in 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 now would have to solve 9 equations with 12 unknowns. To make the system uniquely solvable, usually the diagonal elements of are set to 1
so we get a solvable system of 9 unknowns and 9 equations.
Solving for the other and , we get the following equations:
and for :
We see that there is a calculation pattern, which can be expressed as the following formulas, first for
and then for
We see in the second formula that to get the below the diagonal, we have to divide by the diagonal element (pivot) , so we get problems when is either 0 or very small, which leads to numerical instability.
The solution to this problem is pivoting , which means rearranging the rows of , prior to the decomposition, in a way that the largest element of each column gets onto the diagonal of . Rearranging the rows means to multiply by a permutation matrix :
Example:
The decomposition algorithm is then applied on the rearranged matrix so that
Task description
The task is to implement a routine which will take a square nxn matrix and return a lower triangular matrix , a upper triangular matrix and a permutation matrix , so that the above equation is fullfilled. You should then test it on the following two examples and include your output.
Example 1:
A 1 3 5 2 4 7 1 1 0 L 1.00000 0.00000 0.00000 0.50000 1.00000 0.00000 0.50000 -1.00000 1.00000 U 2.00000 4.00000 7.00000 0.00000 1.00000 1.50000 0.00000 0.00000 -2.00000 P 0 1 0 1 0 0 0 0 1
Example 2:
A 11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 L 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 U 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 P 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1
Ada
decomposition.ads: <lang Ada>with Ada.Numerics.Generic_Real_Arrays; generic
with package Matrix is new Ada.Numerics.Generic_Real_Arrays (<>);
package Decomposition is
-- decompose a square matrix A by PA = LU procedure Decompose (A : Matrix.Real_Matrix; P, L, U : out Matrix.Real_Matrix);
end Decomposition;</lang>
decomposition.adb: <lang Ada>package body Decomposition is
procedure Swap_Rows (M : in out Matrix.Real_Matrix; From, To : Natural) is Temporary : Matrix.Real; begin if From = To then return; end if; for I in M'Range (2) loop Temporary := M (M'First (1) + From, I); M (M'First (1) + From, I) := M (M'First (1) + To, I); M (M'First (1) + To, I) := Temporary; end loop; end Swap_Rows;
function Pivoting_Matrix (M : Matrix.Real_Matrix) return Matrix.Real_Matrix is use type Matrix.Real; Order : constant Positive := M'Length (1); Result : Matrix.Real_Matrix := Matrix.Unit_Matrix (Order); Max : Matrix.Real; Row : Natural; begin for J in 0 .. Order - 1 loop Max := M (M'First (1) + J, M'First (2) + J); Row := J; for I in J .. Order - 1 loop if M (M'First (1) + I, M'First (2) + J) > Max then Max := M (M'First (1) + I, M'First (2) + J); Row := I; end if; end loop; if J /= Row then -- swap rows J and Row Swap_Rows (Result, J, Row); end if; end loop; return Result; end Pivoting_Matrix;
procedure Decompose (A : Matrix.Real_Matrix; P, L, U : out Matrix.Real_Matrix) is use type Matrix.Real_Matrix, Matrix.Real; Order : constant Positive := A'Length (1); A2 : Matrix.Real_Matrix (A'Range (1), A'Range (2)); S : Matrix.Real; begin L := (others => (others => 0.0)); U := (others => (others => 0.0)); P := Pivoting_Matrix (A); A2 := P * A; for J in 0 .. Order - 1 loop L (L'First (1) + J, L'First (2) + J) := 1.0; for I in 0 .. J loop S := 0.0; for K in 0 .. I - 1 loop S := S + U (U'First (1) + K, U'First (2) + J) * L (L'First (1) + I, L'First (2) + K); end loop; U (U'First (1) + I, U'First (2) + J) := A2 (A2'First (1) + I, A2'First (2) + J) - S; end loop; for I in J + 1 .. Order - 1 loop S := 0.0; for K in 0 .. J loop S := S + U (U'First (1) + K, U'First (2) + J) * L (L'First (1) + I, L'First (2) + K); end loop; L (L'First (1) + I, L'First (2) + J) := (A2 (A2'First (1) + I, A2'First (2) + J) - S) / U (U'First (1) + J, U'First (2) + J); end loop; end loop; end Decompose;
end Decomposition;</lang>
Example usage: <lang Ada>with Ada.Numerics.Real_Arrays; with Ada.Text_IO; with Decomposition; procedure Decompose_Example is
package Real_Decomposition is new Decomposition (Matrix => Ada.Numerics.Real_Arrays);
package Real_IO is new Ada.Text_IO.Float_IO (Float); procedure Print (M : Ada.Numerics.Real_Arrays.Real_Matrix) is begin for Row in M'Range (1) loop for Col in M'Range (2) loop Real_IO.Put (M (Row, Col), 3, 2, 0); end loop; Ada.Text_IO.New_Line; end loop; end Print;
Example_1 : constant Ada.Numerics.Real_Arrays.Real_Matrix := ((1.0, 3.0, 5.0), (2.0, 4.0, 7.0), (1.0, 1.0, 0.0)); P_1, L_1, U_1 : Ada.Numerics.Real_Arrays.Real_Matrix (Example_1'Range (1), Example_1'Range (2)); Example_2 : constant Ada.Numerics.Real_Arrays.Real_Matrix := ((11.0, 9.0, 24.0, 2.0), (1.0, 5.0, 2.0, 6.0), (3.0, 17.0, 18.0, 1.0), (2.0, 5.0, 7.0, 1.0)); P_2, L_2, U_2 : Ada.Numerics.Real_Arrays.Real_Matrix (Example_2'Range (1), Example_2'Range (2));
begin
Real_Decomposition.Decompose (A => Example_1, P => P_1, L => L_1, U => U_1); Real_Decomposition.Decompose (A => Example_2, P => P_2, L => L_2, U => U_2); Ada.Text_IO.Put_Line ("Example 1:"); Ada.Text_IO.Put_Line ("A:"); Print (Example_1); Ada.Text_IO.Put_Line ("L:"); Print (L_1); Ada.Text_IO.Put_Line ("U:"); Print (U_1); Ada.Text_IO.Put_Line ("P:"); Print (P_1); Ada.Text_IO.New_Line; Ada.Text_IO.Put_Line ("Example 2:"); Ada.Text_IO.Put_Line ("A:"); Print (Example_2); Ada.Text_IO.Put_Line ("L:"); Print (L_2); Ada.Text_IO.Put_Line ("U:"); Print (U_2); Ada.Text_IO.Put_Line ("P:"); Print (P_2);
end Decompose_Example;</lang>
- Output:
Example 1: A: 1.00 3.00 5.00 2.00 4.00 7.00 1.00 1.00 0.00 L: 1.00 0.00 0.00 0.50 1.00 0.00 0.50 -1.00 1.00 U: 2.00 4.00 7.00 0.00 1.00 1.50 0.00 0.00 -2.00 P: 0.00 1.00 0.00 1.00 0.00 0.00 0.00 0.00 1.00 Example 2: A: 11.00 9.00 24.00 2.00 1.00 5.00 2.00 6.00 3.00 17.00 18.00 1.00 2.00 5.00 7.00 1.00 L: 1.00 0.00 0.00 0.00 0.27 1.00 0.00 0.00 0.09 0.29 1.00 0.00 0.18 0.23 0.00 1.00 U: 11.00 9.00 24.00 2.00 0.00 14.55 11.45 0.45 0.00 0.00 -3.47 5.69 0.00 0.00 0.00 0.51 P: 1.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00
BBC BASIC
<lang bbcbasic> DIM A1(2,2)
A1() = 1, 3, 5, 2, 4, 7, 1, 1, 0 PROCLUdecomposition(A1(), L1(), U1(), P1()) PRINT "L1:" ' FNshowmatrix(L1()) PRINT "U1:" ' FNshowmatrix(U1()) PRINT "P1:" ' FNshowmatrix(P1()) DIM A2(3,3) A2() = 11, 9, 24, 2, 1, 5, 2, 6, 3, 17, 18, 1, 2, 5, 7, 1 PROCLUdecomposition(A2(), L2(), U2(), P2()) PRINT "L2:" ' FNshowmatrix(L2()) PRINT "U2:" ' FNshowmatrix(U2()) PRINT "P2:" ' FNshowmatrix(P2()) END DEF PROCLUdecomposition(a(), RETURN l(), RETURN u(), RETURN p()) LOCAL i%, j%, k%, n%, s, b() : n% = DIM(a(),2) DIM l(n%,n%), u(n%,n%), b(n%,n%) PROCpivot(a(), p()) b() = p() . a() FOR j% = 0 TO n% l(j%,j%) = 1 FOR i% = 0 TO j% s = 0 FOR k% = 0 TO i% : s += u(k%,j%) * l(i%,k%) : NEXT u(i%,j%) = b(i%,j%) - s NEXT FOR i% = j% TO n% s = 0 FOR k% = 0 TO j% : s += u(k%,j%) * l(i%,k%) : NEXT IF i%<>j% l(i%,j%) = (b(i%,j%) - s) / u(j%,j%) NEXT NEXT j% ENDPROC DEF PROCpivot(a(), RETURN p()) LOCAL i%, j%, m%, n%, r% : n% = DIM(a(),2) DIM p(n%,n%) : FOR i% = 0 TO n% : p(i%,i%) = 1 : NEXT FOR i% = 0 TO n% m% = a(i%,i%) r% = i% FOR j% = i% TO n% IF a(j%,i%) > m% m% = a(j%,i%) : r% = j% NEXT IF i%<>r% THEN FOR j% = 0 TO n% : SWAP p(i%,j%),p(r%,j%) : NEXT ENDIF NEXT i% ENDPROC DEF FNshowmatrix(a()) LOCAL @%, i%, j%, a$ @% = &102050A FOR i% = 0 TO DIM(a(),1) FOR j% = 0 TO DIM(a(),2) a$ += STR$(a(i%,j%)) + ", " NEXT a$ = LEFT$(LEFT$(a$)) + CHR$(13) + CHR$(10) NEXT i% = a$</lang>
- Output:
L1: 1.00000, 0.00000, 0.00000 0.50000, 1.00000, 0.00000 0.50000, -1.00000, 1.00000 U1: 2.00000, 4.00000, 7.00000 0.00000, 1.00000, 1.50000 0.00000, 0.00000, -2.00000 P1: 0.00000, 1.00000, 0.00000 1.00000, 0.00000, 0.00000 0.00000, 0.00000, 1.00000 L2: 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 U2: 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 P2: 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
C
Compiled with gcc -std=gnu99 -Wall -lm -pedantic
. Demonstrating how to do LU decomposition, and how (not) to use macros. <lang C>#include <stdio.h>
- include <stdlib.h>
- include <math.h>
- define foreach(a, b, c) for (int a = b; a < c; a++)
- define for_i foreach(i, 0, n)
- define for_j foreach(j, 0, n)
- define for_k foreach(k, 0, n)
- define for_ij for_i for_j
- define for_ijk for_ij for_k
- define _dim int n
- define _swap(x, y) { typeof(x) tmp = x; x = y; y = tmp; }
- define _sum_k(a, b, c, s) { s = 0; foreach(k, a, b) s+= c; }
typedef double **mat;
- define _zero(a) mat_zero(a, n)
void mat_zero(mat x, int n) { for_ij x[i][j] = 0; }
- define _new(a) a = mat_new(n)
mat mat_new(_dim) { mat x = malloc(sizeof(double*) * n); x[0] = malloc(sizeof(double) * n * n);
for_i x[i] = x[0] + n * i; _zero(x);
return x; }
- define _copy(a) mat_copy(a, n)
mat mat_copy(void *s, _dim) { mat x = mat_new(n); for_ij x[i][j] = ((double (*)[n])s)[i][j]; return x; }
- define _del(x) mat_del(x)
void mat_del(mat x) { free(x[0]); free(x); }
- define _QUOT(x) #x
- define QUOTE(x) _QUOT(x)
- define _show(a) printf(QUOTE(a)" =");mat_show(a, 0, n)
void mat_show(mat x, char *fmt, _dim) { if (!fmt) fmt = "%8.4g"; for_i { printf(i ? " " : " [ "); for_j { printf(fmt, x[i][j]); printf(j < n - 1 ? " " : i == n - 1 ? " ]\n" : "\n"); } } }
- define _mul(a, b) mat_mul(a, b, n)
mat mat_mul(mat a, mat b, _dim) { mat c = _new(c); for_ijk c[i][j] += a[i][k] * b[k][j]; return c; }
- define _pivot(a, b) mat_pivot(a, b, n)
void mat_pivot(mat a, mat p, _dim) { for_ij { p[i][j] = (i == j); } for_i { int max_j = i; foreach(j, i, n) if (fabs(a[j][i]) > fabs(a[max_j][i])) max_j = j;
if (max_j != i) for_k { _swap(p[i][k], p[max_j][k]); } } }
- define _LU(a, l, u, p) mat_LU(a, l, u, p, n)
void mat_LU(mat A, mat L, mat U, mat P, _dim) { _zero(L); _zero(U); _pivot(A, P);
mat Aprime = _mul(P, A);
for_i { L[i][i] = 1; } for_ij { double s; if (j <= i) { _sum_k(0, j, L[j][k] * U[k][i], s) U[j][i] = Aprime[j][i] - s; } if (j >= i) { _sum_k(0, i, L[j][k] * U[k][i], s); L[j][i] = (Aprime[j][i] - s) / U[i][i]; } }
_del(Aprime); }
double A3[][3] = {{ 1, 3, 5 }, { 2, 4, 7 }, { 1, 1, 0 }}; double A4[][4] = {{11, 9, 24, 2}, {1, 5, 2, 6}, {3, 17, 18, 1}, {2, 5, 7, 1}};
int main() { int n = 3; mat A, L, P, U;
_new(L); _new(P); _new(U); A = _copy(A3); _LU(A, L, U, P); _show(A); _show(L); _show(U); _show(P); _del(A); _del(L); _del(U); _del(P);
printf("\n");
n = 4;
_new(L); _new(P); _new(U); A = _copy(A4); _LU(A, L, U, P); _show(A); _show(L); _show(U); _show(P); _del(A); _del(L); _del(U); _del(P);
return 0; }</lang>
Common Lisp
Uses the routine (mmul A B) from Matrix multiplication.
<lang lisp>;; Creates a nxn identity matrix. (defun eye (n)
(let ((I (make-array `(,n ,n) :initial-element 0))) (loop for j from 0 to (- n 1) do (setf (aref I j j) 1)) I))
- Swap two rows l and k of a mxn matrix A, which is a 2D array.
(defun swap-rows (A l k)
(let* ((n (cadr (array-dimensions A))) (row (make-array n :initial-element 0))) (loop for j from 0 to (- n 1) do (setf (aref row j) (aref A l j)) (setf (aref A l j) (aref A k j)) (setf (aref A k j) (aref row j)))))
- Creates the pivoting matrix for A.
(defun pivotize (A)
(let* ((n (car (array-dimensions A))) (P (eye n))) (loop for j from 0 to (- n 1) do (let ((max (aref A j j)) (row j)) (loop for i from j to (- n 1) do (if (> (aref A i j) max) (setq max (aref A i j) row i))) (if (not (= j row)) (swap-rows P j row))))
;; Return P. P))
- Decomposes a square matrix A by PA=LU and returns L, U and P.
(defun lu (A)
(let* ((n (car (array-dimensions A))) (L (make-array `(,n ,n) :initial-element 0)) (U (make-array `(,n ,n) :initial-element 0)) (P (pivotize A)) (A (mmul P A)))
(loop for j from 0 to (- n 1) do (setf (aref L j j) 1) (loop for i from 0 to j do (setf (aref U i j) (- (aref A i j) (loop for k from 0 to (- i 1) sum (* (aref U k j) (aref L i k)))))) (loop for i from j to (- n 1) do (setf (aref L i j) (/ (- (aref A i j) (loop for k from 0 to (- j 1) sum (* (aref U k j) (aref L i k)))) (aref U j j)))))
;; Return L, U and P. (values L U P)))</lang>
Example 1:
<lang lisp>(setf g (make-array '(3 3) :initial-contents '((1 3 5) (2 4 7)(1 1 0))))
- 2A((1 3 5) (2 4 7) (1 1 0))
(lu g)
- 2A((1 0 0) (1/2 1 0) (1/2 -1 1))
- 2A((2 4 7) (0 1 3/2) (0 0 -2))
- 2A((0 1 0) (1 0 0) (0 0 1))</lang>
Example 2:
<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((11 9 24 2) (1 5 2 6) (3 17 18 1) (2 5 7 1))
(lup h)
- 2A((1 0 0 0) (3/11 1 0 0) (1/11 23/80 1 0) (2/11 37/160 1/278 1))
- 2A((11 9 24 2) (0 160/11 126/11 5/11) (0 0 -139/40 91/16) (0 0 0 71/139))
- 2A((1 0 0 0) (0 0 1 0) (0 1 0 0) (0 0 0 1))</lang>
D
<lang d>import std.stdio, std.algorithm, std.typecons, std.numeric,
std.array, std.conv, std.string, std.range;
bool isRectangular(T)(in T[][] m) pure nothrow @nogc {
return m.all!(r => r.length == m[0].length);
}
bool isSquare(T)(in T[][] m) pure nothrow @nogc {
return m.isRectangular && m[0].length == m.length;
}
T[][] matrixMul(T)(in T[][] A, in T[][] B) pure nothrow in {
assert(A.isRectangular && B.isRectangular && !A.empty && !B.empty && A[0].length == B.length);
} body {
auto result = new T[][](A.length, B[0].length); auto aux = new T[B.length];
foreach (immutable j; 0 .. B[0].length) { foreach (immutable k, const row; B) aux[k] = row[j]; foreach (immutable i, const ai; A) result[i][j] = dotProduct(ai, aux); }
return result;
}
/// Creates the pivoting matrix for m. T[][] pivotize(T)(immutable T[][] m) pure nothrow in {
assert(m.isSquare);
} body {
immutable n = m.length; auto id = iota(n) .map!((in j) => n.iota.map!(i => T(i == j)).array) .array;
foreach (immutable i; 0 .. n) { // immutable row = iota(i, n).reduce!(max!(j => m[j][i])); T maxm = m[i][i]; size_t row = i; foreach (immutable j; i .. n) if (m[j][i] > maxm) { maxm = m[j][i]; row = j; }
if (i != row) swap(id[i], id[row]); }
return id;
}
/// Decomposes a square matrix A by PA=LU and returns L, U and P. Tuple!(T[][],"L", T[][],"U", const T[][],"P") lu(T)(immutable T[][] A) pure nothrow in {
assert(A.isSquare);
} body {
immutable n = A.length; auto L = new T[][](n, n); auto U = new T[][](n, n); foreach (immutable i; 0 .. n) { L[i][i .. $] = 0; U[i][0 .. i] = 0; }
immutable P = A.pivotize!T; immutable A2 = matrixMul!T(P, A);
foreach (immutable j; 0 .. n) { L[j][j] = 1; foreach (immutable i; 0 .. j+1) { T s1 = 0; foreach (immutable k; 0 .. i) s1 += U[k][j] * L[i][k]; U[i][j] = A2[i][j] - s1; } foreach (immutable i; j .. n) { T s2 = 0; foreach (immutable k; 0 .. j) s2 += U[k][j] * L[i][k]; L[i][j] = (A2[i][j] - s2) / U[j][j]; } }
return typeof(return)(L, U, P);
}
void main() {
immutable a = [[1.0, 3, 5], [2.0, 4, 7], [1.0, 1, 0]]; immutable b = [[11.0, 9, 24, 2], [1.0, 5, 2, 6], [3.0, 17, 18, 1], [2.0, 5, 7, 1]];
auto f = "[%([%(%.1f, %)],\n %)]]\n\n".replicate(3); foreach (immutable m; [a, b]) writefln(f, lu(m).tupleof);
}</lang>
- Output:
[[1.0, 0.0, 0.0], [0.5, 1.0, 0.0], [0.5, -1.0, 1.0]] [[2.0, 4.0, 7.0], [0.0, 1.0, 1.5], [0.0, 0.0, -2.0]] [[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]] [[1.0, 0.0, 0.0, 0.0], [0.3, 1.0, 0.0, 0.0], [0.0, 0.3, 1.0, 0.0], [0.2, 0.2, 0.0, 1.0]] [[11.0, 9.0, 24.0, 2.0], [0.0, 14.5, 11.5, 0.5], [0.0, 0.0, -3.5, 5.7], [0.0, 0.0, 0.0, 0.5]] [[1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
EchoLisp
<lang scheme> (lib 'matrix) ;; the matrix library provides LU-decomposition (decimals 5)
(define A (list->array' (1 3 5 2 4 7 1 1 0 ) 3 3)) (define PLU (matrix-lu-decompose A)) ;; -> list of three matrices, P, Lower, Upper
(array-print (first PLU)) 0 1 0 1 0 0 0 0 1 (array-print (second PLU)) 1 0 0 0.5 1 0 0.5 -1 1 (array-print (caddr PLU)) 2 4 7 0 1 1.5 0 0 -2
(define A (list->array '(11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 ) 4 4)) (define PLU (matrix-lu-decompose A)) ;; -> list of three matrices, P, Lower, Upper (array-print (first PLU)) 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1
(array-print (second PLU)) 1 0 0 0 0.27273 1 0 0 0.09091 0.2875 1 0 0.18182 0.23125 0.0036 1
(array-print (caddr PLU)) 11 9 24 2 0 14.54545 11.45455 0.45455 0 0 -3.475 5.6875 0 0 0 0.51079 </lang>
Fortran
<lang Fortran> program lu1
implicit none call check( reshape([real(8)::1,2,1,3,4,1,5,7,0 ],[3,3]) ) call check( reshape([real(8)::11,1,3,2,9,5,17,5,24,2,18,7,2,6,1,1],[4,4]) )
contains
subroutine check(a) real(8), intent(in) :: a(:,:) integer :: i,j,n real(8), allocatable :: aa(:,:),l(:,:),u(:,:) integer, allocatable :: p(:,:) integer, allocatable :: ipiv(:) n = size(a,1) allocate(aa(n,n),l(n,n),u(n,n),p(n,n),ipiv(n)) forall (j=1:n,i=1:n) aa(i,j) = a(i,j) u (i,j) = 0d0 p (i,j) = merge(1 ,0 ,i.eq.j) l (i,j) = merge(1d0,0d0,i.eq.j) end forall call lu(aa, ipiv) do i = 1,n l(i, :i-1) = aa(ipiv(i), :i-1) u(i,i: ) = aa(ipiv(i),i: ) end do p(ipiv,:) = p call mat_print('a',a) call mat_print('p',p) call mat_print('l',l) call mat_print('u',u) print *, "residual" print *, "|| P.A - L.U || = ", maxval(abs(matmul(p,a)-matmul(l,u))) end subroutine subroutine lu(a,p)
! in situ decomposition, corresponds to LAPACK's dgebtrf
real(8), intent(inout) :: a(:,:) integer, intent(out ) :: p(:) integer :: n, i,j,k,kmax n = size(a,1) p = [ ( i, i=1,n ) ] do k = 1,n-1 kmax = maxloc(abs(a(p(k:),k)),1) + k-1 if (kmax /= k ) p([k, kmax]) = p([kmax, k]) a(p(k+1:),k) = a(p(k+1:),k) / a(p(k),k) forall (j=k+1:n) a(p(k+1:),j) = a(p(k+1:),j) - a(p(k+1:),k) * a(p(k),j) end do end subroutine
subroutine mat_print(amsg,a) character(*), intent(in) :: amsg class (*), intent(in) :: a(:,:) integer :: i print*,' ' print*,amsg do i=1,size(a,1) select type (a) type is (real(8)) ; print'(100f8.2)',a(i,:) type is (integer) ; print'(100i8 )',a(i,:) end select end do print*,' ' end subroutine
end program </lang>
- Output:
a 1.00 3.00 5.00 2.00 4.00 7.00 1.00 1.00 0.00 p 0.00 1.00 0.00 1.00 0.00 0.00 0.00 0.00 1.00 l 1.00 0.00 0.00 0.50 1.00 0.00 0.50 -1.00 1.00 u 2.00 4.00 7.00 0.00 1.00 1.50 0.00 0.00 -2.00 residual || P.A - L.U || = 0.0000000000000000 a 11.00 9.00 24.00 2.00 1.00 5.00 2.00 6.00 3.00 17.00 18.00 1.00 2.00 5.00 7.00 1.00 p 1.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 l 1.00 0.00 0.00 0.00 0.27 1.00 0.00 0.00 0.09 0.29 1.00 0.00 0.18 0.23 0.00 1.00 u 11.00 9.00 24.00 2.00 0.00 14.55 11.45 0.45 0.00 0.00 -3.47 5.69 0.00 0.00 0.00 0.51 residual || P.A - L.U || = 0.0000000000000000
Go
2D representation
<lang go>package main
import "fmt"
type matrix [][]float64
func zero(n int) matrix {
r := make([][]float64, n) a := make([]float64, n*n) for i := range r { r[i] = a[n*i : n*(i+1)] } return r
}
func eye(n int) matrix {
r := zero(n) for i := range r { r[i][i] = 1 } return r
}
func (m matrix) print(label string) {
if label > "" { fmt.Printf("%s:\n", label) } for _, r := range m { for _, e := range r { fmt.Printf(" %9.5f", e) } fmt.Println() }
}
func (a matrix) pivotize() matrix {
p := eye(len(a)) for j, r := range a { max := r[j] row := j for i := j; i < len(a); i++ { if a[i][j] > max { max = a[i][j] row = i } } if j != row { // swap rows p[j], p[row] = p[row], p[j] } } return p
}
func (m1 matrix) mul(m2 matrix) matrix {
r := zero(len(m1)) for i, r1 := range m1 { for j := range m2 { for k := range m1 { r[i][j] += r1[k] * m2[k][j] } } } return r
}
func (a matrix) lu() (l, u, p matrix) {
l = zero(len(a)) u = zero(len(a)) p = a.pivotize() a = p.mul(a) for j := range a { l[j][j] = 1 for i := 0; i <= j; i++ { sum := 0. for k := 0; k < i; k++ { sum += u[k][j] * l[i][k] } u[i][j] = a[i][j] - sum } for i := j; i < len(a); i++ { sum := 0. for k := 0; k < j; k++ { sum += u[k][j] * l[i][k] } l[i][j] = (a[i][j] - sum) / u[j][j] } } return
}
func main() {
showLU(matrix{ {1, 3, 5}, {2, 4, 7}, {1, 1, 0}}) showLU(matrix{ {11, 9, 24, 2}, {1, 5, 2, 6}, {3, 17, 18, 1}, {2, 5, 7, 1}})
}
func showLU(a matrix) {
a.print("\na") l, u, p := a.lu() l.print("l") u.print("u") p.print("p")
}</lang>
- Output:
a: 1.00000 3.00000 5.00000 2.00000 4.00000 7.00000 1.00000 1.00000 0.00000 l: 1.00000 0.00000 0.00000 0.50000 1.00000 0.00000 0.50000 -1.00000 1.00000 u: 2.00000 4.00000 7.00000 0.00000 1.00000 1.50000 0.00000 0.00000 -2.00000 p: 0.00000 1.00000 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 1.00000 a: 11.00000 9.00000 24.00000 2.00000 1.00000 5.00000 2.00000 6.00000 3.00000 17.00000 18.00000 1.00000 2.00000 5.00000 7.00000 1.00000 l: 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 u: 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 p: 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
Flat representation
<lang go>package main
import "fmt"
type matrix struct {
stride int ele []float64
}
func (m *matrix) print(heading string) {
if heading > "" { fmt.Print("\n", heading, "\n") } for e := 0; e < len(m.ele); e += m.stride { fmt.Printf("%8.5f ", m.ele[e:e+m.stride]) fmt.Println() }
}
func (m1 *matrix) mul(m2 *matrix) (m3 *matrix, ok bool) {
if m1.stride*m2.stride != len(m2.ele) { return nil, false } m3 = &matrix{m2.stride, make([]float64, (len(m1.ele)/m1.stride)*m2.stride)} for m1c0, m3x := 0, 0; m1c0 < len(m1.ele); m1c0 += m1.stride { for m2r0 := 0; m2r0 < m2.stride; m2r0++ { for m1x, m2x := m1c0, m2r0; m2x < len(m2.ele); m2x += m2.stride { m3.ele[m3x] += m1.ele[m1x] * m2.ele[m2x] m1x++ } m3x++ } } return m3, true
}
func zero(rows, cols int) *matrix {
return &matrix{cols, make([]float64, rows*cols)}
}
func eye(n int) *matrix {
m := zero(n, n) for ix := 0; ix < len(m.ele); ix += n + 1 { m.ele[ix] = 1 } return m
}
func (a *matrix) pivotize() *matrix {
pv := make([]int, a.stride) for i := range pv { pv[i] = i } for j, dx := 0, 0; j < a.stride; j++ { row := j max := a.ele[dx] for i, ixcj := j, dx; i < a.stride; i++ { if a.ele[ixcj] > max { max = a.ele[ixcj] row = i } ixcj += a.stride } if j != row { pv[row], pv[j] = pv[j], pv[row] } dx += a.stride + 1 } p := zero(a.stride, a.stride) for r, c := range pv { p.ele[r*a.stride+c] = 1 } return p
}
func (a *matrix) lu() (l, u, p *matrix) {
l = zero(a.stride, a.stride) u = zero(a.stride, a.stride) p = a.pivotize() a, _ = p.mul(a) for j, jxc0 := 0, 0; j < a.stride; j++ { l.ele[jxc0+j] = 1 for i, ixc0 := 0, 0; ixc0 <= jxc0; i++ { sum := 0. for k, kxcj := 0, j; k < i; k++ { sum += u.ele[kxcj] * l.ele[ixc0+k] kxcj += a.stride } u.ele[ixc0+j] = a.ele[ixc0+j] - sum ixc0 += a.stride } for ixc0 := jxc0; ixc0 < len(a.ele); ixc0 += a.stride { sum := 0. for k, kxcj := 0, j; k < j; k++ { sum += u.ele[kxcj] * l.ele[ixc0+k] kxcj += a.stride } l.ele[ixc0+j] = (a.ele[ixc0+j] - sum) / u.ele[jxc0+j] } jxc0 += a.stride } return
}
func main() {
showLU(&matrix{3, []float64{ 1, 3, 5, 2, 4, 7, 1, 1, 0}}) showLU(&matrix{4, []float64{ 11, 9, 24, 2, 1, 5, 2, 6, 3, 17, 18, 1, 2, 5, 7, 1}})
}
func showLU(a *matrix) {
a.print("\na") l, u, p := a.lu() l.print("l") u.print("u") p.print("p")
}</lang> Output is same as from 2D solution.
Library gonum/matrix
<lang go>package main
import (
"fmt"
"github.com/gonum/matrix/mat64"
)
func main() {
showLU(mat64.NewDense(3, 3, []float64{ 1, 3, 5, 2, 4, 7, 1, 1, 0, })) fmt.Println() showLU(mat64.NewDense(4, 4, []float64{ 11, 9, 24, 2, 1, 5, 2, 6, 3, 17, 18, 1, 2, 5, 7, 1, }))
}
func showLU(a *mat64.Dense) {
fmt.Printf("a: %v\n\n", mat64.Formatted(a, mat64.Prefix(" "))) var lu mat64.LU lu.Factorize(a) var l, u mat64.TriDense l.LFrom(&lu) u.UFrom(&lu) fmt.Printf("l: %.5f\n\n", mat64.Formatted(&l, mat64.Prefix(" "))) fmt.Printf("u: %.5f\n\n", mat64.Formatted(&u, mat64.Prefix(" "))) fmt.Println("p:", lu.Pivot(nil))
}</lang>
- Output:
Pivot format is a little different here. (But library solutions don't really meet task requirements anyway.)
a: ⎡1 3 5⎤ ⎢2 4 7⎥ ⎣1 1 0⎦ l: ⎡ 1.00000 0.00000 0.00000⎤ ⎢ 0.50000 1.00000 0.00000⎥ ⎣ 0.50000 -1.00000 1.00000⎦ u: ⎡ 2.00000 4.00000 7.00000⎤ ⎢ 0.00000 1.00000 1.50000⎥ ⎣ 0.00000 0.00000 -2.00000⎦ p: [1 0 2] a: ⎡11 9 24 2⎤ ⎢ 1 5 2 6⎥ ⎢ 3 17 18 1⎥ ⎣ 2 5 7 1⎦ l: ⎡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⎦ u: ⎡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⎦ p: [0 2 1 3]
Library go.matrix
<lang go>package main
import (
"fmt"
mat "github.com/skelterjohn/go.matrix"
)
func main() {
showLU(mat.MakeDenseMatrixStacked([][]float64{ {1, 3, 5}, {2, 4, 7}, {1, 1, 0}})) showLU(mat.MakeDenseMatrixStacked([][]float64{ {11, 9, 24, 2}, {1, 5, 2, 6}, {3, 17, 18, 1}, {2, 5, 7, 1}}))
}
func showLU(a *mat.DenseMatrix) {
fmt.Printf("\na:\n%v\n", a) l, u, p := a.LU() fmt.Printf("l:\n%v\n", l) fmt.Printf("u:\n%v\n", u) fmt.Printf("p:\n%v\n", p)
}</lang>
- Output:
a: {1, 3, 5, 2, 4, 7, 1, 1, 0} l: { 1, 0, 0, 0.5, 1, 0, 0.5, -1, 1} u: { 2, 4, 7, 0, 1, 1.5, 0, 0, -2} p: {0, 1, 0, 1, 0, 0, 0, 0, 1} a: {11, 9, 24, 2, 1, 5, 2, 6, 3, 17, 18, 1, 2, 5, 7, 1} l: { 1, 0, 0, 0, 0.272727, 1, 0, 0, 0.090909, 0.2875, 1, 0, 0.181818, 0.23125, 0.003597, 1} u: { 11, 9, 24, 2, 0, 14.545455, 11.454545, 0.454545, 0, 0, -3.475, 5.6875, 0, 0, 0, 0.510791} p: {1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1}
Idris
works with Idris 0.10
Uses The Method Of Partial Pivoting
Solution: <lang Idris> module Main
import Data.Vect
Matrix : Nat -> Nat -> Type -> Type Matrix m n t = Vect m (Vect n t)
-- Creates list from 0 to n (not including n) upTo : (m : Nat) -> Vect m (Fin m) upTo Z = [] upTo (S n) = 0 :: (map FS (upTo n))
-- Creates list from 0 to n-1 (not including n-1) upToM1 : (m : Nat) -> (sz ** Vect sz (Fin m)) upToM1 m = case (upTo m) of
(y::ys) => (_ ** init(y::ys)) [] => (_ ** [])
-- Creates list from i to n (not including n) fromUpTo : {n : Nat} -> Fin n -> (sz ** Vect sz (Fin n)) fromUpTo {n} m = filter (>= m) (upTo n)
-- Creates list from i+1 to n (not including n) fromUpTo1 : {n : Nat} -> Fin n -> (sz ** Vect sz (Fin n)) fromUpTo1 {n} m with (fromUpTo m)
| (_ ** xs) = case xs of (y::ys) => (_ ** ys) [] => (_ ** [])
-- Create Zero Matrix of size m by n
zeros : (m : Nat) -> (n : Nat) -> Matrix m n Double
zeros m n = replicate m (replicate n 0.0)
replaceAtM : (Fin m, Fin n) -> t -> Matrix m n t -> Matrix m n t replaceAtM (i, j) e a = replaceAt i (replaceAt j e (index i a)) a
-- Create Identity Matrix of size m by m eye : (m : Nat) -> Matrix m m Double eye m = map create1Vec (upTo m)
where set1 : Vect m Double -> Fin m -> Vect m Double set1 a n = replaceAt n 1.0 a
create1Vec : Fin m -> Vect m Double create1Vec n = set1 (replicate m 0.0) n
indexM : (Fin m, Fin n) -> Matrix m n t -> t
indexM (i, j) a = index j (index i a)
-- Obtain index for the row containing the
-- largest absolute value for the given column
colAbsMaxIndex : Fin m -> Fin m -> Matrix m m Double -> Fin m
colAbsMaxIndex startRow col a {m} with (fromUpTo startRow)
| (_ ** xs) = snd $ foldl (\(absMax, idx), curIdx => let curAbsVal = abs(indexM (curIdx, col) a) in if (curAbsVal > absMax) then (curAbsVal, curIdx) else (absMax, idx) ) (0.0, startRow) xs
-- Swaps two rows in a given matrix
swapRows : Fin m -> Fin m -> Matrix m n t -> Matrix m n t
swapRows r1 r2 a = replaceAt r2 tempRow $ replaceAt r1 (index r2 a) a
where tempRow = index r1 a
-- Swaps two individual values in a matrix
swapValues : (Fin m, Fin m) -> (Fin m, Fin m) -> Matrix m m Double -> Matrix m m Double
swapValues (i1, j1) (i2, j2) m = replaceAtM (i2, j2) v1 $ replaceAtM (i1, j1) v2 m
where v1 = indexM (i1, j1) m v2 = indexM (i2, j2) m
-- Perform row Swap on Lower Triangular Matrix lSwapRow : Fin m -> Fin m -> Matrix m m Double -> Matrix m m Double lSwapRow row1 row2 l {m} with (filter (< row1) (upTo m))
| (_ ** xs) = foldl (\l',col => swapValues (row1, col) (row2, col) l') l xs
rowSwap : Fin m -> (Matrix m m Double, Matrix m m Double, Matrix m m Double) ->
(Matrix m m Double, Matrix m m Double, Matrix m m Double)
rowSwap col (l,u,p) = (lSwapRow col row l, swapRows col row u, swapRows col row p)
where row = colAbsMaxIndex col col u
calc : (Fin m) -> (Fin m) -> (Matrix m m Double, Matrix m m Double) ->
(Matrix m m Double, Matrix m m Double)
calc i j (l, u) {m} = (l', u')
where l' : Matrix m m Double l' = replaceAtM (j, i) ((indexM (j, i) u) / indexM (i, i) u) l u : (Fin m) -> (Matrix m m Double) -> (Matrix m m Double) u k u = replaceAtM (j, k) ((indexM (j, k) u) - ((indexM (j, i) l') * (indexM (i, k) u))) u u' : (Matrix m m Double) u' with (fromUpTo i) | (_ ** xs) = foldl (\curU, idx => u idx curU) u xs
-- Perform a single iteration of the algorithm for the given column
iteration : Fin m -> (Matrix m m Double, Matrix m m Double, Matrix m m Double) ->
(Matrix m m Double, Matrix m m Double, Matrix m m Double)
iteration i lup {m} = iterate' (rowSwap i lup)
where modify : (Matrix m m Double, Matrix m m Double) -> (Matrix m m Double, Matrix m m Double) modify lu with (fromUpTo1 i) | (_ ** xs) = foldl (\lu',j => calc i j lu') lu xs iterate' : (Matrix m m Double, Matrix m m Double, Matrix m m Double) -> (Matrix m m Double, Matrix m m Double, Matrix m m Double) iterate' (l, u, p) with (modify (l, u)) | (l', u') = (l', u', p)
-- Generate L, U, P matricies from a given square matrix.
-- Where L * U = A, and P is the permutation matrix
luDecompose : Matrix m m Double -> (Matrix m m Double, Matrix m m Double, Matrix m m Double)
luDecompose a {m} with (upToM1 m)
| (_ ** xs) = foldl (\lup,idx => iteration idx lup) (eye m,a,eye m) xs
ex1 : (Matrix 3 3 Double, Matrix 3 3 Double, Matrix 3 3 Double) ex1 = luDecompose [[1, 3, 5], [2, 4, 7], [1, 1, 0]]
ex2 : (Matrix 4 4 Double, Matrix 4 4 Double, Matrix 4 4 Double) ex2 = luDecompose [[11, 9, 24, 2], [1, 5, 2, 6], [3, 17, 18, 1], [2, 5, 7, 1]]
printEx : (Matrix n n Double, Matrix n n Double, Matrix n n Double) -> IO () printEx (l, u, p) = do
putStr "l:" print l putStrLn "\n" putStr "u:" print u putStrLn "\n"
putStr "p:" print p putStrLn "\n"
main : IO() main = do
putStrLn "Solution 1:" printEx ex1 putStrLn "Solution 2:" printEx ex2
</lang>
- Output:
Solution 1: l:[[1, 0, 0], [0.5, 1, 0], [0.5, -1, 1]] u:[[2, 4, 7], [0, 1, 1.5], [0, 0, -2]] p:[[0, 1, 0], [1, 0, 0], [0, 0, 1]] Solution 2: l:[[1, 0, 0, 0], [0.2727272727272727, 1, 0, 0], [0.09090909090909091, 0.2875, 1, 0], [0.1818181818181818, 0.23125, 0.003597122302158069, 1]] u:[[11, 9, 24, 2], [0, 14.54545454545455, 11.45454545454546, 0.4545454545454546], [0, 0, -3.475, 5.6875], [0, 0, 0, 0.510791366906476]] p:[[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]]
J
Taken with slight modification from [1].
Solution: <lang j>mp=: +/ .*
LU=: 3 : 0
'm n'=. $ A=. y if. 1=m do. p ; (=1) ; p{"1 A [ p=. C. (n-1);~.0,(0~:,A)i.1 else. m2=. >.m%2 'p1 L1 U1'=. LU m2{.A D=. (/:p1) {"1 m2}.A F=. m2 {."1 D E=. m2 {."1 U1 FE1=. F mp %. E G=. m2}."1 D - FE1 mp U1 'p2 L2 U2'=. LU G p3=. (i.m2),m2+p2 H=. (/:p3) {"1 U1 (p1{p3) ; (L1,FE1,.L2) ; H,(-n){."1 U2 end.
)
permtomat=: 1 {.~"0 -@>:@:/: LUdecompose=: (permtomat&.>@{. , }.)@:LU</lang>
Example use: <lang j> A=:3 3$1 3 5 2 4 7 1 1 0
LUdecompose A
┌─────┬─────┬───────┐ │1 0 0│1 0 0│1 3 5│ │0 1 0│2 1 0│0 _2 _3│ │0 0 1│1 1 1│0 0 _2│ └─────┴─────┴───────┘
mp/> LUdecompose A
1 3 5 2 4 7 1 1 0
A=:4 4$11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 LUdecompose A
┌───────┬─────────────────────────────┬─────────────────────────────┐ │1 0 0 0│ 1 0 0 0│11 9 24 2│ │0 1 0 0│0.0909091 1 0 0│ 0 4.18182 _0.181818 5.81818│ │0 0 1 0│ 0.272727 3.47826 1 0│ 0 0 12.087 _19.7826│ │0 0 0 1│ 0.181818 0.804348 0.230216 1│ 0 0 0 0.510791│ └───────┴─────────────────────────────┴─────────────────────────────┘
mp/> LUdecompose A
11 9 24 2
1 5 2 6 3 17 18 1 2 5 7 1</lang>
Java
Translation of Common Lisp via D
<lang java>import static java.util.Arrays.stream; import java.util.Locale; import static java.util.stream.IntStream.range;
public class Test {
static double dotProduct(double[] a, double[] b) { return range(0, a.length).mapToDouble(i -> a[i] * b[i]).sum(); }
static double[][] matrixMul(double[][] A, double[][] B) { double[][] result = new double[A.length][B[0].length]; double[] aux = new double[B.length];
for (int j = 0; j < B[0].length; j++) {
for (int k = 0; k < B.length; k++) aux[k] = B[k][j];
for (int i = 0; i < A.length; i++) result[i][j] = dotProduct(A[i], aux); } return result; }
static double[][] pivotize(double[][] m) { int n = m.length; double[][] id = range(0, n).mapToObj(j -> range(0, n) .mapToDouble(i -> i == j ? 1 : 0).toArray()) .toArray(double[][]::new);
for (int i = 0; i < n; i++) { double maxm = m[i][i]; int row = i; for (int j = i; j < n; j++) if (m[j][i] > maxm) { maxm = m[j][i]; row = j; }
if (i != row) { double[] tmp = id[i]; id[i] = id[row]; id[row] = tmp; } } return id; }
static double[][][] lu(double[][] A) { int n = A.length; double[][] L = new double[n][n]; double[][] U = new double[n][n]; double[][] P = pivotize(A); double[][] A2 = matrixMul(P, A);
for (int j = 0; j < n; j++) { L[j][j] = 1; for (int i = 0; i < j + 1; i++) { double s1 = 0; for (int k = 0; k < i; k++) s1 += U[k][j] * L[i][k]; U[i][j] = A2[i][j] - s1; } for (int i = j; i < n; i++) { double s2 = 0; for (int k = 0; k < j; k++) s2 += U[k][j] * L[i][k]; L[i][j] = (A2[i][j] - s2) / U[j][j]; } } return new double[][][]{L, U, P}; }
static void print(double[][] m) { stream(m).forEach(a -> { stream(a).forEach(n -> System.out.printf(Locale.US, "%5.1f ", n)); System.out.println(); }); System.out.println(); }
public static void main(String[] args) { double[][] a = {{1.0, 3, 5}, {2.0, 4, 7}, {1.0, 1, 0}};
double[][] b = {{11.0, 9, 24, 2}, {1.0, 5, 2, 6}, {3.0, 17, 18, 1}, {2.0, 5, 7, 1}};
for (double[][] m : lu(a)) print(m);
System.out.println();
for (double[][] m : lu(b)) print(m); }
}</lang>
1.0 0.0 0.0 0.5 1.0 0.0 0.5 -1.0 1.0 2.0 4.0 7.0 0.0 1.0 1.5 0.0 0.0 -2.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.3 1.0 0.0 0.0 0.1 0.3 1.0 0.0 0.2 0.2 0.0 1.0 11.0 9.0 24.0 2.0 0.0 14.5 11.5 0.5 0.0 0.0 -3.5 5.7 0.0 0.0 0.0 0.5 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0
jq
jq currently does not have builtin support for matrices and therefore some infrastructure is needed to make the following self-contained. Matrices here are represented as arrays of arrays in the usual way.
Infrastructure <lang jq># Create an m x n matrix def matrix(m; n; init):
if m == 0 then [] elif m == 1 then [range(0;n)] | map(init) elif m > 0 then matrix(1;n;init) as $row | [range(0;m)] | map( $row ) else error("matrix\(m);_;_) invalid") end ;
def I(n): matrix(n;n;0) as $m
| reduce range(0;n) as $i ($m; . | setpath( [$i,$i]; 1));
def dot_product(a; b):
reduce range(0;a|length) as $i (0; . + (a[$i] * b[$i]) );
- transpose/0 expects its input to be a rectangular matrix
def transpose:
if (.[0] | length) == 0 then [] else [map(.[0])] + (map(.[1:]) | transpose) end ;
- A and B should both be numeric matrices, A being m by n, and B being n by p.
def multiply(A; B):
(B[0]|length) as $p | (B|transpose) as $BT | reduce range(0; A|length) as $i ([]; reduce range(0; $p) as $j (.; .[$i][$j] = dot_product( A[$i]; $BT[$j] ) ));
def swap_rows(i;j):
if i == j then . else .[i] as $i | .[i] = .[j] | .[j] = $i end ;
- Print a matrix neatly, each cell occupying n spaces, but without truncation
def neatly(n):
def right: tostring | ( " " * (n-length) + .); . as $in | length as $length | reduce range (0;$length) as $i (""; . + reduce range(0;$length) as $j (""; "\(.) \($in[$i][$j] | right )" ) + "\n" ) ;</lang>
LU decomposition <lang jq># Create the pivot matrix for the input matrix.
- Use "range(0;$n) as $i" to handle ill-conditioned cases.
def pivotize:
def abs: if .<0 then -. else . end; length as $n | . as $m | reduce range(0;$n) as $j (I($n); # state: [row; max] (reduce range(0; $n) as $i ([$j, $m[$j][$j]|abs ]; ($m[$i][$j]|abs) as $a | if $a > .[1] then [ $i, $a ] else . end) | .[0]) as $row | swap_rows( $j; $row) ) ;
- Decompose the input nxn matrix A by PA=LU and return [L, U, P].
def lup:
def div(i;j): if j == 0 then if i==0 then 0 else error("\(i)/0") end else i/j end; . as $A | length as $n | I($n) as $L # matrix($n; $n; 0.0) as $L | matrix($n; $n; 0.0) as $U | ($A|pivotize) as $P | multiply($P;$A) as $A2 # state: [L, U] | reduce range(0; $n) as $i ( [$L, $U]; reduce range(0; $n) as $j (.; .[0] as $L | .[1] as $U | if ($j >= $i) then (reduce range(0;$i) as $k (0; . + ($U[$k][$j] * $L[$i][$k] ))) as $s1 | [$L, ($U| setpath([$i,$j]; ($A2[$i][$j] - $s1))) ] else (reduce range(0;$j) as $k (0; . + ($U[$k][$j] * $L[$i][$k]))) as $s2 | [ ($L | setpath([$i,$j]; div(($A2[$i][$j] - $s2) ; $U[$j][$j] ))), $U ] end )) | . + [ $P ]
</lang> Example 1: <lang jq>def a: [[1, 3, 5], [2, 4, 7], [1, 1, 0]]; a | lup[] | neatly(4) </lang>
- Output:
<lang sh> $ /usr/local/bin/jq -M -n -r -f LU.jq
1 0 0 0.5 1 0 0.5 -1 1
2 4 7 0 1 1.5 0 0 -2
0 1 0 1 0 0 0 0 1
</lang> Example 2: <lang jq>def b: [[11,9,24,2],[1,5,2,6],[3,17,18,1],[2,5,7,1]]; b | lup[] | neatly(21)</lang>
- Output:
<lang sh>$ /usr/local/bin/jq -M -n -r -f LU.jq
1 0 0 0 0.2727272727272727 1 0 0 0.09090909090909091 0.2875 1 0 0.18181818181818182 0.23124999999999996 0.0035971223021580693 1
11 9 24 2 0 14.545454545454547 11.454545454545455 0.4545454545454546 0 0 -3.4749999999999996 5.6875 0 0 0 0.510791366906476
1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1</lang>
Example 3: <lang jq>
- A|lup|verify(A) should be true
def verify(A):
.[0] as $L | .[1] as $U | .[2] as $P | multiply($P; A) == multiply($L; $U);
def A:
[[1, 1, 1, 1], [1, 1, -1, -1], [1, -1, 0, 0], [0, 0, 1, -1]];
A|lup|verify(A)</lang>
- Output:
true
Julia
Julia has the predefined functions `lu`, `lufact` and `lufact!` in the standard library to compute the lu decomposition of a matrix.
- Output:
julia> lu([1 3 5 ; 2 4 7 ; 1 1 0]) ( 3x3 Array{Float64,2}: 1.0 0.0 0.0 0.5 1.0 0.0 0.5 -1.0 1.0, 3x3 Array{Float64,2}: 2.0 4.0 7.0 0.0 1.0 1.5 0.0 0.0 -2.0, [2,1,3])
Maple
<lang Maple> A:=<<1.0|3.0|5.0>,<2.0|4.0|7.0>,<1.0|1.0|0.0>>:
LinearAlgebra:-LUDecomposition(A); </lang>
- Output:
[0 1 0] [ 1.0 0. 0.] [2. 4. 7.] [ ] [ ] [ ] [1 0 0], [0.500000000000000 1.0 0.], [0. 1. 1.50000000000000] [ ] [ ] [ ] [0 0 1] [0.500000000000000 -1. 1.0] [0. 0. -2.]
<lang Maple> A:=<<11.0|9.0|24.0|2.0>,<1.0|5.0|2.0|6.0>,
<3.0|17.0|18.0|1.0>,<2.0|5.0|7.0|1.0>>:
with(LinearAlgebra):
LUDecomposition(A); </lang>
- Output:
[1 0 0 0] [ ] [0 0 1 0] [ ], [0 1 0 0] [ ] [0 0 0 1] [ 1.0 0. 0. 0.] [ ] [ 0.272727272727273 1.0 0. 0.] [ ], [0.0909090909090909 0.287500000000000 1.0 0.] [ ] [ 0.181818181818182 0.231250000000000 0.00359712230215807 1.0] [11. 9. 24. 2.] [ ] [ 0. 14.5454545454545 11.4545454545455 0.454545454545455] [ ] [ 0. 0. -3.47500000000000 5.68750000000000] [ ] [ 0. 0. 0. 0.510791366906476]
Mathematica
<lang Mathematica>(*Ex1*)a = {{1, 3, 5}, {2, 4, 7}, {1, 1, 0}}; {lu, p, c} = LUDecomposition[a]; l = LowerTriangularize[lu, -1] + IdentityMatrix[Length[p]]; u = UpperTriangularize[lu]; P = Part[IdentityMatrix[Length[p]], p] ; MatrixForm /@ {P.a , P, l, u, l.u}
(*Ex2*)a = {{11, 9, 24, 2}, {1, 5, 2, 6}, {3, 17, 18, 1}, {2, 5, 7, 1}}; {lu, p, c} = LUDecomposition[a]; l = LowerTriangularize[lu, -1] + IdentityMatrix[Length[p]]; u = UpperTriangularize[lu]; P = Part[IdentityMatrix[Length[p]], p] ; MatrixForm /@ {P.a , P, l, u, l.u} </lang>
- Output:
MATLAB / Octave
LU decomposition is part of language
<lang Matlab> A = [
1 3 5 2 4 7 1 1 0];
[L,U,P] = lu(A)</lang>
- Output:
L = 1.00000 0.00000 0.00000 0.50000 1.00000 0.00000 0.50000 -1.00000 1.00000 U = 2.00000 4.00000 7.00000 0.00000 1.00000 1.50000 0.00000 0.00000 -2.00000 P = 0 1 0 1 0 0 0 0 1
2nd example: <lang Matlab> A = [
11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 ];
[L,U,P] = lu(A)</lang>
- Output:
L = 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 U = 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 P = 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1
Creating a MATLAB function
<lang Matlab> function [ P, L, U ] = LUdecomposition(A)
% Ensures A is n by n sz = size(A); if sz(1)~=sz(2)
fprintf('A is not n by n\n'); clear x; return;
end
n = sz(1); L = eye(n); P = eye(n); U = A;
for i=1:sz(1)
% Row reducing if U(i,i)==0 maximum = max(abs(U(i:end,1))); for k=1:n if maximum == abs(U(k,i)) temp = U(1,:); U(1,:) = U(k,:); U(k,:) = temp; temp = P(:,1); P(1,:) = P(k,:); P(k,:) = temp; end end end if U(i,i)~=1 temp = eye(n); temp(i,i)=U(i,i); L = L * temp; U(i,:) = U(i,:)/U(i,i); %Ensures the pivots are 1. end if i~=sz(1) for j=i+1:length(U) temp = eye(n); temp(j,i) = U(j,i); L = L * temp; U(j,:) = U(j,:)-U(j,i)*U(i,:); end end
end P = P'; end
</lang>
Maxima
<lang maxima>/* LU decomposition is built-in */
a: hilbert_matrix(4)$
/* LU in "packed" form */
lup: lu_factor(a); /* [matrix([1, 1/2, 1/3, 1/4 ],
[1/2, 1/12, 1/12, 3/40 ], [1/3, 1, 1/180, 1/120 ], [1/4, 9/10, 3/2, 1/2800]), [1, 2, 3, 4], generalring] */
/* extract actual factors */
get_lu_factors(lup); /* [matrix([1, 0, 0, 0],
[0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]), matrix([1, 0, 0, 0], [1/2, 1, 0, 0], [1/3, 1, 1, 0], [1/4, 9/10, 3/2, 1]), matrix([1, 1/2, 1/3, 1/4 ], [0, 1/12, 1/12, 3/40 ], [0, 0, 1/180, 1/120 ], [0, 0, 0, 1/2800]) ] */
/* solve for a given right-hand side */
lu_backsub(lup, transpose([1, 1, -1, -1])); /* matrix([-204], [2100], [-4740], [2940]) */</lang>
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:
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]
Perl 6
Translation of Ruby.
<lang perl6>for ( [1, 3, 5], # Test Matrices
[2, 4, 7], [1, 1, 0] ), ( [11, 9, 24, 2], [ 1, 5, 2, 6], [ 3, 17, 18, 1], [ 2, 5, 7, 1] ) -> @test { say-it 'A Matrix', @test; say-it( $_[0], @($_[1]) ) for 'P Matrix', 'Aʼ Matrix', 'L Matrix', 'U Matrix' Z, lu @test;
}
sub lu (@a) {
die unless @a.&is-square; my $n = +@a; my @P = pivotize @a; my @Aʼ = mmult @P, @a; my @L = matrix-ident $n; my @U = matrix-zero $n; for ^$n -> $i { for ^$n -> $j { if $j >= $i { @U[$i][$j] = @Aʼ[$i][$j] - [+] map { @U[$_][$j] * @L[$i][$_] }, ^$i } else { @L[$i][$j] = (@Aʼ[$i][$j] - [+] map { @U[$_][$j] * @L[$i][$_] }, ^$j) / @U[$j][$j]; } }
} return @P, @Aʼ, @L, @U;
}
sub pivotize (@m) {
my $size = +@m; my @id = matrix-ident $size; for ^$size -> $i { my $max = @m[$i][$i]; my $row = $i; for $i ..^ $size -> $j { if @m[$j][$i] > $max { $max = @m[$j][$i]; $row = $j; } } if $row != $i { @id[$row, $i] = @id[$i, $row] } } @id
}
sub is-square (@m) { so @m == all @m[*] }
sub matrix-zero ($n, $m = $n) { map { [ flat 0 xx $n ] }, ^$m }
sub matrix-ident ($n) { map { [ flat 0 xx $_, 1, 0 xx $n - 1 - $_ ] }, ^$n }
sub mmult(@a,@b) {
my @p; for ^@a X ^@b[0] -> ($r, $c) { @p[$r][$c] += @a[$r][$_] * @b[$_][$c] for ^@b; } @p
}
sub rat-int ($num) {
return $num unless $num ~~ Rat; return $num.narrow if $num.narrow.WHAT ~~ Int; $num.nude.join: '/';
}
sub say-it ($message, @array) {
say "\n$message"; $_».&rat-int.fmt("%7s").say for @array;
}</lang>
- Output:
A Matrix 1 3 5 2 4 7 1 1 0 P Matrix 0 1 0 1 0 0 0 0 1 Aʼ Matrix 2 4 7 1 3 5 1 1 0 L Matrix 1 0 0 1/2 1 0 1/2 -1 1 U Matrix 2 4 7 0 1 3/2 0 0 -2 A Matrix 11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 P Matrix 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 Aʼ Matrix 11 9 24 2 3 17 18 1 1 5 2 6 2 5 7 1 L Matrix 1 0 0 0 3/11 1 0 0 1/11 23/80 1 0 2/11 37/160 1/278 1 U Matrix 11 9 24 2 0 160/11 126/11 5/11 0 0 -139/40 91/16 0 0 0 71/139
PL/I
<lang PL/I>(subscriptrange, fofl, size): /* 2 Nov. 2013 */ LU_Decomposition: procedure options (main);
declare a1(3,3) float (18) initial ( 1, 3, 5, 2, 4, 7, 1, 1, 0); declare a2(4,4) float (18) initial (11, 9, 24, 2, 1, 5, 2, 6, 3, 17, 18, 1, 2, 5, 7, 1); call check(a1); call check(a2);
/* In-situ decomposition */
LU: procedure(a, p);
declare a(*,*) float (18); declare p(*) fixed binary; declare (maximum, rtemp) float (18); declare (n, i, j, k, ii, temp) fixed binary;
n = hbound(a,1); do i = 1 to n; p(i) = i; end;
do k = 1 to n-1;
maximum = 0; ii = k; do i = k to n; if maximum < abs(a(p(i),k)) then do; maximum = abs(a(p(i),k)); ii = i; end; end; if ii ^= k then do; temp = p(k); p(k) = p(ii); p(ii) = temp; end;
do i = k+1 to n; a(p(i),k) = a(p(i),k) / a(p(k),k); end;
do j = k+1 to n; do i = k+1 to n; a(p(i),j) = a(p(i),j) - a(p(i),k) * a(p(k),j); end; end;
end;
end LU;
CHECK: procedure(a);
declare a(*,*) float (18) nonassignable;
declare aa(hbound(a,1), hbound(a,2)) float (18); declare L(hbound(a,1), hbound(a,2)) float (18); declare U(hbound(a,1), hbound(a,2)) float (18); declare (p(hbound(a,1), hbound(a,2)), ipiv(hbound(a,1)) ) fixed binary; declare pp(hbound(a,1), hbound(a,2)) fixed binary; declare (i, j, n, temp(hbound(a,1))) fixed binary;
n = hbound(a,1); aa = A; /* work with a copy */ P = 0; L = 0; U = 0; do i = 1 to n; p(i,i) = 1; L(i,i) = 1; /* convert permutation vector to a matrix */ end;
call LU(aa, ipiv);
do i = 1 to n; do j = 1 to i-1; L(i,j) = aa(ipiv(i),j); end; do j = i to n; U(i,j) = aa(ipiv(i),j); end; end;
pp = p; do i = 1 to n; p(ipiv(i), *) = pp(i,*); end;
put skip list ('A'); put edit (A) (skip, (n) f(10,5));
put skip list ('P'); put edit (P) (skip, (n) f(11));
put skip list ('L'); put edit (L) (skip, (n) f(10,5));
put skip list ('U'); put edit (U) (skip, (n) f(10,5));
end CHECK;
end LU_Decomposition; </lang> Derived from Fortran version above. Results:
A 1.00000 3.00000 5.00000 2.00000 4.00000 7.00000 1.00000 1.00000 0.00000 P 0 1 0 1 0 0 0 0 1 L 1.00000 0.00000 0.00000 0.50000 1.00000 0.00000 0.50000 -1.00000 1.00000 U 2.00000 4.00000 7.00000 0.00000 1.00000 1.50000 0.00000 0.00000 -2.00000 A 11.00000 9.00000 24.00000 2.00000 1.00000 5.00000 2.00000 6.00000 3.00000 17.00000 18.00000 1.00000 2.00000 5.00000 7.00000 1.00000 P 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 L 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 U 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
Python
<lang python>from pprint import pprint
def matrixMul(A, B):
TB = zip(*B) return [[sum(ea*eb for ea,eb in zip(a,b)) for b in TB] for a in A]
def pivotize(m):
"""Creates the pivoting matrix for m.""" n = len(m) ID = [[float(i == j) for i in xrange(n)] for j in xrange(n)] for j in xrange(n): row = max(xrange(j, n), key=lambda i: abs(m[i][j])) if j != row: ID[j], ID[row] = ID[row], ID[j] return ID
def lu(A):
"""Decomposes a nxn matrix A by PA=LU and returns L, U and P.""" n = len(A) L = [[0.0] * n for i in xrange(n)] U = [[0.0] * n for i in xrange(n)] P = pivotize(A) A2 = matrixMul(P, A) for j in xrange(n): L[j][j] = 1.0 for i in xrange(j+1): s1 = sum(U[k][j] * L[i][k] for k in xrange(i)) U[i][j] = A2[i][j] - s1 for i in xrange(j, n): s2 = sum(U[k][j] * L[i][k] for k in xrange(j)) L[i][j] = (A2[i][j] - s2) / U[j][j] return (L, U, P)
a = [[1, 3, 5], [2, 4, 7], [1, 1, 0]] for part in lu(a):
pprint(part, width=19) print
print b = [[11,9,24,2],[1,5,2,6],[3,17,18,1],[2,5,7,1]] for part in lu(b):
pprint(part) print</lang>
- Output:
[[1.0, 0.0, 0.0], [0.5, 1.0, 0.0], [0.5, -1.0, 1.0]] [[2.0, 4.0, 7.0], [0.0, 1.0, 1.5], [0.0, 0.0, -2.0]] [[0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]] [[1.0, 0.0, 0.0, 0.0], [0.27272727272727271, 1.0, 0.0, 0.0], [0.090909090909090912, 0.28749999999999998, 1.0, 0.0], [0.18181818181818182, 0.23124999999999996, 0.0035971223021580693, 1.0]] [[11.0, 9.0, 24.0, 2.0], [0.0, 14.545454545454547, 11.454545454545455, 0.45454545454545459], [0.0, 0.0, -3.4749999999999996, 5.6875], [0.0, 0.0, 0.0, 0.51079136690647597]] [[1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
R
- Output:
> A <- c(1, 2, 1, 3, 4, 1, 5, 7, 0) > dim(A) <- c(3, 3) > library(Matrix) > expand(lu(A)) $L 3 x 3 Matrix of class "dtrMatrix" (unitriangular) [,1] [,2] [,3] [1,] 1.0 . . [2,] 0.5 1.0 . [3,] 0.5 -1.0 1.0 $U 3 x 3 Matrix of class "dtrMatrix" [,1] [,2] [,3] [1,] 2.0 4.0 7.0 [2,] . 1.0 1.5 [3,] . . -2.0 $P 3 x 3 sparse Matrix of class "pMatrix" [1,] . | . [2,] | . . [3,] . . |
Racket
<lang racket>
- lang racket
(require math) (define A (matrix
[[1 3 5] [2 4 7] [1 1 0]]))
(matrix-lu A)
- result
- (mutable-array #[#[1 0 0]
- #[2 1 0]
- #[1 1 1]])
- (mutable-array #[#[1 3 5]
- #[0 -2 -3]
- #[0 0 -2]])
</lang>
REXX
<lang rexx>/*REXX program creates a matrix from console input, performs/shows LU decomposition.*/
- =0; P.=0; PA.=0; L.=0; U.=0 /*initialize some variables to zero. */
parse arg x /*obtain matrix elements from the C.L. */ call makeMat /*make the A matrix from the numbers.*/ call showMat 'A', N /*display the A matrix. */ call manPmat /*manufacture P (permutation). */ call showMat 'P', N /*display the P matrix. */ call multMat /*multiply the A and P matrices. */ call showMat 'PA', N /*display the PA matrix. */
do y=1 for N; call manUmat y /*manufacture U matrix, parts. */ call manLmat y /*manufacture L matrix, parts. */ end
call showMat 'L', N /*display the L matrix. */ call showMat 'U', N /*display the U matrix. */ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ er: say; say '***error!***'; say; say arg(1); say; exit 13 /*──────────────────────────────────────────────────────────────────────────────────────*/ makeMat: ?=words(x); do N=1 for ?; if N**2==? then leave; end /*N*/
if N**2\==? then call er 'not correct number of elements entered: ' ?
do r=1 for N /*build the "A" matrix from the input*/ do c=1 for N; #=#+1; _=word(x,#); A.r.c=_ if \datatype(_,'N') then call er "element isn't numeric: " _ end /*c*/ end /*r*/ return
/*──────────────────────────────────────────────────────────────────────────────────────*/ manLmat: parse arg ? /*manufacture L (lower) matrix.*/
do r=1 for N do c=1 for N; if r==c then do; L.r.c=1; iterate; end if c\==? | r==c | c>r then iterate _=PA.r.c do k=1 for c-1; _=_-U.k.c*L.r.k; end /*k*/ L.r.c=_/U.c.c end /*c*/ end /*r*/ return
/*──────────────────────────────────────────────────────────────────────────────────────*/ manPmat: c=N; do r=N by -1 for N /*manufacture P (permutation). */
P.r.c=1; c=c+1; if c>N then c=N%2; if c==N then c=1 end /*r*/ return
/*──────────────────────────────────────────────────────────────────────────────────────*/ manUmat: parse arg ? /*manufacture U (upper) matrix.*/
do r=1 for N; if r\==? then iterate do c=1 for N; if c<r then iterate _=PA.r.c do k=1 for r-1; _=_-U.k.c*L.r.k; end /*k*/ U.r.c=_/1 end /*c*/ end /*r*/ return
/*──────────────────────────────────────────────────────────────────────────────────────*/ multMat: do i=1 for N /*multiply matrix P & A ──► PA */
do j=1 for N do k=1 for N; pa.i.j=(pa.i.j + p.i.k * a.k.j) / 1 end /*k*/ end /*j*/ end /*i*/ return
/*──────────────────────────────────────────────────────────────────────────────────────*/ showMat: parse arg mat,rows,cols; w=0; cols=word(cols rows,1); say
do r=1 for rows do c=1 for cols; w=max(w, length( value( mat'.'r"."c ) ) ) end /*c*/ end /*r*/ say center(mat 'matrix',cols*(w+1)+7,"─") do r=1 for rows; _= do c=1 for cols; _=_ right(value(mat'.'r'.'c),w+1); end /*c*/ say _ end /*r*/ return</lang>
output when using the input of: 1 3 5 2 4 7 1 1 0
──A matrix─── 1 3 5 2 4 7 1 1 0 ──P matrix─── 0 1 0 1 0 0 0 0 1 ──PA matrix── 2 4 7 1 3 5 1 1 0 ─────L matrix────── 1 0 0 0.5 1 0 0.5 -1 1 ─────U matrix────── 2 4 7 0 1 1.5 0 0 -2
output when using the input of: 11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1
─────A matrix────── 11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 ───P matrix──── 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 ─────PA matrix───── 11 9 24 2 3 17 18 1 1 5 2 6 2 5 7 1 ───────────────────────────L matrix──────────────────────────── 1 0 0 0 0.272727273 1 0 0 0.0909090909 0.287500001 1 0 0.181818182 0.23125 0.00359712804 1 ───────────────────────U matrix──────────────────────── 11 9 24 2 0 14.5454545 11.4545455 0.45454545 0 0 -3.47500002 5.6875 0 0 0 0.510791339
Ruby
<lang ruby>require 'matrix'
class Matrix
def lu_decomposition p = get_pivot tmp = p * self u = Matrix.zero(row_size).to_a l = Matrix.identity(row_size).to_a (0 ... row_size).each do |i| (0 ... row_size).each do |j| if j >= i # upper u[i][j] = tmp[i,j] - (0 ... i).inject(0.0) {|sum, k| sum + u[k][j] * l[i][k]} else # lower l[i][j] = (tmp[i,j] - (0 ... j).inject(0.0) {|sum, k| sum + u[k][j] * l[i][k]}) / u[j][j] end end end [ Matrix[*l], Matrix[*u], p ] end def get_pivot raise ArgumentError, "must be square" unless square? id = Matrix.identity(row_size).to_a (0 ... row_size).each do |i| max = self[i,i] row = i (i ... row_size).each do |j| if self[j,i] > max max = self[j,i] row = j end end id[i], id[row] = id[row], id[i] end Matrix[*id] end def pretty_print(format, head=nil) puts head if head puts each_slice(column_size).map{|row| format*row_size % row} end
end
puts "Example 1:" a = Matrix[[1, 3, 5],
[2, 4, 7], [1, 1, 0]]
a.pretty_print(" %2d", "A") l, u, p = a.lu_decomposition l.pretty_print(" %8.5f", "L") u.pretty_print(" %8.5f", "U") p.pretty_print(" %d", "P")
puts "\nExample 2:" a = Matrix[[11, 9,24,2],
[ 1, 5, 2,6], [ 3,17,18,1], [ 2, 5, 7,1]]
a.pretty_print(" %2d", "A") l, u, p = a.lu_decomposition l.pretty_print(" %8.5f", "L") u.pretty_print(" %8.5f", "U") p.pretty_print(" %d", "P")</lang>
- Output:
Example 1: A 1 3 5 2 4 7 1 1 0 L 1.00000 0.00000 0.00000 0.50000 1.00000 0.00000 0.50000 -1.00000 1.00000 U 2.00000 4.00000 7.00000 0.00000 1.00000 1.50000 0.00000 0.00000 -2.00000 P 0 1 0 1 0 0 0 0 1 Example 2: A 11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 L 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 U 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 P 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1
Matrix has a lup_decomposition
built-in method.
<lang ruby>l, u, p = a.lup_decomposition
l.pretty_print(" %8.5f", "L")
u.pretty_print(" %8.5f", "U")
p.pretty_print(" %d", "P")</lang>
Output is the same.
Sidef
<lang ruby>func is_square(m) { m.all { .len == m.len } } func matrix_zero(n, m=n) { m.of { n.of(0) } } func matrix_ident(n) { n.of {|i| [(i-1).of(0)..., 1, (n-i).of(0)...] } }
func pivotize(m) {
var size = m.len var id = matrix_ident(size) for i in ^size { var max = m[i][i] var row = i for j in range(i, size-1) { if (m[j][i] > max) { max = m[j][i] row = j } } if (row != i) { id.swap(row, i) } } return id
}
func mmult(a, b) {
var p = [] for r,c in (^a ~X ^b[0]) { for i in ^b { p[r][c] := 0 += (a[r][i] * b[i][c]) } } return p
}
func lu(a) {
is_square(a) || die "Defined only for square matrices!"; var n = a.len var P = pivotize(a) var Aʼ = mmult(P, a) var L = matrix_ident(n) var U = matrix_zero(n) for i,j in (^n ~X ^n) { if (j >= i) { U[i][j] = (Aʼ[i][j] - (^i->map { U[_][j] * L[i][_] }.sum(0))) } else { L[i][j] = ((Aʼ[i][j] - (^j->map { U[_][j] * L[i][_] }.sum(0))) / U[j][j]) } } return [P, Aʼ, L, U]
}
func say_it(message, array) {
say "\n#{message}" array.each { |row| say row.map{"%7s" % .as_rat}.join(' ') }
}
var t = [[
%n(1 3 5), %n(2 4 7), %n(1 1 0),
],[
%n(11 9 24 2), %n( 1 5 2 6), %n( 3 17 18 1), %n( 2 5 7 1),
]]
t.each { |test|
say_it('A Matrix', test); for a in (['P Matrix', 'Aʼ Matrix', 'L Matrix', 'U Matrix'] ~Z lu(test)) { say_it(a[0], a[1]) }
}</lang>
A Matrix 1 3 5 2 4 7 1 1 0 P Matrix 0 1 0 1 0 0 0 0 1 Aʼ Matrix 2 4 7 1 3 5 1 1 0 L Matrix 1 0 0 1/2 1 0 1/2 -1 1 U Matrix 2 4 7 0 1 3/2 0 0 -2 A Matrix 11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 P Matrix 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 Aʼ Matrix 11 9 24 2 3 17 18 1 1 5 2 6 2 5 7 1 L Matrix 1 0 0 0 3/11 1 0 0 1/11 23/80 1 0 2/11 37/160 1/278 1 U Matrix 11 9 24 2 0 160/11 126/11 5/11 0 0 -139/40 91/16 0 0 0 71/139
Tcl
<lang tcl>package require Tcl 8.5 namespace eval matrix {
namespace path {::tcl::mathfunc ::tcl::mathop}
# Construct an identity matrix of the given size proc identity {order} {
set m [lrepeat $order [lrepeat $order 0]] for {set i 0} {$i < $order} {incr i} { lset m $i $i 1 } return $m
}
# Produce the pivot matrix for a given matrix proc pivotize {matrix} {
set n [llength $matrix] set p [identity $n] for {set j 0} {$j < $n} {incr j} { set max [lindex $matrix $j $j] set row $j for {set i $j} {$i < $n} {incr i} { if {[lindex $matrix $i $j] > $max} { set max [lindex $matrix $i $j] set row $i } } if {$j != $row} { # Row swap inlined; too trivial to have separate procedure set tmp [lindex $p $j] lset p $j [lindex $p $row] lset p $row $tmp } } return $p
}
# Decompose a square matrix A by PA=LU and return L, U and P proc luDecompose {A} {
set n [llength $A] set L [lrepeat $n [lrepeat $n 0]] set U $L set P [pivotize $A] set A [multiply $P $A]
for {set j 0} {$j < $n} {incr j} { lset L $j $j 1 for {set i 0} {$i <= $j} {incr i} { lset U $i $j [- [lindex $A $i $j] [SumMul $L $U $i $j $i]] } for {set i $j} {$i < $n} {incr i} { set sum [SumMul $L $U $i $j $j] lset L $i $j [/ [- [lindex $A $i $j] $sum] [lindex $U $j $j]] } }
return [list $L $U $P]
}
# Helper that makes inner loop nicer; multiplies column and row, # possibly partially... proc SumMul {A B i j kmax} {
set s 0.0 for {set k 0} {$k < $kmax} {incr k} { set s [+ $s [* [lindex $A $i $k] [lindex $B $k $j]]] } return $s
}
}</lang> Support code: <lang tcl># Code adapted from Matrix_multiplication and Matrix_transposition tasks namespace eval matrix {
# Get the size of a matrix; assumes that all rows are the same length, which # is a basic well-formed-ness condition... proc size {m} {
set rows [llength $m] set cols [llength [lindex $m 0]] return [list $rows $cols]
}
# Matrix multiplication implementation proc multiply {a b} {
lassign [size $a] a_rows a_cols lassign [size $b] b_rows b_cols if {$a_cols != $b_rows} { error "incompatible sizes: a($a_rows, $a_cols), b($b_rows, $b_cols)" } set temp [lrepeat $a_rows [lrepeat $b_cols 0]] for {set i 0} {$i < $a_rows} {incr i} { for {set j 0} {$j < $b_cols} {incr j} { lset temp $i $j [SumMul $a $b $i $j $a_cols] } } return $temp
}
# Pretty printer for matrices proc print {matrix {fmt "%g"}} {
set max [Widest $matrix $fmt] lassign [size $matrix] rows cols foreach row $matrix { foreach val $row width $max { puts -nonewline [format "%*s " $width [format $fmt $val]] } puts "" }
} proc Widest {m fmt} {
lassign [size $m] rows cols set max [lrepeat $cols 0] foreach row $m { for {set j 0} {$j < $cols} {incr j} { set s [format $fmt [lindex $row $j]] lset max $j [max [lindex $max $j] [string length $s]] } } return $max
}
}</lang> Demonstrating: <lang tcl># This does the decomposition and prints it out nicely proc demo {A} {
lassign $A L U P foreach v {A L U P} {
upvar 0 $v matrix puts "${v}:" matrix::print $matrix %.5g if {$v ne "P"} {puts "---------------------------------"}
}
} demo {{1 3 5} {2 4 7} {1 1 0}} puts "=================================" demo {{11 9 24 2} {1 5 2 6} {3 17 18 1} {2 5 7 1}}</lang>
- Output:
A: 1 3 5 2 4 7 1 1 0 --------------------------------- L: 1 0 0 0.5 1 0 0.5 -1 1 --------------------------------- U: 2 4 7 0 1 1.5 0 0 -2 --------------------------------- P: 0 1 0 1 0 0 0 0 1 ================================= A: 11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 --------------------------------- L: 1 0 0 0 0.27273 1 0 0 0.090909 0.2875 1 0 0.18182 0.23125 0.0035971 1 --------------------------------- U: 11 9 24 2 0 14.545 11.455 0.45455 0 0 -3.475 5.6875 0 0 0 0.51079 --------------------------------- P: 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1
zkl
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:=eye(n); U:=make_array(n,n,0.0); P:=pivotize(A); A=matMult(P,A);
foreach j in (n){ 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>
- Output:
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))
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).
- Output:
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