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 {m v} "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 such a matrix by a 3D column vector in homogeneous coordinates, represented by a list containing three numbers. The fourth row is assumed to be 1.
* Affine transformation matrices can scale, shear, rotate, mirror, and translate. Parallel lines remain parallel after affine transformations.
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.