LU decomposition: Difference between revisions

 

 

The task is to implement a routine which will take a square nxn matrix <math>A</math> and return a lower triangular matrix <math>L</math>, a upper triangular matrix <math>U</math> and a permutation matrix <math>P</math>,

You should then test it on the following two examples and include your output.

 

<pre>

A

</pre>

 

<pre>

A

0 1 0 0

0 0 0 1

</pre>

 

{{works with|Ada 2005}}

decomposition.ads:

generic

with package Matrix is new Ada.Numerics.Generic_Real_Arrays (<>);

procedure Decompose (A : Matrix.Real_Matrix; P, L, U : out Matrix.Real_Matrix);

 

 

decomposition.adb:

 

procedure Swap_Rows (M : in out Matrix.Real_Matrix; From, To : Natural) is

end Decompose;

 

 

Example usage:

with Ada.Text_IO;

with Decomposition;

Ada.Text_IO.Put_Line ("U:"); Print (U_2);

Ada.Text_IO.Put_Line ("P:"); Print (P_2);

 

{{out}}

0.00 1.00 0.00 0.00

0.00 0.00 0.00 1.00</pre>

 

=={{header|BBC BASIC}}==

A1() = 1, 3, 5, 2, 4, 7, 1, 1, 0

PROCLUdecomposition(A1(), L1(), U1(), P1())

a$ = LEFT$(LEFT$(a$)) + CHR$(13) + CHR$(10)

NEXT i%

{{out}}

<pre>

 

=={{header|C}}==

#include <stdlib.h>

#include <math.h>

 

return 0;

 

=={{header|Common Lisp}}==

Uses the routine (mmul A B) from [[Matrix multiplication]].

 

