LU decomposition: Difference between revisions

=={{header|AutoHotkey}}==

<lang AutoHotkey>;--------------------------

LU_decomposition(A){

P := Pivot(A)

A_ := Multiply_Matrix(P, A)

U := [], L := [], n := A_.Count()

loop % n {

i := A_Index

loop % n {

j := A_Index

Sigma := 0, k := 1

while (k <= i-1)

Sigma += (U[k, j] * L[i, k]), k++

U[i, j] := A_[i, j] - Sigma

Sigma := 0, k := 1

while (k <= j-1)

Sigma += (U[k, j] * L[i, k]), k++

L[i, j] := (A_[i, j] - Sigma) / U[j, j]

}

}

return [L, U, P]

}

;--------------------------

Pivot(M){

n := M.Count(), P := [], i := 0

while (i++ < n){

P.push([])

j := 0

while (j++ < n)

P[i].push(i=j ? 1 : 0)

}

i := 0

while (i++ < n){

maxm := M[i, i], row := i, j := i

while (j++ < n)

if (M[j, i] > maxm)

maxm := M[j, i], row := j

if (i != row)

tmp := P[i], P[i] := P[row], P[row] := tmp

}

return P

}

;--------------------------

Multiply_Matrix(A,B){

if (A[1].Count() <> B.Count())

return

RCols := A[1].Count()>B[1].Count()?A[1].Count():B[1].Count()

RRows := A.Count()>B.Count()?A.Count():B.Count(), R := []

Loop, % RRows {

RRow:=A_Index

loop, % RCols {

RCol:=A_Index, v := 0

loop % A[1].Count()

col := A_Index, v += A[RRow, col] * B[col,RCol]

R[RRow,RCol] := v

}

}

return R

}

;--------------------------

ShowMatrix(L, f:=3){

for r, obj in L{

row := ""

for c, v in obj

row .= Format("{:." f "f}", v) ", "

output .= "[" trim(row, ", ") "]`n,"

}

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

}

;--------------------------</lang>

Examples:<lang AutoHotkey>A1 := [[1, 3, 5]

, [2, 4, 7]

, [1, 1, 0]]

A2 := [[11, 9, 24, 2]

,[1, 5, 2, 6]

,[3, 17, 18, 1]

,[2, 5, 7, 1]]

loop 2 {

L := LU_Decomposition(A%A_Index%)

result .= ""

. "A:=`n" ShowMatrix(A%A_Index%, 4)

. "`n`nL:=`n" ShowMatrix(L.1)

. "`n`nU:=`n" ShowMatrix(L.2)

. "`n`nP:=`n" ShowMatrix(L.3)

. "`n--------------------------------`n"

}

MsgBox, 262144, , % result

return</lang>

{{out}}

<pre>A:=

[[1.0000, 3.0000, 5.0000]

,[2.0000, 4.0000, 7.0000]

,[1.0000, 1.0000, 0.0000]]

L:=

[[1.000, 0.000, 0.000]

,[0.500, 1.000, 0.000]

,[0.500, -1.000, 1.000]]

U:=

[[2.000, 4.000, 7.000]

,[0.000, 1.000, 1.500]

,[0.000, 0.000, -2.000]]

P:=

[[0.000, 1.000, 0.000]

,[1.000, 0.000, 0.000]

,[0.000, 0.000, 1.000]]

--------------------------------

A:=

[[11.0000, 9.0000, 24.0000, 2.0000]

,[1.0000, 5.0000, 2.0000, 6.0000]

,[3.0000, 17.0000, 18.0000, 1.0000]

,[2.0000, 5.0000, 7.0000, 1.0000]]

L:=

[[1.000, 0.000, 0.000, 0.000]

,[0.273, 1.000, 0.000, 0.000]

,[0.091, 0.287, 1.000, 0.000]

,[0.182, 0.231, 0.004, 1.000]]

U:=

[[11.000, 9.000, 24.000, 2.000]

,[0.000, 14.545, 11.455, 0.455]

,[0.000, 0.000, -3.475, 5.688]

,[0.000, 0.000, 0.000, 0.511]]

P:=

[[1.000, 0.000, 0.000, 0.000]

,[0.000, 0.000, 1.000, 0.000]

,[0.000, 1.000, 0.000, 0.000]

,[0.000, 0.000, 0.000, 1.000]]

--------------------------------</pre>