[Arjen Markus] (1 april 2004 and dead-serious!) I wanted to get some feeling for the ease by which you can do matrix operations in Tcl (I mean: the kind of things that you deal with in Linear Algebra). So the script below explores three methods of representation and the performance you get with them. Since two of these methods are very similar I used a different approach in the matrix*vector computation for the third. Warning: the results here are for a single matrix operation only, so it is much too early to draw any conclusions. If, and only if, these results are representative, then the first implementation, using lists of lists, is the fastest and easiest to program. However, the matrix*vector computation used here is very suitable for that structure. Matrix*matrix computations would require more jumping in the two levels of lists, so probably the other two would be faster then (Still to be measured!) ---- # testmat.tcl -- # Test the speed of several matrix implementations: # The code below is meant to investigate the consequences # of three different implementations of a matrix data type: # 1. As a list of lists # 2. As a flat list with the first few elements representing the size # 3. As a pair of lists, one the matrix sizes and the other the matrix # as a flat list # The algorithm is simple: multiple a vector by the matrix # The matrix is simple: # 1 2 3 ... N-1 0 # 2 3 ..... 0 1 # 3 4 ..... 1 2 # ... # The vector is too: (1,2,3,4,...) # Note: all numbers are double precision reals! # catch {console show} proc makemat1 {N} { set result {} for { set j 0 } { $j < $N } { incr j } { set row {} for { set i 0 } { $i < $N } { incr i } { lappend row [expr {double(($i+$j+1)%$N)}] } lappend result $row } return $result } proc makemat2 {N} { set result [list 2 $N $N] ;# Two-dimensional for { set j 0 } { $j < $N } { incr j } { for { set i 0 } { $i < $N } { incr i } { lappend result [expr {double(($i+$j+1)%$N)}] } } return $result } proc makemat3 {N} { set result [list [list $N $N]] set values {} for { set j 0 } { $j < $N } { incr j } { for { set i 0 } { $i < $N } { incr i } { lappend values [expr {double(($i+$j+1)%$N)}] } } lappend result $values } proc makevect {N} { set result {} for { set i 0 } { $i < $N } { incr i } { lappend result [expr {double($i+1)}] } return $result } proc matmul1 { matrix vector } { set newvect {} foreach row $matrix { set sum 0.0 foreach v $vector c $row { set sum [expr {$sum+$v*$c}] } lappend newvect $sum } return $newvect } proc matmul2 { matrix vector } { foreach {dummy ncol nrow} $matrix {break} set newvect {} set idx 3 for { set j 0 } { $j < $nrow } { incr j } { set sum 0.0 foreach v $vector { set c [lindex $matrix $idx] set sum [expr {$sum+$v*$c}] incr idx } lappend newvect $sum } return $newvect } proc matmul3 { matrix vector } { foreach {ncol nrow} [lindex $matrix 0] {break} set values [lindex $matrix 1] set newvect {} for { set j 0 } { $j < $nrow } { incr j } { set idx1 [expr {$j*$ncol}] set idx2 [expr {$idx1+$ncol-1}] set sum 0.0 foreach v $vector c [lrange $values $idx1 $idx2] { set sum [expr {$sum+$v*$c}] } lappend newvect $sum } return $newvect } # # Test and measure # puts "Test the matrix construction" puts [set matrix1 [makemat1 4]] puts [set matrix2 [makemat2 4]] puts [set matrix3 [makemat3 4]] puts "Test the matrix multiplication" set matrix1 [makemat1 10] set matrix2 [makemat2 10] set matrix3 [makemat3 10] set vector [makevect 10] puts [set newvect1 [matmul1 $matrix1 $vector]] puts [set newvect2 [matmul2 $matrix2 $vector]] puts [set newvect3 [matmul3 $matrix3 $vector]] puts "Measure the relative speeds" foreach size {3 10 30 100} { puts "Size: $size" set matrix1 [makemat1 $size] set matrix2 [makemat2 $size] set matrix3 [makemat3 $size] set vector [makevect $size] puts [time {set newvect1 [matmul1 $matrix1 $vector]} 100] puts [time {set newvect2 [matmul2 $matrix2 $vector]} 100] puts [time {set newvect3 [matmul3 $matrix3 $vector]} 100] } ---- The result (on a Windows XP machine with Tcl 8.4.1): Test the matrix construction {1.0 2.0 3.0 0.0} {2.0 3.0 0.0 1.0} {3.0 0.0 1.0 2.0} {0.0 1.0 2.0 3.0} 2 4 4 1.0 2.0 3.0 0.0 2.0 3.0 0.0 1.0 3.0 0.0 1.0 2.0 0.0 1.0 2.0 3.0 {4 4} {1.0 2.0 3.0 0.0 2.0 3.0 0.0 1.0 3.0 0.0 1.0 2.0 0.0 1.0 2.0 3.0} Test the matrix multiplication 285.0 250.0 225.0 210.0 205.0 210.0 225.0 250.0 285.0 330.0 285.0 250.0 225.0 210.0 205.0 210.0 225.0 250.0 285.0 330.0 285.0 250.0 225.0 210.0 205.0 210.0 225.0 250.0 285.0 330.0 Measure the relative speeds Size: 3 14 microseconds per iteration 16 microseconds per iteration 21 microseconds per iteration Size: 10 85 microseconds per iteration 102 microseconds per iteration 105 microseconds per iteration Size: 30 650 microseconds per iteration 864 microseconds per iteration 728 microseconds per iteration Size: 100 6998 microseconds per iteration 9193 microseconds per iteration 7340 microseconds per iteration ---- [Lars H]: Here is an alternative matmul3: proc matmul3b { matrix vector } { foreach {ncol nrow} [lindex $matrix 0] {break} set values [lindex $matrix 1] set newvect [list] for {set i 0; set i_flat 0} {$i < $nrow} {incr i} { set sum 0.0 foreach v $vector { set sum [expr {$sum + $v*[lindex $values $i_flat]}] incr i_flat } lappend newvect $sum } return $newvect } Interestingly enough, I find that this is faster than matmul3 for matrices of side 3 and 10, but slower for matrices of side 30 and 100. Apparently [lrange]+[foreach] outperforms [lindex]+[incr] on long lists. That's curious (and would be an argument in favour of the first matrix format, if the operations are supposed to be implemented in Tcl rather than by a compiled extension). [AM] Yes, it goes against the somewhat anecdotal evidence I have. So I think it is important to do these measurements systematically .... ---- [[ [Category Numerical Analysis] ]]