Numerical Analysis in Tcl

Numerical Analysis in Tcl provides information on tools available for numerical analysis, as well as discussion about development of such tools.

Reference edit

Numerical Analysis, Wikipedia

Tcl Tools edit

Geometry: a work-in-progress by Christian Heide Damm, Arjen Markus, and others. The idea is simple - define points, vectors, planes, lines, line segments and so on via "qualified" lists. Define operations on these lists like, calculate the distance between points, project a line on a plane.
L: provides more algebraic syntax for mathematical expressions
infix: more algebraic syntax for mathematical expressions
MATLAB
Octave: various interfaces to Tcl/Tk are available
NAP: 0-based indexing, deals with large collections of data in a way similar to MATLAB or Fortran's array operations
Tcllib's math: provides a lot of functionality
Vkit is a vector engine: a different approach to vectors in Tcl. It uses a descriptor to describe vectors, and can mix Tcl lists as well as 1/2/4/8/16/32-bit int and 32/64-bit floats transparently. There's a pure-Tcl version of the underlying "mvec" primitive, as well as an efficient C-based on which also supports direct reading off memory-mapped files. Please note that VKIT is work-in-progress.
A little math language
ASCEND IV
Reading and writing NetCDF files, AM
Simple implementation of Real Number
scilab
VecTcl: An implementation of a proposed vector math engine
Tensor
narray
la
HQP: (Huge Quadratic Programming) is a solver for nonlinearly constrained large-scale optimization written in C++. It is delivered as an extended Tcl interpreter. Tcl is used for selecting solver modules, configuring parameters and for controlling the execution.
Hmd Mathematics: a software package for numerical solution of ordinary differential equations (ODE) and finite elements (FE), featuring a Fortran/C style language interpreter for entering math statements into a console window. It has native array processing capability, like Matlab and Scilab. For the number theory researcher it has whole number, rational number processing capabilty (see Gnumbers). It has a Tcl/Tk graphical display and C/Tcl interpreter integration. The graphical display features a console interpreter window, a table (appears like a spreadsheet) window, and a plot window (line plots). A Run Script window provides a convenient for quickly executing script files.

Algorithms edit

Runge-Kutta: Methods of solving ordinary differential equations.

Description edit

As of 8.5, expr works with integers of arbitrary size.

Notes on Syntax edit

GWM 2006-02-04: Some discussion below suggested commands like "a = 2*3" should be implemented. I produced a prodding code (prodding others to do better, improve or contribute) and a good solution is now to be found in C-like syntax for numbers - this has uncovered some alternative implementations as referenced in the next paragraph.

slebetman notes that others (actually, mostly Richard Suchenwirth) have already done better. See: Gadgets, Radical Language Modification and Let unknown know.

GWM: On the other hand if you can't find it, it may as well not be there. The gadgets implementation is nice; the others are clever but too many self modifying parts for such a simple task make the code hard to verify.

Morten Skaarup Jensen writes on clt:

It's wonderful to see some people with an interest in a public math extension. I'd like to support the opinions that prime candidates for such an extension would be support for complex numbers, vectors and matrices (preferably n-dimensional) of both real and complex values. Lapack, Linpack etc would also greatly enhace the value of such a package. We also need to be able to evaluate expressions involving these new datatypes (slicing, arithmetic operators, mathematical functions, linear algebra functions etc) and I'm sure we need some new syntax for this (lindex $list [expr {[llength $list]}]-5 is just not good enough if we want people to prefer this over Python or MATLAB), but I think we should be constructive and thankful for any work done on such an extension rather than constantly picking on minor details of the syntax. I do however hope that the syntax will be well thought through from the beginning.

As I see it there are two possible styles:

A style similar to expr. In this case I hope it doesn't inherit all the peculiarities of C (e.g. bitwise operators could be replaced by functions; do we really need a?b:c).
A completely different style - in this case I would hope one draws inspiration from languages that are good at numerics (e.g. MATLAB). One might also get a better result by not allowing obnoxious variable/function names.

Joe English replies:

