Keith Vetter 2005-11-02 - 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
}DKF: Here's an alternate version of _Cofactor3 that should be a bit (~30%) faster (array setup costs are a little high and lassign is faster than foreach+break for the case it handles):
proc _Cofactor3 {matrix i j} {
set COLS {{1 2} {0 2} {0 1}}
lassign [lindex $COLS $j] row1 row2
lassign [lindex $COLS $i] col1 col2
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
}