(defun eye (n)

(let ((I (make-array `(,n ,n) :initial-element 0)))

 

;; Return L, U and P.

 

Example 1:

 

#2A((1 3 5) (2 4 7) (1 1 0))

 

#2A((1 0 0) (1/2 1 0) (1/2 -1 1))

#2A((2 4 7) (0 1 3/2) (0 0 -2))

 

Example 2:

 

#2A((11 9 24 2) (1 5 2 6) (3 17 18 1) (2 5 7 1))

 

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

 

=={{header|D}}==

{{trans|Common Lisp}}

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

 

foreach (immutable m; [a, b])

writefln(f, lu(m).tupleof);

{{out}}

<pre>[[1.0, 0.0, 0.0],

 

=={{header|EchoLisp}}==

(lib 'matrix) ;; the matrix library provides LU-decomposition

(decimals 5)

0 0 -3.475 5.6875

0 0 0 0.51079

 

=={{header|Fortran}}==

implicit none

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

call lu(aa, ipiv)

do i = 1,n

end do

p(ipiv,:) = p

do k = 1,n-1

kmax = maxloc(abs(a(p(k:),k)),1) + k-1

end do

end subroutine

end subroutine

 

{{out}}

<pre>

|| P.A - L.U || = 0.0000000000000000

</pre>

 

=={{header|Go}}==

===2D representation===

{{trans|Common Lisp}}

 

import "fmt"

u.print("u")

p.print("p")

{{out}}

<pre>

</pre>

===Flat representation===

 

import "fmt"

u.print("u")

p.print("p")

Output is same as from 2D solution.

 

 

import (

"fmt"

 

)

 

func main() {

1, 3, 5,

2, 4, 7,

}))

fmt.Println()

11, 9, 24, 2,

1, 5, 2, 6,

}

 

lu.Factorize(a)

fmt.Println("p:", lu.Pivot(nil))

{{out}}

Pivot format is a little different here. (But library solutions don't really meet task requirements anyway.)

 

===Library go.matrix===

 

import (

fmt.Printf("u:\n%v\n", u)

fmt.Printf("p:\n%v\n", p)

{{out}}

<pre>

</pre>

 

 

=={{header|Idris}}==

 

'''Solution:'''

module Main

 

putStrLn "Solution 2:"

printEx ex2

 

{{out}}

 

'''Solution:'''

 

LU=: 3 : 0

 

permtomat=: 1 {.~"0 -@>:@:/:

 

'''Example use:'''

LUdecompose A

┌─────┬─────┬───────┐

1 5 2 6

3 17 18 1

 

=={{header|Java}}==

Translation of [[#Common_Lisp|Common Lisp]] via [[#D|D]]

{{works with|Java|8}}

import java.util.Locale;

import static java.util.stream.IntStream.range;

print(m);

}

<pre> 1.0 0.0 0.0

0.5 1.0 0.0

0.0 1.0 0.0 0.0

0.0 0.0 0.0 1.0 </pre>

 

=={{header|jq}}==

 

'''Infrastructure'''

def matrix(m; n; init):

if m == 0 then []

| reduce range (0;$length) as $i

(""; . + reduce range(0;$length) as $j

'''LU decomposition'''

# Use "range(0;$n) as $i" to handle ill-conditioned cases.

def pivotize:

| . + [ $P ]

;

'''Example 1''':

a | lup[] | neatly(4)

{{Out}}

1 0 0

0.5 1 0

1 0 0

0 0 1

'''Example 2''':

{{Out}}

1 0 0 0

0.2727272727272727 1 0 0

0 0 1 0

0 1 0 0

 

'''Example 3''':

# A|lup|verify(A) should be true

def verify(A):

[0, 0, 1, -1]];

 

{{out}}

true

Julia has the predefined functions `lu`, `lufact` and `lufact!` in the standard library to compute the lu decomposition of a matrix.

{{Out}}

(

3x3 Array{Float64,2}:

 

=={{header|Kotlin}}==

 

typealias Vector = DoubleArray

printMatrix("U:", u2, "%8.5f")

printMatrix("P:", p2, "%1.0f")

 

{{out}}

0 1 0 0

0 0 0 1

</pre>

 

=={{header|Maple}}==

A:=<<1.0|3.0|5.0>,<2.0|4.0|7.0>,<1.0|1.0|0.0>>:

 

LinearAlgebra:-LUDecomposition(A);

{{out}}

<pre>

[0 0 1] [0.500000000000000 -1. 1.0] [0. 0. -2.]

</pre>

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

 

LUDecomposition(A);

{{out}}

<pre>

</pre>

 

{lu, p, c} = LUDecomposition[a];

l = LowerTriangularize[lu, -1] + IdentityMatrix[Length[p]];

P = Part[IdentityMatrix[Length[p]], p] ;

MatrixForm /@ {P.a , P, l, u, l.u}

{{out}}

[[File:LUex1.png]]

[[File:LUex2.png]]

 

=={{header|MATLAB}} / {{header|Octave}}==

LU decomposition is part of language

 

1 3 5

2 4 7

1 1 0];

 

{{out}}

<pre>

</pre>

2nd example:

11 9 24 2

1 5 2 6

2 5 7 1 ];

 

{{out}}

<pre>

 

===Creating a MATLAB function===

function [ P, L, U ] = LUdecomposition(A)

 

end

 

=={{header|Maxima}}==

 

 

a: hilbert_matrix(4)$

 

lu_backsub(lup, transpose([1, 1, -1, -1]));

 

=={{header|PARI/GP}}==

 

{

my (L = matid(#M), U = M, P = L);

 

[L,U,P] \\ return L,U,P triple matrix

 

Output:

</pre>

 

 

}

 

my $n = +@a;

} else {

}

}

}

}

 

my $size = +@m;

my $row = $i;

$row = $j;

}

}

}

@id

}

 

 

}

 

}

{{out}}

 

 

 

 

 

 

 

 

 

</pre>

 

=={{header|PL/I}}==

LU_Decomposition: procedure options (main);

declare a1(3,3) float (18) initial ( 1, 3, 5,

 

end LU_Decomposition;

Derived from Fortran version above.

Results:

=={{header|Python}}==

{{trans|D}}

 

def matrixMul(A, B):

for part in lu(b):

pprint(part)

{{out}}

<pre>[[1.0, 0.0, 0.0],

[0.0, 1.0, 0.0, 0.0],

[0.0, 0.0, 0.0, 1.0]]</pre>

 

=={{header|R}}==

{{Out}}

<pre>

$L

3 x 3 Matrix of class "dtrMatrix" (unitriangular)

[3,] . . |

</pre>

 

=={{header|Racket}}==

#lang racket

(require math)

; #[0 -2 -3]

; #[0 0 -2]])

 

=={{header|REXX}}==

parse arg x /*obtain matrix elements from the C.L. */

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

/*──────────────────────────────────────────────────────────────────────────────────────*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

<pre>

──A matrix───

0 0 -2

</pre>

<pre>

─────A matrix──────

 

=={{header|Ruby}}==

 

class Matrix

l.pretty_print(" %8.5f", "L")

u.pretty_print(" %8.5f", "U")

 

{{out}}

 

Matrix has a <code>lup_decomposition</code> built-in method.

l.pretty_print(" %8.5f", "L")

u.pretty_print(" %8.5f", "U")

Output is the same.

 

=={{header|Sidef}}==

func matrix_zero(n, m=n) { m.of { n.of(0) } }

func pivotize(m) {

var size = m.len

var id = matrix_ident(size)

var max = m[i][i]

var row = i

if (m[j][i] > max) {

max = m[j][i]

}

}

id.swap(row, i)

}

return id

}

func mmult(a, b) {

var p = []

}

return p

}

func lu(a) {

is_square(a) || die "Defined only for square matrices!";

var L = matrix_ident(n)

var U = matrix_zero(n)

if (j >= i) {

} else {

}

}

return [P, Aʼ, L, U]

}

func say_it(message, array) {

say "\n#{message}"

}

}

var t = [[

%n(1 3 5),

%n( 2 5 7 1),

]]

say_it(a, b)

}

<pre style="height: 40ex; overflow: scroll">

A Matrix

See [http://www.stata.com/help.cgi?mf_lud LU decomposition] in Stata help.

 

: lud(a=(1,3,5\2,4,7\1,1,0),l=.,u=.,p=.)

 

2 | 1 |

3 | 3 |

 

=== Implementation ===

real scalar i,j,n,s

real vector js

u = uppertriangle(l)

l = lowertriangle(l, 1)

 

'''Example''':

 

: a

2 | 1 |

3 | 3 |

 

=={{header|Tcl}}==

namespace eval matrix {

namespace path {::tcl::mathfunc ::tcl::mathop}

return $s

}

Support code:

namespace eval matrix {

# Get the size of a matrix; assumes that all rows are the same length, which

return $max

}

Demonstrating:

proc demo {A} {

lassign [matrix::luDecompose $A] L U P

demo {{1 3 5} {2 4 7} {1 1 0}}

puts "================================="

{{out}}

<pre>

0 1 0 0

0 0 0 1

</pre>

 

=={{header|zkl}}==

Using the GNU Scientific Library, which does the decomposition without returning the permutations:

fcn luTask(A){

A.LUDecompose(); // in place, contains L & U

2.0, 5.0, 7.0, 1.0);

L,U:=luTask(A);

{{out}}

<pre>

 

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

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

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

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

ans

Example 1

{{out}}

<pre>

</pre>

Example 2

L( 1.0, 5.0, 2.0, 6.0),

L( 3.0, 17.0, 18.0, 1.0),

 

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

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

{{out}}