>As I see it there are two possible styles:

How would people react to an APL-like syntax, specifically one that doesn't have any notion of operator precedence and just evaluates everything from right to left? E.g., a * b + c would be interpreted as a * (b + c), not (a * b) + c.

Larry Smith This is actually a pretty cool idea since the whole thing could be had for the effort of bolting NARS2000 [1] into the interpreter via a load command that gives us an "apl" command, which could just pass it off to NARS2000 for evaluation and return the result. This interpreter is relatively small, not blindingly fast but reasonable, gets cutting-edge features added regularly, and has pretty good support.

Or, we could stick with Tcl syntax like math::* [math::+ $a $b $c] $d] and add a few new commands for subscripting and slicing.

Larry Smith Of course. If someone did it. But I wished to point out, incorporating the actual apl interpreter requires far less effort.

Morten replies:

I can only talk for myself and the people that I know who use numerics a lot: The further from normal math textbook syntax one gets, the less likely it will be that people will use it. Both of your proposals would be acceptable as a temporary step so that one at least can use the new datatypes, but the ultimate goal needs to be something close to textbook syntax.

Larry Smith Yes and no. APL is the oldest and still the most popular of the vector languages and very much forms a standard in itself. It also strikes me that it makes far more sense to add APL to Tcl, rather than try to make APL a visual language itself. It certainly fits the "little language" paradigm, even if it is an entire programming language. Putting the two together makes a lot of sense to me, adding the capabilities not in one to the other seems like re-inventing the wheel.

Larry Smith Just a few quick thoughts - a quick look at NARS2000 suggests it might be fairly trivial to make it run as a powerful math back end for Tcl, but one issue will need to be addressed. Although you can retrieve the value of a variable within the workspace by simply typing the var name and parsing what the interpreter spits back out at you, it would like be quite painful to do this for tightly coupled code, like in a simulation of some sort. Ideally, all the vars in the workspace should accessible from Tcl and vice versa, suggesting that the variable look-up code in NARS would need to be merged with Tcl's somehow - likely by adding an extension to Tcl_Obj with access information in the NARS interpreter, so that NARS variables could be accessed as if they were normal tcl vars in the same namespace as the interpreter (and thereby permitting sharing of NARS workspace vars with other namespaces using upvar). When NARS failed to find a var in the workspace rather than dumping out the DOMAIN ERROR msg, it should try to look it up in the Tcl namespace the workspace lives in and construct the link, so that access back and forth is transparent. Given the amount of math-grinding code written in apl, this might be a big win for Tcl in general, making it+NARS the logical successor to a lot of proprietary apl systems. It might encourage a whole new bunch of users to come in to the Tcl party. Each language specializes in making different programming needs very simple.

NARS2000 is an apl testbed for lots of new apl features, too. It can handle a large number of dimensions, new data types (quaternions and other esoterica) and, as an interpreter in itself, it can off-load a lot of grunt work that expr does in Tcl, and do it far more concisely - even more so than various Tcl9 suggestions to make it shorter. The special characters it uses can be accessed easily using \u notation even from an ascii file, and not much of an editor or input processor would be needed to support them directly.

Donal Fellows replies: If you've got an LR or LL parser in there, putting in precedence is not hard.

