Arjen Markus Just a little experiment with a common form of matrix multiplication: transform coordinates in a 3D space. This type of operation arises when you want to do 3D graphics for instance. It can easily be adjusted to any fixed number of dimensions, but if you have serious needs for matrix multiplication or matrix manipulation in general, I suggest you look at Ed Hume's LA package.
# matrices.tcl --
# Simple demonstration of matrix-vector multiplications
#
# The idea:
# In 3D space:
# - A vector is represented as a list of 3 numbers
# - A matrix is represented as a list of 3 lists of 3 numbers
# - A series of vectors is also a list of lists
# For instance:
# The identity matrix is { {1 0 0} {0 1 0} {0 0 1}}
# Rotation over X degrees around the z-axis is implemented
# via:
# { {C -S 0}
# {S C 0}
# {0 0 1}}
# where C = cos X, S = sin X
#
# Rotating the vector {0 1 1} thus gives:
# { [expr {C*0-S*1+0*1}] [expr {S*0+C*1+0*1}] [expr {0*0+0*1+1*1}] }
# = {-S C 1}
#
namespace eval ::matrices {
namespace export matmul
}
# matmul --
# Multiply a vector (or list of vectors) by a 3x3 matrix
#
# Arguments:
# matrix Matrix to multiply with
# vector Vector operand (or list of vectors)
# Result:
# Transformed vector (vector list)
# Note:
# A vector list of 3 vectors is identical to a 3x3 matrix,
# hence this procedure can be used for matrix multiplications too
# Transposition is implied by the implementation (just to remind
# myself :)
#
proc ::matrices::matmul {matrix vector} {
set n1 11
set n2 12
set n3 13
foreach row $matrix {
foreach [list m$n1 m$n2 m$n3] $row { break }
incr n1 10
incr n2 10
incr n3 10
}
#
# Add an extra level of list if there is only one vector
#
if { [llength [lindex $vector 0]] == 1 } {
set vector [list $vector]
}
set r [list]
foreach v $vector {
foreach {x y z} $v { break }
lappend r [list [expr {$m11*$x+$m12*$y+$m13*$z}] \
[expr {$m21*$x+$m22*$y+$m23*$z}] \
[expr {$m31*$x+$m32*$y+$m33*$z}] ]
}
return $r
}
# Test code
#
namespace import ::matrices::*
set vectors { {1.0 0.0 0.0} {0.0 1.0 0.0} {0.0 0.0 1.0} }
set matrix { {1.0 2.0 3.0} {4.0 5.0 6.0} {7.0 8.0 9.0} }
set expected { {1.0 4.0 7.0} {2.0 5.0 8.0} {3.0 6.0 9.0} }
puts "Per vector:"
foreach vector $vectors exp $expected {
set result [matmul $matrix $vector]
puts "$result -- $exp"
}
puts "For a list of vectors:"
set result [matmul $matrix $vectors]
puts "$result -- $expected"
puts "Matrix multiplication (matrix**2 is identity)"
set matrix { {0.0 1.0 0.0} {1.0 0.0 0.0} {0.0 0.0 1.0} }
set result [matmul $matrix $matrix]
puts "$result"AMG: For my computer graphics class I developed some 3D routines and a nifty model display application. (I might post it here after the class is done.) Concerned about speed, and not wanting to depend on anything other than Tcl and Tk, I wrote the following. This unreadable hunk of metaprogramming produces matrix multipliers with fully unrolled loops.
# Assemble a fast matrix multiply procedure. All matrices are 3x4 and are
# treated as if they had a fourth row of [0 0 0 1].
set code ""
for {set y 0} {$y < 3} {incr y} {
set row ""
for {set x 0} {$x < 4} {incr x} {
set cell ""
for {set i 0} {$i < 3} {incr i} {
if {$i != 0} {
append cell " + "
}
append cell "\[lindex \$m1 $y $i\] * \[lindex \$m2 $i $x\]"
}
if {$x == 3} {
append cell " + \[lindex \$m1 $y $x\]"
}
if {$x != 0} {
append row " "
}
append row "\[expr [list $cell]\]"
}
if {$y != 0} {
append code " "
}
append code "\[list $row\]"
}
proc 3dmtx_compose {m args} [string map [list %CODE% $code] {
if {[llength $args] == 0} {
return $m
} else {
set m2 [lindex $args end]
set matrices [concat [list $m] [lrange $args 0 end-1]]
for {set i 0} {$i < [llength $matrices]} {incr i} {
set m1 [lindex $matrices end-$i]
set m2 [list %CODE%]
}
return $m2
}
}]
# Put together a fast matrix-vector multiply procedure. The matrix is 3x4
# and is combined with a fourth row of [0 0 0 1], and the vector is 1x3 and
# is transposed and combined with a fourth row of [1].
set code ""
for {set y 0} {$y < 3} {incr y} {
set cell ""
for {set i 0} {$i < 3} {incr i} {
if {$i != 0} {
append cell " + "
}
append cell "\[lindex \$m $y $i\] * \[lindex \$v $i\]"
}
append cell " + \[lindex \$m $y $i\]"
if {$y != 0} {
append code " "
}
append code "\[expr [list $cell]\]"
}
proc 3dmtx_apply {args} [string map [list %CODE% $code] {
switch -- [llength $args] {
0 {error "wrong # args: should be \"3dmtx_apply ?matrix? ?...? vector\""}
1 {return [lindex $args 0]}
2 {set m [lindex $args 0]}
default {set m [3dmtx_compose {*}[lrange $args 0 end-1]]}
}
set v [lindex $args end]
return [list %CODE%]
}][3dmtx_compose] multiplies any number of 4x4 affine transformation matrices*, each of which is stored as 3x4 list (three elements in the outer list; each element is itself a list containing four numbers). The fourth row is taken to be (0 0 0 1).
[3dmtx_apply] multiplies (the composition of) any number of such matrices by a 3D column vector in homogeneous coordinates, represented by a list of three numbers. The fourth row is implicitly (1).
* Affine transformation matrices can scale, shear, rotate, mirror, and translate. Parallel lines remain parallel after affine transformations. 3x3 matrices, as used by Arjen's code, cannot perform translation. 4x4 transformation matrices are constructed by adding an extra translation column (x y z) to the right of an Arjen-style 3x3 matrix, then adding a new row (0 0 0 1) for the sake of homogeneous coordinates.
Hmm, I probably should have profiled the "loopy" version first. I plan to update this section to show a comparision between loopy multiplies and unrolled multiplies. Maybe there's no benefit. I don't know yet.
So you can play around with my multipliers, here's code to generate useful transformation matrices.
# Returns a matrix which translates by ($x,$y,$z).
proc 3dmtx_trans {x y z} {
return [list [list 1 0 0 $x]\
[list 0 1 0 $y]\
[list 0 0 1 $z]]
}
# Returns a matrix which rotates by $t radians about vector ($x,$y,$z).
proc 3dmtx_rotate {x y z t} {
if {$x == 0 && $y == 0 && $z == 0} {
error "Cannot rotate around zero vector"
}
set l [expr {sqrt($x ** 2 + $y ** 2 + $z ** 2)}]
foreach dim {x y z} {
set $dim [expr {[set $dim] / $l}]
}
set result [list]
foreach row_expr {
{{$x * $x + (1 - $x ** 2) * cos($t)}
{$x * $y * (1 - cos($t)) - $z * sin($t)}
{$x * $z * (1 - cos($t)) + $y * sin($t)}}
{{$y * $x * (1 - cos($t)) + $z * sin($t)}
{$y * $y + (1 - $y ** 2) * cos($t)}
{$y * $z * (1 - cos($t)) - $x * sin($t)}}
{{$z * $x * (1 - cos($t)) - $y * sin($t)}
{$z * $y * (1 - cos($t)) + $x * sin($t)}
{$z * $z + (1 - $z ** 2) * cos($t)}}
} {
set row [list]
foreach cell_expr [concat $row_expr [list 0]] {
lappend row [expr $cell_expr]
}
lappend result $row
}
return $result
}
# Returns a matrix which scales by ($x,$y,$z). Negative scalars mirror.
proc 3dmtx_scale {x y z} {
return [list [list $x 0 0 0]\
[list 0 $y 0 0]\
[list 0 0 $z 0]]
}And now an example.
% set pi_2 [expr {acos(-1) / 2}]
% 3dmtx_apply [3dmtx_trans -1 0 0] \
[3dmtx_rotate 1 0 0 $pi_2]\
[3dmtx_scale 2 2 2] \
{1 2 3}
1.0 -6.0 4.0Let's trace that. (1,2,3) scaled by 2 becomes (2,4,6) rotated by 90° about X becomes (2,-6,4) translated by -X becomes (1,-6,4). Nifty. And what's the matrix we used?
% 3dmtx_compose [3dmtx_trans -1 0 0] \
[3dmtx_rotate 1 0 0 $pi_2]\
[3dmtx_scale 2 2 2]
{{ 2.0 0.0 0.0 -1.0}
{ 0.0 0.0 -2.0 0.0}
{ 0.0 2.0 0.0 0.0}}I took the liberty of breaking the list into three lines and replacing 1.2246063538223773e-16 with 0.0 so you can more easily see what's going on. And remember, all matrices have an implicit fourth row of (0 0 0 1).
On my 500MHz P3 laptop, the above composition takes 1097.3528 microseconds, and that includes time spent generating the inputs. Applying that matrix to a vector takes 48.9898 microseconds. Composing and applying together take 1149.24 microseconds.
That's enough for now. It's time for me to go to class.