# matrix.tcl - matrix math routines # Jonas Beskow 2000 # the internal representation of a matrix is an array structured as below: # m(rows) - number of rows # m(cols) - number of columns # m(0,0) m(0,1) ... m(\$i,\$j) - elements of the matrix # # matrices are passed as lists (i.e. array get/set is used to # serialize/restore the matrix array) namespace eval matrix { } # mult - multiply two matrices proc matrix::mult {m1l m2l} { array set m1 \$m1l array set m2 \$m2l if {\$m1(cols) != \$m2(rows)} { error "matrix inner dimensions must agree" } for {set i 0} {\$i<\$m1(rows)} {incr i} { for {set j 0} {\$j<\$m2(cols)} {incr j} { set iprod 0 for {set k 0} {\$k<\$m1(cols)} {incr k} { set iprod [expr {\$iprod + \$m1(\$i,\$k)*\$m2(\$k,\$j)}] } set m(\$i,\$j) \$iprod } } set m(rows) \$m1(rows) set m(cols) \$m2(cols) array get m } # transpose - transpose a matrix proc matrix::transpose {ml} { array set m \$ml for {set i 0} {\$i<\$m(rows)} {incr i} { for {set j 0} {\$j<\$m(cols)} {incr j} { set n(\$j,\$i) \$m(\$i,\$j) } } set n(rows) \$m(cols) set n(cols) \$m(rows) array get n } # eye - generate a unit matrix of a given dimension proc matrix::eye {dim} { set m(rows) \$dim set m(cols) \$dim for {set i 0} {\$i<\$dim} {incr i} { for {set j 0} {\$j<\$dim} {incr j} { if {\$i==\$j} {set m(\$i,\$j) 1.0} {set m(\$i,\$j) 0.0} } } array get m } # det - calculate the determinant of a matrix proc matrix::det {ml} { array set m \$ml if {\$m(rows)!=\$m(cols)} {error "matrix must be square"} set n \$m(rows) set det 0 for {set i 0} {\$i<\$n} {incr i} {lappend colperm0 \$i} foreach {colperm alpha} [_permute2 \$colperm0] { set row 0 set prod 1 foreach col \$colperm { set prod [expr {\$prod*\$m(\$row,\$col)}] incr row } if {\$alpha%2==0} {set sign 1} {set sign -1} set det [expr {\$det+\$sign*\$prod}] } return \$det } # det2 - calculate the determinant of a matrix (recursive algorithm) # algorithm: # Choose a fixed row value i. # The determinant can be calculated emanating from the ith row. # |A| = Ai,1 * ai,1 + Ai,2 * ai,2 + Ai,3 * ai,3 + ... Ai,n * ai,n # The value of Ai,j = (-1)^(i+j) * (the determinant of the sub-matrix of A, # obtained from A by crossing out the ith row and the jth column) proc matrix::det2 {ml} { array set m \$ml set self [lindex [info level 0] 0] if {\$m(rows)!=\$m(cols)} {error "matrix must be square"} set n \$m(rows) if {\$n==1} {return \$m(0,0)} set i 0 set det 0 for {set j 0} {\$j<\$n} {incr j} { if {(\$i+\$j)%2==0} {set sign 1} {set sign -1} set subm [_wipe \$ml \$i \$j] set subdet [\$self \$subm] set det [expr {\$det + \$sign * \$subdet * \$m(\$i,\$j)}] } return \$det } # lset - convert from list to internal representation proc matrix::lset {list} { set i 0 foreach row \$list { puts row=\$row set j 0 foreach elem \$row { set m(\$i,\$j) \$elem incr j } incr i } set m(rows) \$i set m(cols) \$j array get m } # lget - convert from internal to list representation proc matrix::lget {ml} { array set m \$ml set list "" for {set i 0} {\$i<\$m(rows)} {incr i} { set row {} for {set j 0} {\$j<\$m(cols)} {incr j} { lappend row \$m(\$i,\$j) } puts row=\$row lappend list \$row } set list } # cset - convert from string to internal representation proc matrix::cset {s} { set list "" regsub -all {([[:blank:]])*(\n)+([[:blank:]])*} [string trim \$s] \n s puts s=\$s regsub -all {[[:blank:]]+} \$s " " s puts s=\$s foreach row [split \$s \n] { lappend list [split \$row " "] } lset \$list } # lget - convert from internal to string representation proc matrix::cget {ml} { array set m \$ml set str "" for {set i 0} {\$i<\$m(rows)} {incr i} { set row {} for {set j 0} {\$j<\$m(cols)} {incr j} { append str \$m(\$i,\$j) if {\$j!=[expr \$m(cols)-1]} {append str \t} } if {\$i!=[expr \$m(rows)-1]} {append str \n} } string trim \$str } # # utility routines to generate common 4x4 matrices for 3D transformation # proc matrix::translate {x y z} { lset [list \ "1 0 0 \$x" \ "0 1 0 \$y" \ "0 0 1 \$z" \ "0 0 0 1"] } proc matrix::scale {x y z} { lset [list \ "\$x 0 0 0" \ "0 \$y 0 0" \ "0 0 \$z 0" \ "0 0 0 1"] } proc matrix::rotX {a} { set ca [expr {cos(\$a)}] set sa [expr {sin(\$a)}] set nsa [expr {-\$sa}] lset [list \ "1 0 0 0" \ "0 \$ca \$nsa 0" \ "0 \$sa \$ca 0" \ "0 0 0 1"] } proc matrix::rotY {a} { set ca [expr {cos(\$a)}] set sa [expr {sin(\$a)}] set nsa [expr {-\$sa}] lset [list \ "\$ca 0 \$sa 0" \ "0 1 0 0" \ "\$nsa 0 \$ca 0" \ "0 0 0 1"] } proc matrix::rotZ {a} { set ca [expr {cos(\$a)}] set sa [expr {sin(\$a)}] set nsa [expr {-\$sa}] lset [list \ "\$ca \$nsa 0 0" \ "\$sa \$ca 0 0" \ "0 0 1 0" \ "0 0 0 1"] } # _permute2 # - recursive help function to calculate all possible permutations of a list # returns a list of the form {perm1 n1 perm2 n2 ... permN nN} where # permK is a sub-list containing a permutation and nK is the number of # exchanges necessary to bring the original sequence to permK proc matrix::_permute2 {list} { set self [lindex [info level 0] 0] if {[llength \$list]==1} {return [list \$list 0]} set result "" set first [lindex \$list 0] set tail [lrange \$list 1 end] foreach {p n} [\$self \$tail] { for {set i 0} {\$i<=[llength \$p]} {incr i} { lappend result [linsert \$p \$i \$first] [expr \$i+\$n] } } return \$result } # _wipe - help func to det2 proc matrix::_wipe {ml rr cc} { array set m \$ml set m2(rows) [expr {\$m(rows)-1}] set m2(cols) [expr {\$m(cols)-1}] for {set r 0} {\$r<\$m(rows)} {incr r} { if {\$r==\$rr} continue if {\$r>\$rr} {set r2 [expr {\$r-1}]} else {set r2 \$r} for {set c 0} {\$c<\$m(cols)} {incr c} { if {\$c==\$cc} continue if {\$c>\$cc} {set c2 [expr {\$c-1}]} else {set c2 \$c} set m2(\$r2,\$c2) \$m(\$r,\$c) } } array get m2 } return # Testing below # ----------------------------------------------------------------------------- # # magic squares set m10 [m::cset " 92 99 1 8 15 67 74 51 58 40 98 80 7 14 16 73 55 57 64 41 4 81 88 20 22 54 56 63 70 47 85 87 19 21 3 60 62 69 71 28 86 93 25 2 9 61 68 75 52 34 17 24 76 83 90 42 49 26 33 65 23 5 82 89 91 48 30 32 39 66 79 6 13 95 97 29 31 38 45 72 10 12 94 96 78 35 37 44 46 53 11 18 100 77 84 36 43 50 27 59"] set m7 [m::cset " 30 39 48 1 10 19 28 38 47 7 9 18 27 29 46 6 8 17 26 35 37 5 14 16 25 34 36 45 13 15 24 33 42 44 4 21 23 32 41 43 3 12 22 31 40 49 2 11 20"] set m5 [m::cset " 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9"] set m3 [m::cset " 8 1 6 3 5 7 4 9 2"] set a [m::lset [list {1 2} {3 4}]] set a [m::lset [list {1 2 3 4} {5 6 7 8} {9 10 11 12} {13 14 15 16}]] set b [m::lset [list {9 2 1 6} {6 4 7 23} {23 1 3 6} {98 5 3 1}]] puts [m::lget [m::mult \$a \$b]] set v [m::lset {1 0 0 1}] set T [m::translate 0 .5 0] set R [m::rotZ [expr 2*atan(1)]] puts [m::lget [m::mult \$R \$v]]