[Steve West] comments:

 > an extension would be support for complex numbers, vectors and matrices
 > (preferably n-dimensional) of both real and complex values. Lapack, Linpack
 > etc would also greatly enhance the value of such a package. We also need to
 > be able to evaluate expressions involving these new datatypes (slicing,

I'm jumping in the middle of this without seeing the whole thread, but I'm confused. People seem to be advocating re-inventing the wheel. I'll ask a naive question: why aren't we discussing registering the commands of Octave with tcl (treat it like an extension package)? Since the mid 70s, astronomers have plotted with a program called pgplot

. Some have registered all the pgplot commands within Tcl. This philosophy seems like a more straightforward approach (to integrate something which already exists). Since I'm oblivious as to just how to register outside routines and make them share memory and operate in the Tcl environment, I may be way off base here.

The main person developing Octave is stepping down, and there has been a lot of discussion about integrating some kind of GUI environment with Octave. Seems to me, tcl would a good fit...

OTOH, one has access to octave via Tcl through the former's command line interface. Piping or socketing information between the two scripting languages would seem to suffice for most things. The only problem with this approach is that you have to learn and maintain two different syntaxes.

bobs5, comp.lang.tcl, 2001-02:

I'm not sure if this suggestion qualifies as 'numerical support', but it seems to me that tcl expressions 'got ugly' when the need for bracing arrived. It just makes them hard to read and I have the hardest time remembering to put in the braces.

For normal scripting (not numerical analysis) I think it might be good to have a new command, call it assign as an abbreviation for assignment expression. For example:

assign {a  =  $b  +  $c}

assign takes only one argument, so you can't forget to put in the braces. If you use double quotes you get two substitution passes. Some people may want that sometimes. It might be possible to make assign handle multiple lines of expressions. Not sure.

This isn't a substitute for the math command mentioned elsewhere. That is more directed for numerical analysis.

Of course this can be done with uplevel in a proc, but I'm not sure if the compiler will optimize it.

Just my .02 cents.

bob

Donal Fellows replies:

I've no objection to standard Tcl expressions having an assignment operator, but much experience with C and other languages tells me that the operator must not be a single "=" sign (too easy to mix up with equality comparison, the source of many C bugs) but rather something like ":=". IMHO, this stuff with the syntax of the equality operator is probably the only really serious fault in the C language...

bobs5 writes:

You may be right about using =, but that is what the most common languages use. I would choose = because I think it will sell better. The one language that I can think of that uses := is pascal and I don't think it would be good establish a connection to pascal. I haven't seen any trouble people have had using =. What would that be? If its a problem learning to use the language then maybe it goes away after you've learned.

GWM 2007-02: The 'problem' is a commonly-written bug which compiles but does a different operation from planned. The following C/C++ lines show that a finger which only hits the = key once may cause an error.

if (c=1) { ....} // always runs the {} AND sets c to the value 1
if (c==1) { ....} // sometimes runs the {} if c equals 1

:= is one way of differentiating = from == by ensuring that bad typing should not result in the wrong operator. Tcl handles this by the set z ... and FORTRAN by if (c .eq. 1) then... . Perhaps we should look for a replacement for C's == rather than throw out the long established c=abc syntax. (End of gwm's insert)

Actually, I suggested adding an assign command, not operator. That way you don't change the expressions that are currently in use. Of course if you did add the operator, you could write:

if {a = $b + $c} {

instead of:

if {[set a [expr {$a + $b}]]} {

It looks a lot better to me.

I still think an assign command is needed because it accepts only one parameter, requiring people to use braces (or quotes).

If assign can do multi-line, then

set y [expr {$x + [f [expr {pow($x,$z)}]]}]

becomes

assign { exp = pow($x,$z)
         y = $x + [f $exp]
       }

Still not perfect, but somewhat better.

Of course the math command can have much better syntax since it doesn't have to keep the one expr uses.

Bob Techentin, comp.lang.tcl, 2001-02:

[email protected] wrote:

 > What type of numerical support would be useful to the community in general?
 > Perhaps someone could prototype an extension that would demonstrate what
 > is needed?

I, personally, would like n-dimensional matrices of complex numbers with a complete complement of matrix mathematics and linear algebra. Is that too tall an order? :-)

In my case, I use Rogue Wave's Math.h++ and Lapack.h++ libraries to develop application objects in C++. Then I use SWIG to wrap my application specific objects. So all my numerical code is essentially a specialized extension. I thought about directly wrapping the Rogue Wave library, but that would be a pretty flagrant violation of their license. :-(

A different approach is to simply script-enable a solid mathematics application like Octave or Matlab. Oops. Those are already scriptable. (Not Tcl, though.)

Another approach might be to build a good math tool or library and make it scriptable. Math Wizards ( http://www.mathwizards.com/

) takes that approach.

[Mumit Khan] replies:

In article <[email protected]>, Bob Techentin <[email protected]> wrote:

 >I, personally, would like n-dimensional matrices of complex numbers with
 >a complete complement of matrix mathematics and linear algebra.  Is that
 >too tall an order? :-)

