[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 - '''b'''asic '''l'''inear 
'''a'''lgebra '''s'''ystem, if I am not mistaken) and the Fortran wrapper (''wrapfort'') you can find in my ftcl library (see [http://ftcl.sf.net]).
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:

   * Is that flexibility something we want in the Tcl command as well?
   * Do we want to update a vector or should the Tcl command simpy return the result?
   * Do we want to expose the full interface? That is the most flexible, but it is also distinctly non-Tcl

Right now, my answers are:

   * We merely want to pass two vectors X and Y, not columns or rows of a matrix
   * We do not want to update any arguments - so the Tcl command will return the resulting vector and it is up to the user to put it into the original Y argument, if he/she wants to.
   * The interface should be as simple as:

   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:
   * For the routines dsgesv and dlarrv a number of arguments are input/output and the program does not support that yet
   * If there are more than one output arguments, the result on the Tcl side should be a list of these arguments. Again something I still need to implement.

Despite these shortcomings I am making progress :).

<<categories>>Category Mathematics|Category Numerical Analysis