| On Fortran linear algebra "codes" |
[Feb. 12th, 2009|11:25 pm] |
There are a lot of solid, useful libraries out there for scientific computing such as ARPACK, CHOLMOD, LAPACK, and others. These all seem old (with some pages last modified over a decade ago) and all seem to work reliably, but all have lousy APIs with cryptic documentation lacking any examples; none of these Just Work.
For example, using ARPACK to compute the n smallest eigenvalue–eigenvector pairs of a large sparse matrix requires calls to two functions from ARPACK as well as a call to your own linear solver. (LU decomposition so the repeated x = A\b solves are efficient is left as an exercise.) First, you call DSAUPD in a loop.
What is DSAUPD? I think it's something like Double-Precision Symmetric Arnoldi U? Positive Definite. This function takes sixteen (24) arguments, all by non-const pointer. Some of these are magic numbers (as opposed to enums), and some of these are magic strings(!) I've seen this in other Fortran APIs; what are people thinking passing strings around? Sure for solving large-scale eigenproblems, the overhead is small, but it's just the Wrong Way.
It returns a magic number and you the calling code are supposed to perform different operations depending on what is returned until it says it is done.
As if that's not enough, you then have to call DSEUPD with 22 arguments(!!!); that actually gives you the resulting eigenvalue–eigenvector pairs you were looking for.
It's frusterating that this good code is out there tangled behind a lousy, archaic API. Boost.uBLAS provides a decent general-purpose generic C++ classes for linear algebra, but for some reason, Boost (which usually has libraries that Just Work) appears to provide no built-in solvers or wrappers to these old Fortran solvers. Am I missing something or are we really missing a good API to these tremendously useful functions?
The bright light in the midst of all of this appears to be SciPy which provides Python bindings to a lot of this stuff, as well as providing Python bindings for some nice nonlinear optimization libraries (which again wind up with kludgy Fortran calls at their core, but hide those gory details). Still, it seems like (short of rewriting these libraries in a more expressive language than Fortran) these things should be wrapped in a clean, efficient C++ interface and then wrapped in a nice friendly interpreted language like Python (or an ugly kludgy interpreted language like Matlab).
PS
Why must Fortran people call everything a "code"? 0xDEADBEEF is a code; ROT13 is a code (although I guess pedants would beg to differ); I'd even call MOV and NOP codes. A source file contains code and a library is a library, not a code. Come on people, talk like software engineers. </rant> |
|
|
| Reading hideously-long argument lists |
[Feb. 15th, 2007|12:19 pm] |
Dear Lazyweb,
In a C++ function with C linkage, are the arguments stored in a predictable way such that, e.g., I could cast the location of the first argument to a pointer to a struct containing the argument list?
Background
There is a commercial finite element package, ABAQUS, which has the nice feature of letting you programmatically define arbitrary materials. You just provide it with an object file containing the function umat_.
Unfortunately, since this interface was designed by Fortran programmers a while ago, umat_ takes 38 arguments, many of which are two-dimensional arrays. Since it is the year 2007 and some languages actually provide meaningful abstraction, I quickly made a umat_arguments class so my umat_ function just calls a C++ function, int umat(umat_arguments& in). Of course, I have to construct the umat_arguments structure which means copy-paste coding – exactly what I'm trying to avoid.
While the huge argument list is dirty, it is nice that ABAQUS has presumably put a lot more thought than I want to into exactly what a finite-element package would need to ask of a material point to be completely general. |
|
|
![[info]](http://l-stat.livejournal.com/img/userinfo.gif){kind=small}
benfrantzdale's journal