[SciPy-dev] fblas tests

eric eric at scipy.org
Mon Apr 15 08:33:37 CDT 2002

> > I was looking over the behavior of the new fblas wrappers and think we
> > make some changes back to how the previous version worked (for reasons
> > below).  The old test suite is pretty extensive and covers many of the
> > calling approaches that can occur.  We should definitely try to get all of
> > tests to pass, as they exercise the interface fairly well.
> I hope you don't want the old interface back just because to get the old
> tests working quickly again. In this particular case, the interface tests
> should be changed, not the interfaces, I think.

The old tests were pretty much correct, though I did find one error in a strided
array test.  I wasn't aware of the added overwrite keyword, adding it would have
corrected all the old tests.  Still, I feel these guys are explicitly low level
enough that overwrite should be the default behaviour.  If we want to put a
higher level wrapper around them (such as matrix_multipy, or matmult), then that
could default (or even only allow) copy.

> And the main reason is that the new interfaces are definitely better
> (and not just because I wrote them;-).
> The new interfaces are better because
>  a) ... they avoid side-effects (making the intent(inout) feature
>     obsolete).
>     In C or Fortran, intent(inout) type of side-effects are widely used
>     and often the only way to get job done. But this is not the case
>     for Python. I think that we should adopt the interfaces to Python
>     style as much as possible as they are mostly used by Python users,
>     not C or Fortran programmers.
>  b) ... unless power users request side-effects explicitely (the
>     overwrite_* = 1 feature) for performance reasons.
>     Note that the new interfaces do not sacrifice performance or memory
>     consumption if they are used properly, in fact, they are more optimal
>     compared to the old interfaces. Of cource, the proper usage must be
>     documented somewhere.

I think "power users", or at least people that understand explicitly what they
are doing, are the only ones who will ever touch these things.  They have goofy
names, and were designed for use in Fortran to copy data from one array to
another.  Preserving this behavior is pretty much the only way to get any
useful, efficient work out of them.  Otherwise using pure Python will gets you
reasonably close to the same speed.

I've only looked at 5 of the functions so far, but the interfaces looked pretty
much the same on all of them with the exception of the default "copy" behavior
and 'a' not having a default value of 1. for ?axpy (I added it back in).  The
others are different and improved, then that is great.

>  c) ... it is safe and simple for casual users to use them and getting
>     some speed-up compared to not using them at all (scopy is a rather
>     extreme counter-example that shows interface behaving poorly with
>     default values but that does not prove that all other interfaces are
>     poor because these default values are optimal for these other cases.

I guess I feel there shouldn't be a "casual" user of these routines.  They
should never have to know that gemm is a matrix multiply.  We should have
matmult or something like that that they call.

> The old interface could not do all that because f2py did not had required
> features earlier.


> In CVS log you also mention reverting gemv to the previous interface
> style. What do you mean? The task of gemv is
>   y <- alpha*op(a)*x + beta*y
> gemv old signature:
>   gemv(trans,a,x,y,[m,n,alpha,lda,incx,beta,incy])
> gemv current signature:
>   y = gemv(alpha,a,x,[beta,y,offx,incx,offy,incy,trans,overwrite_y])
> Notes:
>  1) in old gemv, arguments m,n,lda are redudant. Why would you want to
>     restore them?

The original purpose of these wrappers was to experts writing efficient linear
algebra routines, such as LAPACK, a tool for experimenting with the low-level
BLAS routines from Python.  I occasionally use BLAS, but never for such
purposes, so I don't know the needs of these users.  As a result, the original
wrappers tried to follow the API as faithfully as possible, but adding keyword
variables when they were helpful.

It looks like m, lda, and n are now hard coded to be the shape of a.  I am not
sure that linear algebra experts always want this behavior.  That is why they
were left as they were.  I think these have been replaced in some sense by offx
and offy which aren't even a part of the Fortran interface.  As a linear algebra
non-expert, I wasn't confident enough to go through and change a long standing
interface.  The current one may indeed by superior to the original, but I think
we need to ask some longtime and expert users of BLAS if our interfaces limit
the capabilities of the original functions in some way.  Pearu, you may be an
expert in the field.  If so, re-assure me that the interfaces are not limiting,
and we'll move on.  Otherwise, we should probably ask someone like Clint Whaley
(if he is willing) to look them over and make comments before we change them.

>  2) trans in old gemv expects a character from string 'ntc' while
>     trans in new gemv expects an int from a list [0,1,2].
>     I see little difference in using either of them, in the latter case
>     the interface code is a bit simpler.

This is fine.  This is not the kind of change I objected to.

>  3) In new gemv, alpha is required. I agree that it should be made
>     optional with the default value 1.


>  4) New gemv has additional offx and offy arguments so that using gemv
>     from Python covers all possible senarious that one would have when
>     using gemv directly from C or Fortran.

OK.  If this is true, then maybe they are a good addition.

>  5) When changing the style here, you should do it also in all
>     other wrappers of blas1,blas2,blas3,lapack routines. I can ensure you,
>     it is not a simple replace-str1-with-str2 job.
> To sum up, I would suggest the following signature for gemv:
>   y = gemv(a,x,[alpha,beta,y,offx,incx,offy,incy,trans,overwrite_y])
> that in the simplest case, y = gemv(a,x), corresponds to
> matrix multiplication of a matrix `a' with a vector `x'.
> And gradually using other optional arguments the task of gemv is extended
> to more specific cases.
> What do you think?

I like overwrite as the default the blas stuff.  And the offx and offy still
need some discussion.  If there are any other heavy duty blas users, please
speak up now.


More information about the Scipy-dev mailing list