Updated 2012-05-13 05:12:43 by RLE

Here's another weekend fun project from Richard Suchenwirth - a gadget for matrix operations. I call a gadget a "poor man's object" - it has a clear, friendly API, but doesn't require an OO system, not even namespaces. Sort of like widgets in Tk. One might also call a gadget a value-added variable.

Matrices are represented as lists of lists, each sublist standing for a row. Row/col indices begin with 1. Operations ("methods") include so far:

  • constructor: matrix name ?value?
  • get dimensions: name dim
  • get the whole matrix: name ($M and [M] are equivalent)
  • set the whole matrix: name set value (see below)
  • get one column/one row: name col|row n
  • modify a row: name row n add|mult|swap n2
  • get or set an element: name elem row col ?value?
  • transpose (swap rows and columns): name transpose
  • formatting: name format ;# dirt simple, just join with \n
  • destructor: name delete (or unset name)

See usage examples at end. Assignment can be done by both
 M set {1 2 3} {4 5 6}
 set M {{1 2 3} {4 5 6}}

but in the first case you can omit the outer braces, and the matrix name is returned, while in the second you get your value back. Pretty minimal but fun. Enjoy ;-)
 proc matrix {_M args} {
        if {$_M=="-"} {
                upvar 2 [lindex $args 0] M
                set argl [lindex $args 1]
                set cmd [lindex  $argl 0]
                set para [lrange $argl 1 end]
                switch -glob -- $cmd {
                col*   {return [matrix:col M [lindex $para 0]]}
                delete {unset M}
                dim*   {return [list [llength $M] [llength [lindex $M 0]]]}
                el*    {return [matrix:element M $para]}
                format {return [join $M \n]}
                row    {return [eval matrix:row M $para]}
                set     {set M $para; return [lindex $args 0]}
                transp* {return [matrix:transpose M]}
                ""      {return $M}
                default {return -code error "unknown command $cmd"}
                }
        } else {
                #---------------- constructor
                upvar $_M M
                if [info exists M] {return -code error "$_M exists"}
                proc $_M args "matrix - $_M \$args" ;# helper proc
                trace variable M u "rename $_M {};#"
                if [llength $args]==1 {set args [lindex $args 0]}
                set M $args ;# assume only data coming..
                return $_M ;# own name, for stacking 
        }
 }

In order to call a matrix by its name, a one-liner helper proc is created for each (and deleted on unset) that calls the dispatcher.
 proc matrix:col {_M col} {
        upvar $_M M
        incr col -1
        if {$col<0 || $col>[llength [lindex $M 0]]} {
                return -code error "bad column $col in $_M"
        }
        set res [list]
        foreach row $M {lappend res [lindex $row $col]}
        set res
 }
 proc matrix:element {_M dim} {
        upvar $_M M
        foreach {rows cols value} $dim break
        set el [lindex [lindex $M [incr rows -1]] [incr cols -1]]
        if {$el==""} {return -code error "bad indices $dim in $_M"}
        if {$value!=""} {
                set row [lindex $M $rows]
                set row [lreplace $row $cols $cols [set el $value]]
                set tmp [lreplace $M $rows $rows $row]
                set M $tmp ;# direct assignment crashes Tcl 8.1a2
        }
        set el
 }
 proc matrix:row {_M args} {
    #extended row handler, allows swap/add/mult operations
    upvar $_M M
    if [llength $args]==1 {
        set res [lindex $M [incr args -1]]
        if {$res==""} {return -code error "bad row [incr args] in $_M"}
        return $res
    } else {
        foreach {rowno cmd value} $args break
        set row1 [lindex $M [incr rowno -1]]
        switch -- $cmd {
        add {
                set row2 [lindex $M [incr value -1]]
                set row [list]
                foreach x $row1 y $row2  {
                        lappend row [expr {$x+$y}]
                }
                set tmp [lreplace $M $rowno $rowno $row]
                set M $tmp ;# direct assignment crashes Tcl 8.1a2                               leted on unset) that calls the dispatcher.}
        }
        mult {
                set row [list]
                foreach el $row1 {lappend row [expr {$el*$value}]}
                set tmp [lreplace $M $rowno $rowno $row]
                set M $tmp ;# direct assignment crashes Tcl 8.1a2                               leted on unset) that calls the dispatcher.}
        }
        swap {
                set row2 [lindex $M [incr value -1]]
                set tmp1 [lreplace $M $rowno $rowno $row2]
                set M [lreplace $tmp1 $value $value $row1]
        }
        default {return -code error "$cmd? Must be add, mult, or swap"}
        } 
    }
 }

In transposition I make use of Tcl's unique feature that numbers can be variable names, so incr can generate them...
 proc matrix:transpose _M {
        upvar $_M M
        foreach row $M {
                set c 0
                foreach element $row {lappend $c $element; incr c}
        }
        set res [list]
        for {set c 0} {$c<[llength $row]} {incr c} {
                lappend res [set $c]
        }
        set res
 }
 ##################################### test code, usage examples
 catch {unset M; unset MT}        ;# for repeated sourcing
 matrix M {1 2 3} {4 5 6} {7 8 9} ;# construct and initialize
 puts dim=[M dim]                 ;# list of {nrows ncols}
 puts row2=[M row 2]
 puts r3*10=[M row 3 mult 10]     ;# multiply a row by a scalar
 puts r1Xr2=[M row 1 swap 2]       ;# swap rows 1 and 2
 puts r3+r1=[M row 3 add 1]        ;# add row 1 to row 3
 M element 2 3 5.5                ;# changing an element
 [matrix MT [M transpose]] format ;# stacking three operations

See also tally: a string counter gadget - Gadgets for a generalized approach - lexpr for simple arithmetics on pure-value matrices

FF 2007-06-06 : Nice! If there were also a method to get a submatrix (often called the minor matrix around ij, where i the row, j the column), then it would be possible to do compute determinant using Cofactor Expansion [1] -- here's my implementation:
 proc matrix:minor {_M args} {
    upvar $_M M
    if {[llength $args] == 2} {
        lassign $args i j
        set ci 0
        set c [list]
        foreach row $M {
                incr ci
                if {$i == $ci} {continue}
                set tmprow [list]
                set cj 0
                foreach element $row {
                        incr cj
                        if {$j == $cj} {continue}
                        lappend tmprow $element
                }
                if {[llength $tmprow] > 0} {lappend c $tmprow}
        }
        set c
    }
 }

and here's my implementation of determinant:
 proc matrix:determinant _M {
    # compute along the first column, but it can be done better
    upvar $_M M
    lassign [M dim] rows cols
    if {$rows == $cols} {
        if {$rows == 1} {return $M}
        set res 0
        set i 0
        set sign 1
        foreach row $M {
                incr i
                incr res [expr $sign*[lindex $row 0]*[[matrix TEMP [M minor $i 1]] det]]
                set sign [expr ($sign==1)?-1:1]
                unset TEMP
        }
        return $res
    } else {
        puts "works only for NxN matrices!"
    }
 }

beware, it has j! (factorial) complexity (where j the number of rows = sqrt(n)... really considering also the 'minor' method, complexity is (n^(3/2))!