[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 [http://local.wasp.uwa.edu.au/~pbourke/miscellaneous/determinant/index.html%|%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 } ====== <> Mathematics