LU decomposition: Difference between revisions

{{trans|Python}}

 

L(row) m

print(row)

L(part) lu(b)

pprint(part)

 

{{out}}

{{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|AutoHotkey}}==

LU_decomposition(A){

P := Pivot(A)

return "[" Trim(output, "`n,") "]"

}

, [2, 4, 7]

, [1, 1, 0]]

}

MsgBox, 262144, , % result

{{out}}

<pre>A:=

 

=={{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|C++}}==

#include <cmath>

#include <iomanip>

 

return 0;

 

{{out}}

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

end subroutine

 

{{out}}

<pre>

===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.

 

===Library gonum/mat===

 

import (

fmt.Printf("u: %.5f\n\n", mat.Formatted(u, mat.Prefix(" ")))

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>

=={{header|Haskell}}==

''Without elem-at-index modifications; doesn't find maximum but any non-zero element''

import Data.List

import Data.Maybe

main = putStrLn "Task1: \n" >> solveTask mY1 >>

putStrLn "Task2: \n" >> solveTask mY2

{{out}}

<pre>

===With Numeric.LinearAlgebra===

 

 

a1, a2 :: Matrix R

main = do

print $ lu a1

{{out}}

<pre>((3><3)

 

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

=={{header|Javascript}}==

{{works with|ES5 ES6}}

const mult=(a, b)=>{

let res = new Array(a.length);

return [...lu(A),P];

}

 

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

 

=={{header|Kotlin}}==

 

typealias Vector = DoubleArray

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

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

 

{{out}}

 

=={{header|Lobster}}==

 

// derived from JAMA v1.03

print_L A

print_U A

{{out}}

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

 

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

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

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|Nim}}==

 

For display, we use the third party module "strfmt" which allows to specify dynamically the format.

import strfmt

 

l2.print("L:", "8.5f")

u2.print("U:", "8.5f")

 

{{out}}

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

 

{

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

 

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

 

Output:

=={{header|Perl}}==

{{trans|Raku}}

 

for $test (

print "$line\n";

}

{{out}}

<pre style="height:35ex">A matrix

=={{header|Phix}}==

{{trans|Kotlin}}

<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>

<span style="color: #008080;">function</span> <span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span>

<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lu</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]),{</span><span style="color: #004600;">pp_Nest</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_Pause</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">})</span>

<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>

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

 

=={{header|R}}==

A <- matrix(c(1, 3, 5, 2, 4, 7, 1, 1, 0), 3, 3, byrow=T)

dim(A) <- c(3, 3)

 

{{Out}}

 

=={{header|Racket}}==

#lang racket

(require math)

; #[0 -2 -3]

; #[0 0 -2]])

 

=={{header|Raku}}==

 

Translation of Ruby.

[2, 4, 7],

[1, 1, 0]

say "\n$message";

$_».&rat-int.fmt("%7s").say for @array;

{{out}}

<pre>A Matrix

 

=={{header|REXX}}==

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

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

end /*c*/

say _

{{out|output|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 1 3 5 &nbsp; &nbsp; 2 4 7 &nbsp; &nbsp; 1 1 0 </tt>}}

<pre>

 

=={{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|Rust}}==

{{libheader| ndarray}}

#![allow(non_snake_case)]

use ndarray::{Array, Axis, Array2, arr2, Zip, NdFloat, s};

(L, U, P)

}

 

 

[0, 0, 0, 1]]

</pre>

 

=={{header|Sidef}}==

{{trans|Raku}}

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>

=={{header|VBA}}==

{{trans|Phix}}

Private Function pivotize(m As Variant) As Variant

Dim n As Integer: n = UBound(m)

Debug.Print

Next i

<pre>== a,l,u,p: ==

1 3 5

{{libheader|Wren-matrix}}

{{libheader|Wren-fmt}}

 

var arrays = [

Fmt.mprint(lup[2], 2, 0)

System.print()

 

{{out}}

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