[SciPy-dev] New clapack and flapack

[eric]

Hey Pearu,

Jeeze.  You work fast.  I haven't even had a chance to grab the new f2py
yet...  I'm guessing I need to update f2py from CVS to try all this out.
I'll try and spend some time with this today, but it may be tomorrow before
I have a chance to take a look.

> 1.
> $ python interface_gen.py   # this creates clapack.pyf and flapack.pyf
>
> 2.
> $ f2py -c
lapack.pyf -L/usr/local/lib/atlas -llapack -lf77blas -lcblas  -latlas -lg2c
>
> $ f2py -c
flapack.pyf  -L/usr/local/lib/atlas -llapack -lf77blas -lcblas -latlas -lg2c
>
> This builds clapack.so and flapack.so.

Looks easy enough.  f2py is really gotten things pretty automated.

> 3.
>   LU,piv,X,info = sgesv(A,B,rowmajor=1,copy_A=1,copy_B=1)
>   LU,piv,info = sgetrf(A,rowmajor=1,copy_A=1)

A few points about naming.

1.
SciPy has settled on all lower case convention for consistency across
interfaces.  So using:

    lu,piv,x,info = sgesv(a,b,rowmajor=1,copy_a=1,copy_b=1)

will fit the scheme.

2.
Second, is "copy_A" a f2py generated name or can we specify it?  I prefer
overwrite = 0 to copy=1 because it is more explicit about what is happening.
It really warns people that they are clobbering data when they choose to set
the flag.  Also, the "A" part of copy_A will be ambiguous during usage.  In
the worst case for the sgesv example, imagine someone does the following:

   >>> a = rhs_vector
   >>> b = lhs_matrix
   >>> lu,piv,x,info = sgesv(b,a,copy_a=1,copy_b=0)

If the person is trying to save their variable "a" here, they didn't do it,
because the "copy_a" is refering to the function's argument "a", not the
variable bound to it which is "b" in this example.  I ran into this
conundrum once when writing an interface generator in the past.

It be better if the names somehow were more descriptive.  For sgetrf, the
following seems fine:

    lu,piv,info = sgetrf(a,overwrite=0,rowmajor=1)

For sgesv, the naming isn't quite as obvious.  Maybe,

   lu,piv,x,info = sgesv(a,b,overwrite_rhs=0,overwrite_lhs=0)

3.
In the C API's, I think the overwrite argument will be used much more than
the rowmajor flag.  It should come first in the argument list.  In fact, in
the "old" wrappers, rowmajor was just hardcoded to 1 because the flalpack
routines were considered the colummajor functions.  Keeping it is fine, it
just should go to the end.

    >>> lu,piv,info = sgetrf(a, overwrite=0, rowmajor=1)

Not sure how hard these changes are -- I expect they might require some
enhancements to f2py.  I'm not opposed to having Python wrappers around the
routines that help with this stuff.  We used to have to do that for the
clapack wrappers (until your additions to f2py's capabilities made this
unnecessary).  If it is easy to get f2py to do this stuff, then, of course,
that is better.

Comment on testing.
There was a fairly complete set of tests written for the *interfaces* of the
blas routines.  These guys only tested how the interfaces worked and didn't
do much testing of the functions accuracy.  We need both.  I don't remember
if the LAPACK tests were as complete -- I certainly didn't write them, but
Travis O. may have.

Based on the errors we saw a week or two ago in the CLAPACK functions, we
really do need to get some tests implemented for this stuff.  You and I
probably handle writing the interface tests (very tedious, but maybe not as
bad with all your changes to f2py), but we really should find an expert to
put the linear algebra routines through their paces.  We can write some
"smoke tests" though that make sure LU and solvers, etc. work for general
problems.  I'd stick this way up there as one of the most important packages
in SciPy, so we should really try and make sure it works right!

Maybe there is a test set of matrices laying around somewhere to help us out
in all this.  Anyone know of such a thing?  -- perhaps in the LAPACK
distribution itself.

> Few remarks:
> 1) The sources to these extension modules are quite big (more than 5000
> lines in both). Should we factor these modules to
>    cslapack.so cclapack.so cdlapack.so czlapack.so
>    fslapack.so fclapack.so fdlapack.so fzlapack.so
> (any better name convention?). What do you think?

Well, I think in practice, they'll rarely be needed separately.  Having big
wrappers isn't that big a deal (heck, I've had SWIG wrappers that were 60000
lines...).  I'd vote for keeping them all together unless there is a
compelling reason besides size to split them.  This keeps things a little
simpler (import once, get everything).

> 2) How shall I commit all this to the scipy CVS repository? In a
> subdirectory of linalg? Eventually it should be in linalg directory, I
> think.

We'll since we don't know how to do CVS branches reliably... :| Why don't
you make a directory called linalg2.  We'll work in there until it gets
stable.

Is there a CVS expert out there that can help us out on this?

Pearu, this all looking very promising.  Thanks a ton.

eric