LU decomposition

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:

A={\begin{pmatrix}a_{11}&a_{21}&a_{31}\\a_{21}&a_{22}&a_{32}\\a_{31}&a_{32}&a_{33}\\\end{pmatrix}}={\begin{pmatrix}l_{11}&0&0\\l_{21}&l_{22}&0\\l_{31}&l_{32}&l_{33}\\\end{pmatrix}}{\begin{pmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{pmatrix}}=LU

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

l_{11}=1

l_{22}=1

l_{33}=1

so we get a solvable system of 9 unknowns and 9 equations.

A={\begin{pmatrix}a_{11}&a_{21}&a_{31}\\a_{21}&a_{22}&a_{32}\\a_{31}&a_{32}&a_{33}\\\end{pmatrix}}={\begin{pmatrix}1&0&0\\l_{21}&1&0\\l_{31}&l_{32}&1\\\end{pmatrix}}{\begin{pmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{pmatrix}}={\begin{pmatrix}u_{11}&u_{12}&u_{13}\\u_{11}l_{21}&u_{12}l_{21}+u_{22}&u_{13}l_{21}+u_{23}\\u_{11}l_{31}&u_{12}l_{31}+u_{22}l_{32}&u_{13}l_{31}+u_{23}l_{32}+u_{33}\end{pmatrix}}=LU

Solving for the other $l$ and $u$ , we get the following equations:

u_{11}=a_{11}

u_{12}=a_{12}

u_{13}=a_{13}

u_{22}=a_{22}-u_{12}l_{21}

u_{23}=a_{23}-u_{13}l_{21}

u_{33}=a_{33}-(u_{13}l_{31}+u_{23}l_{32})

and for $l$ :

l_{21}={\frac {1}{u_{11}}}a_{21}

l_{31}={\frac {1}{u_{11}}}a_{31}

l_{32}={\frac {1}{u_{22}}}(a_{32}-u_{12}l_{31})

We see that there is a calculation pattern, which can be expressed as the following formulas, first for $U$

u_{ij}=a_{ij}-\sum _{k=1}^{i-1}u_{kj}l_{ik}

and then for $L$

l_{ij}={\frac {1}{u_{jj}}}(a_{ij}-\sum _{k=1}^{j-1}u_{kj}l_{ik})

We see in the second formula that to get the $l_{ij}$ below the diagonal, we have to divide by the diagonal element (pivot) $u_{ij}$ , so we get problems when $u_{ij}$ 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$ :

PA\Rightarrow A'

Example:

{\begin{pmatrix}0&1\\1&0\end{pmatrix}}{\begin{pmatrix}1&4\\2&3\end{pmatrix}}\Rightarrow {\begin{pmatrix}2&3\\1&4\end{pmatrix}}

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

Ada

Works with: Ada 2005

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

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

Translation of: Common Lisp

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

      std.array, std.conv;

bool isRectangular(T)(in T[][] m) pure /*nothrow*/ {

   return !canFind!((r){ return r.length != m[0].length; })(m);

}

bool isSquare(T)(in T[][] m) pure /*nothrow*/ {

   return isRectangular(m) && m[0].length == m.length;

}

immutable(T[][]) matrixMul(T)(immutable T[][] A,

                             immutable T[][] B) pure nothrow
   in {
       assert(A.length && B.length &&
              isRectangular(A) && isRectangular(B) &&
              A[0].length == B.length);
   } body {
       auto result = new T[][](A.length, B[0].length);
       auto aux = new T[B.length];
       foreach (j; 0 .. B[0].length) {
           foreach (k, row; B)
               aux[k] = row[j];
           foreach (i, ai; A)
               result[i][j] = dotProduct(ai, aux);
       }
       return result;
   }

/// Creates a nxn identity matrix. T[][] identityMatrix(T)(in int n) pure nothrow

   in {
       assert(n >= 0);
   } out(result) {
       assert(isSquare(result));
   } body {
       auto m = new typeof(return)(n, n);
       foreach (i, row; m) {
           row[] = 0;
           row[i] = 1;
       }
       return m;
   }

/// Creates the pivoting matrix for m. immutable(T[][]) pivotize(T)(immutable T[][] m) pure nothrow

   in {
       assert(isSquare(m));
   } body {
       immutable n = m.length;
       auto P = identityMatrix!T(n);

       foreach (i; 0 .. n) {
           T max = m[i][i];
           auto row = i;
           foreach (j; i .. n)
               if (m[j][i] > max) {
                   max = m[j][i];
                   row = j;
               }
           if (i != row)
               swap(P[i], P[row]);
       }

       return P;
   }

/// 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 (i; 0 .. n) {
           L[i][i .. $] = 0;
           U[i][0 .. i] = 0;
       }

       immutable P = pivotize!T(A);
       immutable A2 = matrixMul!T(P, A);

       foreach (j; 0 .. n) {
           L[j][j] = 1;
           foreach (i; 0 .. j+1) {
               T s1 = 0;
               foreach (k; 0 .. i)
                   s1 += U[k][j] * L[i][k];
               U[i][j] = A2[i][j] - s1;
           }
           foreach (i; j .. n) {
               T s2 = 0;
               foreach (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., 3, 5],
                  [2., 4, 7],
                  [1., 1, 0]];
   foreach (part; lu(a))
       writeln("[", join(map!text(part), ",\n "), "]\n");
   writeln();

   immutable b = [[11., 9, 24, 2],
                  [1.,  5,  2, 6],
                  [3., 17, 18, 1],
                  [2.,  5,  7, 1]];
   foreach (part; lu(b))
       writeln("[", join(map!text(part), ",\n "), "]\n");

}</lang> Output:

[[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]]


[[1, 0, 0, 0],
 [0.272727, 1, 0, 0],
 [0.090909, 0.2875, 1, 0],
 [0.181818, 0.23125, 0.00359712, 1]]

[[11, 9, 24, 2],
 [0, 14.5455, 11.4545, 0.454545],
 [0, 0, -3.475, 5.6875],
 [0, 0, 0, 0.510791]]

[[1, 0, 0, 0],
 [0, 0, 1, 0],
 [0, 1, 0, 0],
 [0, 0, 0, 1]]

Go

2D representation

Translation of: Common Lisp

<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 {

   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")

}</lang>

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>

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> Outputs:

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>

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: <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>

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

Python

Translation of: D

<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: 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]]

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-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")</lang>

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

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