[SciPy-dev] lapack and blas

[Pearu Peterson]

On Fri, 18 Jan 2002, eric wrote:

> Hey Pearu,
> 
> > Now I would like to re-work with lapack and blas modules in
> > scipy. Currently they are included in linalg package but I am thinking of
> > separating lapack and blas from there to stand-alone packages into scipy/
> > tree. Here are my reasons:
> > 1) I don't want to break linalg package.
> > 2) lapack and blas packages are valuable also outside scipy environment.
> 
> I think some reworking would be beneficial.  Can you give a brief outline of
> your plans?  There is hope that the current wrapper generator will also be
> able to generate Scalapack wrappers from the same set of generic_lapack.pyf
> definitions.  I want to keep this possibility open.  We also need to keep
> both c and fortran interfaces so that people using the respective languages
> will have a "native" ordering package to use for the fastest possible
> computations.

I was thinking of creating two extension modules: flapack.so and
clapack.so. One would be efficient for column-major ordered arrays, the
other for row-major ordered arrays, respectively. Both should have
identical interfaces for Python: if one is not available then other should
serve the same needs, just with less efficency due to ordering issues.
Now, flapack.so and clapack.so are used from a Python module lapack.py
that is interface for lapack users. Basically, lapack.py
implements the selection mechanism for choosing Lapack functions from
flapack.so or clapack.so based on the array types and ordering (btw, f2py
provides an efficent way how to decide if the array is in column-major
orderd).
Also, if one of the extension modules is not available (as it may happen 
with flapack.so on Windows), then it would automatically use clapack.so
instead. (Btw, clapack has also functions for column-major arrays, we
could also consider supporting those functions).

The same scheme applies also for fblas and cblas. If scalapack has also a 
C version available, then we can wrap it in the similar way.

Note that I do see advantages in the current wrapper generator, I was
thinking that some of its code should re-used. Basically the
<c,s,d,z>-stuff and may be info-stuff also. Many things like reordering
axes are dealt by f2py generated modules.

In order to deal with C versions of lapack subroutines, that are actually
C functions, I will introduce intent(nowrap) attribute for function names
so that f2py will not generate wrappers for them. In this way we get rid
of global --no-wrap-functions flag. (Or if subroutine has already
intent(c) attribute, then f2py should not create a wrapper for it?)

> Also, I have a small set of routines that overlap with some in linalg that
> we use on some high performance (but serial) applications.  They're
> interface is necessarily differnt from the current linalg stuff so that
> uncessary copying, etc. can be avoided (they have and overwrite=1 flag and
> things like that).  I'd like us to revisit the linalg interfaces in the
> hopes that we can come up with a design that provides convenience and still
> allows for high performance computations.

This issue is also dealt by f2py. There is now intent(copy) attribute for
array variables. By default, a copy of the array is made before passing it
(acctually its copy) to the wrapped function. f2py generates also keyword
arguments for such cases that can be used to disable the copying for
efficency. For example, consider

subroutine foo (a,..)
   integer intent(in,out,copy):: a(m,n) 
   ! do nothing
end

Its Python signature is

   a = foo(a,copy_a=1)

If using the proper type arrays for arguments then the following holds:

   id(foo(a))  != id(a)         # a copy of `a' is passed to foo and
                                # returned to Python

   id(foo(a,copy_a=0)) == id(a) # no copy is made, the same array is
                                # return to Python

So, I hope this solves the problems of efficient memory usage in a
flexible way.

> > Later, linalg would start to use these new lapack and blas packages as
> > they should become easier to maintain and more efficient when using the
> > latest features in f2py.
> 
> Is it possible to just branch the current package in the CVS? I've never
> done this (the one time I tried, I screwed everything up).  We can add an
> alias in the modules file so that cvs co linalg just grabs the linear
> algebra package.  Also, the current build infrastructure makes it possible
> build sub-packages in scipy either within SciPy or standalone.  So, if there
> are no scipy dependencies in the linalg module, then it should work fine
> elsewhere.

I haven't tried to branch CVS repositories either.

> If we can find a way to do it, I'd like to do it within the current package.
> That keeps from having two of them floating around.  Also, it will keep
> everybod interested in the re-write involved in it.

It is fine for me to do it in inside linalg as scipy already has so many
directories with different stuff. However, I would like to work without
breaking current linalg until we have wrappers to all these functions that
linalg currently uses. It shouldn't take too long but it would be safer
for now. I would suggest separate sub-directories for
lapack,blas,scalapack,etc in linalg/.

Also, I would like keep the set of wrapped functions minimal as a first
instance. We should choose some representatives and test the re-write on
them. After some experience we can wrap other functions as well, but then
with less pain.

> I'm all for the re-write though, and am currently interested in the same
> topic.  Lets get to it.

Very good, lets do it.
	Pearu