Keith Vetter 20051102 - today I needed to do some linear programming on equations of three unknowns and hence, wanted to invert some 3x3 matrices. For matrices this small, it's easiest just to compute the inverse as the
adjoint matrix divided by the determinant
.
Here's the code--it works fine for my toy application, but the code ignores all the thorny issues of numerical stability that plague general matrix code.
proc Inverse3 {matrix} {
if {[llength $matrix] != 3 ||
[llength [lindex $matrix 0]] != 3 ||
[llength [lindex $matrix 1]] != 3 ||
[llength [lindex $matrix 2]] != 3} {
error "wrong sized matrix"
}
set inv {{? ? ?} {? ? ?} {? ? ?}}
# Get adjoint matrix : transpose of cofactor matrix
for {set i 0} {$i < 3} {incr i} {
for {set j 0} {$j < 3} {incr j} {
lset inv $i $j [_Cofactor3 $matrix $i $j]
}
}
# Now divide by the determinant
set det [expr {double([lindex $matrix 0 0] * [lindex $inv 0 0]
+ [lindex $matrix 0 1] * [lindex $inv 1 0]
+ [lindex $matrix 0 2] * [lindex $inv 2 0])}]
if {$det == 0} {
error "non-invertable matrix"
}
for {set i 0} {$i < 3} {incr i} {
for {set j 0} {$j < 3} {incr j} {
lset inv $i $j [expr {[lindex $inv $i $j] / $det}]
}
}
return $inv
}
proc _Cofactor3 {matrix i j} {
array set COLS {0 {1 2} 1 {0 2} 2 {0 1}}
foreach {row1 row2} $COLS($j) break
foreach {col1 col2} $COLS($i) break
set a [lindex $matrix $row1 $col1]
set b [lindex $matrix $row1 $col2]
set c [lindex $matrix $row2 $col1]
set d [lindex $matrix $row2 $col2]
set det [expr {$a*$d - $b*$c}]
if {($i+$j) & 1} { set det [expr {-$det}]}
return $det
}
