LU decomposition: Difference between revisions

=={{header|Tcl}}==

<lang tcl>package require Tcl 8.5

namespace eval matrix {

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

proc I {order} {

set m [lrepeat $order [lrepeat $order 0]]

for {set i 0} {$i < $order} {incr i} {

lset m $i $i 1

}

return $m

}

proc pivotize {matrix} {

set n [llength $matrix]

set p [I $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} {

set tmp [lindex $p $j]

lset p $j [lindex $p $row]

lset p $row $tmp

}

}

return $p

}

# Decomposes a square matrix A by PA=LU and returns 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

}

# Code adapted from Matrix_multiplication and Matrix_transposition tasks

proc size {m} {

set rows [llength $m]

set cols [llength [lindex $m 0]]

return [list $rows $cols]

}

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

}

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

}

proc print {matrix {fmt "%s"}} {

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

}

}

}

# Demonstration helper

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}}</lang>

<pre>

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

</pre>