You missed out on a few things, but I guess we have to start somewhere ;-)

 >In my case, I use Rogue Wave's Math.h++ and Lapack.h++ libraries to
 >develop application objects in C++.  Then I use SWIG to wrap my
 >application specific objects.

My experiment, many years ago (mostly during long flights when DVD wasn't pervasive on laptops ;-), involved the following: creating the basic matrix type, and then hooking up BLAS and LAPACK with it. I believe I used a bits of narray and perhaps even blt::vector as a template, and then built the matrix on top. Had lots of issues with slicing syntax, but mostly because I didn't want to do work too hard. BLAS and LAPACK came trivially after I had the backend done (in C++, not that it really matters). I believe I used cephes

math library for the usual stuff. Wrapping external libraries was somewhat of a pain, but only because automated tools like SWIG weren't readily available back then.

Lots of things were missing of course - good iterator support, generalized slicing, limited math library support for these containers. And there were the inevitable glaring bugs. Wonder I still have it all somewhere in my CVS repo ... not that I want to embarrass myself by looking at it again.

This was before Tcl_Obj came about, so the implementation lacked good structure, and more importantly, good aesthetics.

Bottom line is that I'm going to keep on enjoying embedding Tcl in most of our simulation systems which benefit from it, while using other scripting languages such as Python for prototyping new (numerically intensive) models that I work with on a daily basis. Tcl is not about to lose its advantage just because it doesn't try to the be all end all, and that in my book is a good thing.

[Mumit Khan], in February, 2001, writes:

What type of numerical support would be useful to the community in general?

Good question, with no easy answers. There has been on and off discussion regarding this topic, and these tended to die out right away due to lack of interest (I plead guilty to not contributing). The last significant posting on this topic was by Kevin Kenny, who did a great job of summarizing the issues involved; the silence that followed after one or two followups, yours included, is enough to cause those numerically- sensitive among us to observe a moment of silence (ok, that's way overboard, but you get my point). Kevin points out the fundamentals, and provides a short roadmap of what needs to be done, and I doubt if I can improve on that. Please search for the Numerical Analysis in Tcl

, comp.lang.tcl, 2000-10-02.

 >Perhaps someone could prototype an extension that would demonstrate what
 >is needed?

Kevin's first point is the lack of a good numeric container in Tcl. The Tcl arrays and lists are obviously not that useful, and so that's the first step. BLT's vector is actually partly there, especially with its support for "slices"; however, I'd take out a good part of its guts for performance reasons. Matrix and multi-dimensional arrays built similarly would be a great boost, along with consistent syntax for generalized slices that these extended types would handle internally (hence no changes to core in that regard at least).

Kevin then points out the issue of type shimmering, and it's not just a Tcl issue. Python in a slightly different way needs to handle this as well, eg., shimmering between Python lists, array in bundled library and array in Numeric extension. Not the same, not similar.

Extensions are actually not that hard, but the hard part is getting the core support to make both the performance worthwhile and the syntactical issues usable/tolerable. Tools like SWIG can wrap a reasonably complex numerical library, and as Kevin points out, there are lots of high quality and free ones in wide use, with no big shakes; but of course only after we settle on a reasonably usable basic types.

There was also some discussion on math functions in tcllib, which will also be very useful.

Even with all this, there's always the syntactical issue:

set y [expr {$x + [f [expr {pow($x,$z)}]]}]

y = x + f(x**z)

or even

y = x + f(pow(x,z))

Kristoffer Lawson replied:

: Even with all this, there's always the syntactical issue: : set y [expr{$x + [f [expr {pow($x,$z)}]]}

I'd still rather see mathematical functions used as actual commands, so:

set y [expr {$x + [f [pow $x $z]]}]

or even

set y [add $x [f [pow $x $z]]]

There are some hacks that do it but they use expr in the end so no doubt they're not as speedy as is possible. I've always wondered why expr has in-built mathematical functions anyway? Why are they not all simply Tcl commands? Performance?

Anyhow, I'd definitely approve the above to the one suggested below:

:   y = x + f(x**z)

This is just too un-Tclish.

[Mumit Khan]'s response:

 >This is just too un-Tclish.

Perhaps I didn't make myself clear -- I'm not trying to advocate this syntax for Tcl, just pointing out what a "typical" numerical code may look like in a different language as opposed to in Tcl.

Tom Poindexter jumps in with:

 >look like in a different language as opposed to in Tcl.

But Tcl can handle this syntax in a Tcl fashion -

math {
    y = x + f(x**z)
}

Where 'math' is some not-yet-written extension to parse and execute it's own flavor of expressions. Tcl already has mini-languages already - the printf type specifiers of scan and format, regexp/regsub expressions, even the existing expr. There's no reason why the hypothetical 'math' shouldn't allow expressions more natural to its domain.

[email protected] wrote:

 > I disagree that the _core_ of Tcl needs stronger numerical calculation
 > support.   I would _agree_ however that stronger numerical calculation
 > support is needed.
 >
 > It is my personal opinion that the Tcl core needs to be strong enough
 > to enable people to write 80% or more of any particular application -
 > and stronger so that extensions can be written to fit the other 20% or so.
 >
 > Adding large amounts of code to the core for specialized applications -
 > whether it is numerical analysis, SAP, security, WWW, etc. only adds the
 > possibility of additional fragility and dependencies that could result
 > in a collapse of the core.

Hear, hear!

Let me add another voice calling for stronger numerical support without compromising the integrity of the core. What we need, though is a roadmap for getting there from here.

Let's examine one straw-man:

Virtually all the numeric codes already in existence require single- and multi-dimensional arrays of integers and of single- and double-precision floating point numbers. For the success of numerical calculations in the Tcl core, such arrays must have natural representations as Tcl objects. George Howlett has done some work on Tcl representations of vectors and matrices.
One problem with inter-converting vectors with Tcl lists and strings is type shimmering -- the fact that the Tcl library can represent two types for an object but not three leads to difficult cases where types are converted inside inner loops. For this reason, the Tcl numeric community should consider Feather a high priority for integration in the core.
Once representations of vectors and matrices (in column-major order, please!) are available, efforts can begin on adapting available numeric codes to Tcl language bindings. Much of this work initially could be semi-automated conversion of Fortran calls to Tcl commands. This could be handled by a package like SWIG, extended to support passing array dimensions as explicit parameters. (A few years ago, I made a start at such a thing, directed at wrapping CGI scripts around available Fortran codes. It's still out on the Web at http://camnet.ge.com/ftncgi/ and I'm darned if I know whether it still works.)
There are a number of reasonably high-quality public codes available as candidates for inclusion. Comprehensive libraries such as SLATEC might be one starting point. It would be interesting, too, to adapt some of the more complex libraries to use Tcl conventions. An adapter so that an optimization package like MINPACK could do its function (and Jacobian?) evaluations with Tcl callbacks would be an incredibly powerful framework for setting up all kinds of solutions.
Tacking something like VTK ( http://www.kitware.com/ ) onto the back end would then add a terrific visualization environment to the mix. Alas, VTK would have to be made TEA-compatible; I've been out of touch with the VTK guys, so I don't know if there's any intention to do so; last I heard, they were still satisfied with Tcl/Tk 8.0.2 and saw no reason to move forward.

The only item that requires patching the core is the Feather support, to make it easier for extensions to add new types of first-class objects and avoid shimmering. Once that's done, and we put together Tcl bindings for the BLAS, then we can turn hordes of developers loose with SWIG or something similar to bring in the libraries.

Of course, this is only a straw-man proposal; I eagerly solicit comments!

 --
 73 de ke9tv/2, Kevin KENNY
 GE Corporate R&D, Niskayuna, New York, USA

Donal Fellows replied:

In article <[email protected]>, Kevin Kenny <[email protected]> wrote:

 >   * One problem with inter-converting vectors with Tcl lists and
 >     strings is type shimmering -- the fact that the Tcl library can
 >     represent two types for an object but not three leads to
 >     difficult cases where types are converted inside inner loops.
 >     For this reason, the Tcl numeric community should consider
 >     Feather a high priority for integration in the core.

Tcl doesn't have support for more than one (non-UTF8) type at a time because of a trade-off between speed of execution and memory usage. While it would be possible to extend Tcl_Objs with a second internalRep-like field, this would expand a significant proportion of the structures present in *every Tcl application* by 50% (out from 24 bytes to 36 bytes on a 32-bit architecture.) The code to manage these values would also become significantly more complex; maybe the result would be faster, but the number of bugs present would undoubtedly increase by a lot too. And finally, the gain from doing this wouldn't be too spectacular in the case where you have updates from the PoV of a non-string rep, since then you may well have additional non-string reps to discard in addition to the string rep.

All in all, the current implementation is good enough for 99% of its uses, and the remaining 1% would add a lot more work (as well as being incompatible with lots of existing code in extensions.)

FWIW, I just reckon there ought to be a nice easy way to link Tcl to things like C and FORTRAN for this sort of thing (SWIG and/or Ffidl.) Why try to be all things to all people when we can enlist others to help out instead?

Donal.

Hear hear! It is remarkably easy to wrap external C and FORTRAN libs using SWIG. This is what the Python and Perl communities have been doing for years to quickly add new functionality without modifying the core language.

 From: Kevin Kenny ([email protected])
 Subject: Re: Numerical analysis in Tcl (Was: Re: RFC: Dramatic overhaul
 of lists, eval, and subst)
 Newsgroups: comp.lang.tcl
 Date: 2000-10-03 10:57:29 PST

"Donal K. Fellows" wrote:

 > Tcl doesn't have support for more than one (non-UTF8) type at a time
 > because of a trade-off between speed of execution and memory usage.
 > While it would be possible to extend Tcl_Objs with a second
 > internalRep-like field, this would expand a significant proportion of
 > the structures present in *every Tcl application* by 50% (out from 24
 > bytes to 36 bytes on a 32-bit architecture.)  The code to manage these
 > values would also become significantly more complex; maybe the result
 > would be faster, but the number of bugs present would undoubtedly
 > increase by a lot too.  And finally, the gain from doing this wouldn't
 > be too spectacular in the case where you have updates from the PoV of
 > a non-string rep, since then you may well have additional non-string
 > reps to discard in addition to the string rep.

OK, I gave a bit of a wrong impression here. Feather doesn't multiply internal representations; instead, it lets an internal rep have multiple interfaces associated with it. With Feather, the Tcl_* functions that operate on list objects get delegated to a container interface, which any Tcl_Obj type is free to implement. If an object's internal rep doesn't support the container interface, then a call to Tcl_ListObj* will result in coercing the object to a list, which is the default container. In this way, an object that represents a structure, a vector, or another data type can function as a list and not shimmer.

Paul's paper from the Tcl2K conference in Austin [2] discusses all this.

 --
 73 de ke9tv/2, Kevin KENNY
 GE Corporate R&D, Niskayuna, New York, USA

 From: Donal K. Fellows ([email protected])
 Subject: Numerical analysis in Tcl (Was: Re: RFC: Dramatic overhaul of
 lists, eval, and subst)
 Newsgroups: comp.lang.tcl
 Date: 2000-10-04 14:10:04 PST

In article <[email protected]>, Tom Krehbiel <[email protected]> writes

 >2) A second and even more important consequence of not having a clear
 >vision of what "tcl" really means (i.e. core plus what??) is that the
 >development community doesn't know how to proceed. The development
 >of many reasonable extensions requires the existence of other  extensions
 >(i.e. packages are coupled) without guidance about what direction the
 >language is going the result at best is confusion and wasted time.  This
 >has been particularly apparent with regard to the state of OO in tcl.

There are some extremely different ideas about what OO in Tcl should look like. About the only thing you'll get reasonably rapid agreement on is that methods should be called something like this:

  $object $methodname arg1 arg2 ...

After that, you're getting onto contentious ground... :^(

 >3) I think that tcl is in desperate need of a good numerical package but
 >a developer would be a fool to create one unless the OO issue is  resolved.
 >A package like NumPy is by its very nature OO because the entities
 >that are being manipulated are not simple data types.

There's so many different groups wanting extended numerical/mathematical packages, and they all want something different. HPC people want native complex numbers, native matrices and tensors, and interoperability with FORTRAN. Financial people want exact control over precision and rounding, but aren't too bothered about trigonometric functions[*]. Security hackers want fixed-field arithmetic. And so on.

I don't think that supporting all of these in the core (or even at the SUMO level) is very practical. This sort of code is often hard to get right and there aren't that many people working on the core. It is much better to provide a generic interface and let people put their own functionality in there behind it. But that is exactly what Tcl is all about; not providing advanced mathematical and numerical handling is precisely within the spirit of the language. Not that we can't discuss what the common interface would look like. :^)

For the record, I have no need for this sort of stuff in my own code. Networking and Tk are much more useful to me, and regexps are a good thing too.

 >If you interprete the above to mean that I don't like tcl, you would be
 >wrong, I just want to make it better.

I don't interpret it that way at all. I also want to make Tcl better, but I don't want to lose what makes Tcl good now either. It's all a matter of balancing Gain against the Gain/Pain ratio...

Donal.

[* Banking is complex enough without square root of minus one! ]

przemek writes:

[email protected] (Mumit Khan) writes on embedding Octave in Tcl:

 > Before even considering something like Octave, the first order of business
 > is creating a common set of "primitives" -- vector/matrix with usable
 > syntax for slices. FYI, Octave is rather hard to embed in other
 > applications; the liboctave and libinterp libraries are a start, but it
 > was never really designed to be easily embeddable,

Yup, reusing liboctave in Tcl looks scary, plus currently octave does 2D matrices rather than full N-dim; the proper job has to be a complete redesign. I notice that even the matter of proper Tcl syntax for math has been hotly debated in this thread, not to mention the design of internal data structures, etc.

I have a modest proposal for the meantime: leave everything math-related to octave (syntax/parsing of expressions, data storage, even I/O), and drive it like an external process. Octave is pretty good in this stuff (nice math-friendly syntax a la [s, v, d] = svd(a), high-performance I/O using HDF libraries, its own autoloading, etc.), so it's smart to let it use them.

Of course we need a way to exchange data efficiently: so there's a need for some IPC. At the simplest, it could be through ASCII lists:

set ta  { 1 2 3 4 5}
set tr [octave  "a = $ta; vander(a)"]
# tr = list of lists of rows of the Vandermonde([1 2 3 4 5]) matrix

Of course for any significant numerical processing there would have to be some fast 'unformatted' method of passing data; possibly via shared memory, unless there is a more portable method. The sanest way would probably be directly from Octave's array variables to BLT vectors. This is really not a very general design, but for most practical uses Tcl would need the data for plotting and such, where the BLT vectors do nicely:

vector va vr
va set { 1 2 3 4 5}
octave  "blt([vechandle va], a);
         z=vander(a);
         blt(z(3,:),         [vecHandle vr])"

where vecHandle is the Tcl procedure that returns a vector handle that octave can use (presumably something like a shared memory ID), and blt() is an octave extension that knows enough about BLT vectors to be able to transform data between them and octave arrays.

Note that this would work well doing the opposite, i.e. for using Tcl from within octave.

Arjen Markus: I will not comment (too much) on the previous text (though I know a few other faults in C's design that continue to cause great grief - arrays for instance, null-terminated strings and the impossibility to check pointers and I do know at least two other programming languages that use a different assignment operator than a single equal sign - Ada and Eiffel, possibly Algol did/does as well).

My contribution is the following:

A suitable data structure that does not require adjustments to the core is one that employs binary strings as an opaque data type that is very suitable for transferring contiguous data to and from extensions written in C or Fortran.
The documentation on the C interface for BLAS sheds some light on the issue of how to interface with Fortran in the presence of large arrays and computationally intensive algorithms. http://www.netlib.org/ is good starting point (I forgot the rest).

To illustrate the internal representation of a two-dimensional array that can be passed to C or Fortran:

set data_array [list 10 20 [binary format f* $list_of_floats]]

The third element of this list (for Tcl and C programmers, number 2) is the array that maps directly into a C float[x][y] array or a FORTRAN REAL, DIMENSION(x,y) array, where x and y are 10 and 20 or vice versa.

JCW: That is precisely what VKIT supports: access to opaque binary data, with efficient indexing for different formats. Vectors can have arbitrary attributes, so adding a "dimension" with value "10 20" would be easy. The descriptor-based design means that the data is described, not forced into a straight-jacket. Furthermore, VKIT's design gives you the option to leave all/some data on disk and use memory-mapped files instead of loading it into a buffer.

Binding CBLAS into Tcl would be a perfect little task for Critcl. In fact, so would binding to the "netCDF" datafile exchange standard lib. I'd be interested to do all this - if netCDF is indeed popular (is it?). Please let me know if you want this - I'm keen to tie into numerical data handling, but don't really know what people use these days.

Well... CBLAS now exists, and the Tcl + Critcl combo made it fun, see [3] and the generated C code [4]. The cblas.h header is parsed into a Tcl script, and the cblas.tcl package uses this script to generate some 140 C bindings on the fly. There's a demo at the end, and the test output shows that it all works. Warning: the current data/array handling approach has plenty of rope to hang yourself. -- JCW

Latest news, Nov 23 2001: a wrapper for CLAPACK also appears to be within easy reach. It's not there yet (arg handling, in/out arrays, ptrs), but as the output below indicates, it's pretty much within reach now - JCW

This script:

package require clapack
namespace eval clapack {
    catch cbdsqr err
    set v {}
    foreach x [info commands ::clapack::*] { lappend v [namespace tail $x] }
    puts [lsort $v]
    puts "[llength $v] commands defined in ::clapack:: namespace"
    puts $err
}

Generates this output:

 cbdsqr cgbbrd cgbcon cgbequ cgbrfs cgbsv cgbsvx cgbtf2 cgbtrf cgbtrs cgebak
 [...]
 zunmtr zupgtr zupmtr
 1217 commands defined in ::clapack:: namespace
 wrong # args: should be "cbdsqr uplo n ncvt nru ncc d e vt ldvt u ldu c ldc rwork info"

Arjen Markus: Changes I made a couple of days ago seem to have gone lost (I did have trouble saving the text). The point I do want to make here, is that the Fortran 90 interface to BLAS seems to be more natural to Tcl than the Fortran 77 or C interface: Tcl will already have information on the array dimensions, so these arguments to the underlying routines should not be visible to the Tcl user.

JCW: Agreed 100%. IMO that should be done at the Tcl level, fitting in with comments I made elsewhere, see a critical mindset about policy.

AM (28 december 2009) I have made a start with a wrapper extension for LAPACK. As there are a number of design choices to be made, I welcome discussion.

gold 2014-05: Part of the reason that early Tcl was so successful was that really good & clean working programs were furnished as examples with the Tcl distribution. For myself, I'd rather see statements from a math library in a working Tcl program than tons of algebraic definitions. Maybe there are better sources or distributions of working examples for the math libraries. If you are planning for more numerical analysis libraries, companion example programs should be either part of the library effort or the work of examples partially farmed out to supporters.

Page Authors edit

AM
Arjen Markus
Lars H

arjen - 2014-10-13 08:36:40

I have completed a first version of a wrapper for the Cephes library - see the Chiselapp site

With this wrapper you can access, say, the Bessel function of the first kind of order 1, J1, as:

    package require tcl_cephes
    puts [::Cephes::J1 $x]

Category Mathematics