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 (wrapfort) 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] AM (6 january 2010) Here is a small program (rather incomplete, but still) that scans the LAPACK source files and produces a first version of an interface definition that can be used by wrapfort
# autowrap.tcl --
# Program to automatically generate the wrapfort-style
# wrapper definitions for LAPACK
#
# TODO:
# - string arguments
# - input/output arguments
# - functions
# - recognising matrix arguments (dlarrv: matrix z for instance)
#
# extractInterface --
# Extract the interface for a LAPACK/BLAS routine
#
# Arguments:
# filename Name of the source file to examine
#
# Returns:
# Description of the interface as a wrapfort command
#
proc extractInterface {filename} {
set infile [open $filename r]
#
# First part: the subroutine/function header
#
set header ""
while {1} {
gets $infile line
append header $line
if { [string first ) $line] >= 0 } {
break
}
}
set header [string tolower $header]
#
# Second part: the description of the arguments
#
while { [gets $infile line] >= 0 } {
if { [string first "* Arguments" $line] >= 0 } {
gets $infile line
break
}
}
set arguments {}
set count 0
while { [gets $infile line] >= 0 } {
if { [string first "======" $line] >= 0 } {
break
}
set line [string tolower $line]
switch -- [lindex $line 2] {
"(input)" {
incr count
set name [lindex $line 1]
set type [extractType $line]
lappend arguments [list $type $name input]
}
"(output)" {
incr count
set name [lindex $line 1]
set type [extractType $line]
set size [extractSize $line]
if { $name != "info" } {
lappend arguments [list $type $name result]
} else {
lappend arguments [list $type $name {assign 0}]
}
}
"(workspace)" {
incr count
set name [lindex $line 1]
set type [extractType $line]
set size [extractSize $line]
lappend arguments [list $type $name [list allocate $size]]
}
default {
# Ignore
}
}
}
foreach {routine numberArgs call rtype} [convertHeader $header $arguments] {break}
if { $count != $numberArgs } {
puts "File: $filename -- not all arguments analysed ($count -- $numberArgs)"
}
lappend arguments [list code {} $call]
return [list [string toupper $routine] [string tolower $routine] $arguments $rtype]
}
# extractType --
# Extract the data type of an argument
#
# Arguments:
# line Line of text describing the argument
#
# Returns:
# Keyword for wrapfort
#
proc extractType {line} {
set type "?"
switch -glob -- $line {
"*integer array*" {
set type "integer-array"
}
"*integer*" {
set type "integer"
}
"*double precision array*" {
set type "double-array"
}
"*double precision*" {
set type "double"
}
"*real array*" {
set type "real-array"
}
"*real*" {
set type "real"
}
"*character*" {
set type "string"
}
"*logical*" {
set type "logical"
}
}
return $type
}
# extractSize --
# Extract the size of an array argument
#
# Arguments:
# line Line of text describing the argument
#
# Returns:
# Size information useable by wrapfort
#
proc extractSize {line} {
set size "?"
regexp {dimension \((.*)\)} $line ==> size
return $size
}
# convertHeader --
# Convert the subroutine/function header
#
# Arguments:
# header Complete Fortran header
# arguments Information about the arguments (types notably)
#
# Returns:
# Size information useable by wrapfort
#
proc convertHeader {header arguments} {
regexp { ([a-z0-9]+)\(} $header ==> routine
set count [llength [split $header ,]]
# For now
set rtype ""
#
# Check the individual arguments
#
set call "$routine\("
set wrapup ""
set listargs {}
foreach arg $arguments {
foreach {type name size} $arg {break}
if { [string match *-array $type] } {
lappend listargs "$name"
} else {
lappend listargs "&$name"
}
if { $name == "info" } {
set wrapup "
if ( info != 0 ) {
Tcl_SetResult( interp, \"Error in $routine\", NULL );
return TCL_ERROR;
}"
}
}
set call "\n$call [join $listargs ", "] );\n$wrapup"
return [list $routine $count $call $rtype]
}
# main --
# Get it all going
#
set outfile [open "autowrap_lapack.tcl" w]
foreach file [glob d*.f] {
puts "File: $file ..."
foreach {routine cname arguments rtype} [extractInterface $file] {break}
puts $outfile "\nWrapfort::fproc ::Lapack::$routine $cname {"
foreach arg $arguments {
puts $outfile " $arg"
}
puts $outfile "}"
}The result - for a few, arbitrary, routines is this:
Wrapfort::fproc ::Lapack::DLAED6 dlaed6 {
integer kniter input
logical orgati input
double rho input
double-array d input
double-array z input
double finit input
double tau result
integer info {assign 0}
code {} {
dlaed6( &kniter, &orgati, &rho, d, z, &finit, &tau, &info );
if ( info != 0 ) {
Tcl_SetResult( interp, "Error in dlaed6", NULL );
return TCL_ERROR;
}}
}
Wrapfort::fproc ::Lapack::DLARRV dlarrv {
integer n input
double vl input
double vu input
integer-array isplit input
integer m input
integer dol input
integer dou input
double minrgp input
double rtol1 input
double rtol2 input
integer-array iblock input
integer-array indexw input
double-array gers input
double-array z result
integer ldz input
integer-array isuppz result
double-array work {allocate 12*n}
integer-array iwork {allocate 7*n}
integer info {assign 0}
code {} {
dlarrv( &n, &vl, &vu, isplit, &m, &dol, &dou, &minrgp, &rtol1, &rtol2, iblock, indexw, gers, z, &ldz, isuppz, work, iwork, &info );
if ( info != 0 ) {
Tcl_SetResult( interp, "Error in dlarrv", NULL );
return TCL_ERROR;
}}
}
Wrapfort::fproc ::Lapack::DSGESV dsgesv {
integer n input
integer nrhs input
integer lda input
integer-array ipiv result
double-array b input
integer ldb input
double-array x result
integer ldx input
double-array work {allocate n*nrhs}
real-array swork {allocate n*(n+nrhs)}
integer iter result
integer info {assign 0}
code {} {
dsgesv( &n, &nrhs, &lda, ipiv, b, &ldb, x, &ldx, work, swork, &iter, &info );
if ( info != 0 ) {
Tcl_SetResult( interp, "Error in dsgesv", NULL );
return TCL_ERROR;
}}
}These interfaces are not useable as yet:
Despite these shortcomings I am making progress :).