LU decomposition
From Rosetta Code
You are encouraged to solve this task according to the task description, using any language you may know.
Every square matrix A can be decomposed into a product of a lower triangular matrix L and a upper triangular matrix U, as described in LU decomposition.
- A = LU
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 L are set to 1
- l11 = 1
- l22 = 1
- l33 = 1
so we get a solvable system of 9 unknowns and 9 equations.
Solving for the other l and u, we get the following equations:
- u11 = a11
- u12 = a12
- u13 = a13
- u22 = a22 − u12l21
- u23 = a23 − u13l21
- u33 = a33 − (u13l31 + u23l32)
and for l:
We see that there is a calculation pattern, which can be expressed as the following formulas, first for U
and then for L
We see in the second formula that to get the lij below the diagonal, we have to divide by the diagonal element (pivot) uij, so we get problems when uij is either 0 or very small, which leads to numerical instability.
The solution to this problem is pivoting A, which means rearranging the rows of A, prior to the LU decomposition, in a way that the largest element of each column gets onto the diagonal of A. Rearranging the columns means to multiply A by a permutation matrix P:
Example:
The decomposition algorithm is then applied on the rearranged matrix so that
- PA = LU
Task description
The task is to implement a routine which will take a square nxn matrix A and return a lower triangular matrix L, a upper triangular matrix U and a permutation matrix P, 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
Contents |
[edit] Ada
decomposition.ads:
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;
decomposition.adb:
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;
Example usage:
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;
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
[edit] BBC BASIC
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$
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
[edit] C
Compiled withgcc -std=gnu99 -Wall -lm -pedantic. Demonstrating how to do LU decomposition, and how (not) to use macros. #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;
}
[edit] Common Lisp
Uses the routine (mmul A B) from Matrix multiplication.
;; 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)))
Example 1:
(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))
Example 2:
(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))
[edit] 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*/ {
return m.all!(r => r.length == m[0].length);
}
bool isSquare(T)(in T[][] m) /*pure nothrow*/ {
return isRectangular(m) && 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(isSquare(m));
} body {
immutable n = m.length;
auto id = iota(n)
.map!(j=> n.iota.map!(i => cast(T)(i == j)).array)
.array;
foreach (immutable i; 0 .. n) {
// immutable row = iota(i, n).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(isSquare(A));
} 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;
}
const P = pivotize!T(A);
const 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 = std.array.replicate("[%([%(%.1f, %)],\n %)]]\n\n", 3);
foreach (m; [a, b])
writefln(f, lu(m).tupleof);
}
- 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]]
[edit] Fortran
program lu1
implicit none
real*8 :: a1(3,3), a2(4,4)
a1 = reshape((/real*8::1,2,1,3,4,1,5,7,0/),(/3,3/))
call check(a1)
a2 = reshape((/real*8::11,1,3,2,9,5,17,5,24,2,18,7,2,6,1,1/),(/4,4/))
call check(a2)
contains
subroutine lu(a,p)
! in situ decomposition, correspondes to LAPACK's dgebtrf
implicit none
real*8, intent(inout) :: a(:,:)
integer, intent(out) :: p(:)
integer :: n, i,j,k,ii
n = ubound(a,1)
p = (/ ( i, i=1,n ) /)
do k = 1,n-1
ii = k-1+maxloc(abs(a(p(k:),k)),1)
if (ii /= k ) then
p((/k, ii/)) = p((/ii, k/))
end if
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 forall
end do
end subroutine
subroutine check(a)
implicit none
real*8, intent(in) :: a(:,:)
real*8 :: aa(ubound(a,1), ubound(a,2))
real*8 :: l(ubound(a,1), ubound(a,2))
real*8 :: u(ubound(a,1), ubound(a,2))
integer :: p(ubound(a,1), ubound(a,2)), ipiv(ubound(a,1))
integer :: i, n
character(len=100) :: fmt
n = ubound(a,1)
aa = a ! work with a copy
p = 0; l=0; u = 0
forall (i=1:n)
p(i,i) = 1d0; l(i,i) = 1d0 ! convert permutation vector a matrix
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
write (fmt,"(a,i1,a)") "(",n,"(f8.2,1x))"
print *, "a"
print fmt, transpose(a)
print *, "p"
print fmt, transpose(dble(p))
print *, "l"
print fmt, transpose(l)
print *, "u"
print fmt, transpose(u)
print *, "residual"
print *, "|| P.A - L.U || = ", maxval(abs(matmul(p,a)-matmul(l,u)))
end subroutine
end program>
- 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
[edit] Go
[edit] 2D representation
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")
}
- 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
[edit] Flat representation
package main
import "fmt"
type matrix struct {
ele []float64
stride int
}
func matrixFromRows(rows [][]float64) *matrix {
if len(rows) == 0 {
return &matrix{nil, 0}
}
m := &matrix{make([]float64, len(rows)*len(rows[0])), len(rows[0])}
for rx, row := range rows {
copy(m.ele[rx*m.stride:(rx+1)*m.stride], row)
}
return m
}
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{make([]float64, (len(m1.ele)/m1.stride)*m2.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{make([]float64, rows*cols), 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(matrixFromRows([][]float64{
{1, 3, 5},
{2, 4, 7},
{1, 1, 0}}))
showLU(matrixFromRows([][]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")
}
Output is same as from 2D solution.
[edit] Library
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)
}
- 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}
[edit] J
Taken with slight modification from [1].
Solution:
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
Example use:
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
[edit] 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}
Outputs:
[edit] MATLAB / Octave
LU decomposition is part of language
A = [
1 3 5
2 4 7
1 1 0];
[L,U,P] = lu(A)
gives the 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:
A = [
11 9 24 2
1 5 2 6
3 17 18 1
2 5 7 1 ];
[L,U,P] = lu(A)
gives the 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
[edit] Creating a MATLAB function
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
[edit] 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]) */
[edit] 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: 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)
b = [[11,9,24,2],[1,5,2,6],[3,17,18,1],[2,5,7,1]]
for part in lu(b):
pprint(part)
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]]
[edit] 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]])
[edit] REXX
/*REXX pgm makes a matrix from input, performs/shows LU decomposition. */
#=0; P.=0; PA.=0; L.=0; U.=0 /*initialize some variables to 0.*/
parse arg x /*get the matrix elements from CL*/
call makeMat /*make the A matrix from 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 done.*/
/*──────────────────────────────────er subroutine───────────────────────*/
er: say; say '***error!***'; say; say arg(1); say; exit 13
/*──────────────────────────────────makeMat subroutine──────────────────*/
makeMat: ?=words(x); do N=1 for ?; if N**2==? then leave; end
if N**2\==? then call er 'not correct number of elements entered: ' ?
do r=1 for N /*build the "A" matrix from 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 subroutine──────────────────*/
manLmat: 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 subroutine──────────────────*/
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 subroutine──────────────────*/
manUmat: 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 subroutine──────────────────*/
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
end /*j*/
end /*i*/
return
/*──────────────────────────────────showMat subroutine──────────────────*/
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
end
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; say _
end
return
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
[edit] 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-1).inject(0.0) {|sum, k| sum + u[k][j] * l[i][k]}
else
# lower
l[i][j] = (tmp[i,j] - (0 .. j-1).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)
each_with_index do |val, i, j|
print "#{format} " % val
puts "" if j==column_size-1
end
end
end
a = Matrix[[1, 3, 5],
[2, 4, 7],
[1, 1, 0]]
puts "A"; a.pretty_print("%2d")
l, u, p = a.lu_decomposition
puts "U"; u.pretty_print("%8.5f")
puts "L"; l.pretty_print("%8.5f")
puts "P"; p.pretty_print("%d")
a = Matrix[[11, 9,24,2],
[ 1, 5, 2,6],
[ 3,17,18,1],
[ 2, 5, 7,1]]
puts "A"; a.pretty_print("%2d")
l, u, p = a.lu_decomposition
puts "U"; u.pretty_print("%8.5f")
puts "L"; l.pretty_print("%8.5f")
puts "P"; p.pretty_print("%d")
output
A 1 3 5 2 4 7 1 1 0 U 2.00000 4.00000 7.00000 0.00000 1.00000 1.50000 0.00000 0.00000 -2.00000 L 1.00000 0.00000 0.00000 0.50000 1.00000 0.00000 0.50000 -1.00000 1.00000 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 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 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 P 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1
[edit] 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
}
}
Support code:
# 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
}
}
Demonstrating:
# This does the decomposition and prints it out nicely
proc demo {A} {
lassign [matrix::luDecompose $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}}
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
