Arjen Markus (28 december 2009) I have been contemplating this for several years now, but I finally made a start with a wrapper extension to get the functionality of LAPACK into Tcl.
I started off with two simple BLAS routines (BLAS is the more basic library that underlies LAPACK - basic linear algebra system, if I am not mistaken) and the Fortran wrapper you can find in my ftcl library (see [L1 ]). Immediately I ran into design issues:
Take the routine dapxy, which computes a*X + Y (a a scalar, X and Y vectors) and returns the result in Y. The interface is such that X and Y could also be two-dimensional arrays (matrices) so that you could use this routine to update the rows and columns of a matrix Y using a vector X o a matrix X.
Questions that arise:
Right now, my answers are:
set result [daxpy $a $x $y]
(The Fortran interface is:
call dapxy( n, a, dx, incx, dy, incy )
Because of incx and incy, you can pass two-dimensional arrays such that only a part is updated.)
Another question is: should we wrap all routines (some 1500!) or just the more commonly used ones? (If so, which ones?)
AM (31 december 2009) A discussion in the Tcl chatroom brought back an old idea: use byte arrays (binary strings) to store the vector or matrix data for a Fortran or C routine directly. Here is a small Tcl program that illustrates the basic storage mechanism (for vectors only at the moment):
# vecmat.tcl --
# Utilities to manipulate Tcl lists that hold
# binary data that can directly act as Fortran or C arrays
#
# TODO: error checking
#
# Vecmat --
# Namespace to hold these utilities
#
namespace eval ::Vecmat {
namespace export vector matrix
}
# vector --
# Utility to manipulate vectors
#
# Arguments:
# method Method in question
# args Any arguments for the method
#
# Returns:
# Depends on the method
#
proc ::Vecmat::vector {method args} {
switch -- $method {
"create" {
if { [llength $args] != 1 } {
return -code error "Invalid number of arguments - should be: vector create number"
}
return [::Vecmat::Vcreate {*}$args]
}
"set" {
if { [llength $args] != 3 } {
return -code error "Invalid number of arguments - should be: vector set name index number"
}
return [::Vecmat::Vset {*}$args]
}
"tolist" {
if { [llength $args] != 1 } {
return -code error "Invalid number of arguments - should be: vector tolist vector-value"
}
return [::Vecmat::Vtolist {*}$args]
}
"fromlist" {
if { [llength $args] != 1 } {
return -code error "Invalid number of arguments - should be: vector fromlist list"
}
return [::Vecmat::Vfromlist {*}$args]
}
default {
return -code error "Unknown method: $method"
}
}
}
# Vcreate --
# Create a vector of doubles of a given length
#
# Arguments:
# number Number of elements
#
# Returns:
# Vector with values initialised to 0.0
#
proc ::Vecmat::Vcreate {number} {
set data [binary format d* [lrepeat $number 0.0]]
return [list VECTOR $number $data]
}
# Vset --
# Set a single element in a vector
#
# Arguments:
# name Name of the vector
# index Index of the element
# value The new value
#
# Returns:
# Vector with new value at position index
#
# Note:
# This updates the vector variable whose name is $name
#
proc ::Vecmat::Vset {name index value} {
upvar 2 $name name_
set data [lindex $name_ 2]
set data [binary format a*@[expr {8*$index}]d $data $value]
lset name_ 2 $data
return $name_
}
# Vtolist --
# Convert a vector into a regular Tcl list
#
# Arguments:
# vector value of the vector
#
# Returns:
# List containing the values of the vector
#
proc ::Vecmat::Vtolist {vector} {
binary scan [lindex $vector 2] d* data
return $data
}
# Vfromlist --
# Create a vector from a regular Tcl list of doubles
#
# Arguments:
# list list of values to be converted
#
# Returns:
# Vector containing the values of the list
#
proc ::Vecmat::Vfromlist {list} {
set data [binary format d* $list]
return [list VECTOR [llength $list] $data]
}
# main --
# Test the vector and matrix utilities
#
set v [::Vecmat::vector create 10]
foreach index {0 4} {
::Vecmat::vector set v $index [expr {cos($index)}]
}
puts [::Vecmat::vector tolist $v]
set w [::Vecmat::vector fromlist {1.0 .2 3.0}]
puts [::Vecmat::vector tolist $w]