[SciPy-dev] fblas tests

[pearu at scipy.org]

Hi,

On Mon, 15 Apr 2002, eric wrote:

> 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.

So you prefer living on an edge ;)

Say, for a function foo we can have two signatures:

1)  x = foo(x,overwrite_x=1)
PLUS: it is convinient for power user not to specify overwrite_a if
in-situ change of x is allowed.
MINUS: Explicit is better than implicit. The above is implicit.
MINUS: So, if a careless/tired/etc power user defines

  def bar(x):
    x1 = foo(x)
    return x1 + x

then the result can be unpredictable: if foo increases its argument with
type 'i' by one, then we have
  bar(1) -> 3
  bar(array(1,"i")) -> 4

2) x = foo(x,overwrite_x=0)
MINUS: if in-situ change of x is desired, then user must specify
overwrite_x = 1 argument.
PLUS: if in-situ change of x is desired, then user must specify
overwrite_x = 1 argument. Explicit is better than implicit.
PLUS: Careless/tired/etc power user can safely define bar as above.

Purely because 'Explicit is better than implicit' I would prefer
default overwrite to be 0, even for low level functions such as in
fblas.
It is more typing for the developer, but later when looking the code, it
is clear what developer expected to happen. In-situ changes are not
Pythonic and therefore such things should show up explictely when used in 
Python.
Even a C developer while using Python would not assume in-situ changes
by default just because it is not Pythonic, IMHO.

Personally, I am not against intent(overwrite). But I just imagine that
intent(copy) is safer for other developers. Either they use overwrite=1
if it is really needed or they don't need to bother about setting
overwrite=0 if overwrite is not expected or desired.

And note that using overwrite_x = 1 makes sense only if you don't care
what would be the contents of x after the call. Such a situation can be
useful only for temporary variables with initially filled with input
data. You should _not_ use overwrite_x = 1 to get the same effect as if
the argument would be intent(inout). Arguments that are changed in-situ
and filled with output data, should always be defined as intent(in,out),
that is,
  x = foo(x)
is better than
  foo(x)
if x is going to be changed in-situ.
Do you agree?

> 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.

OK, let's assume only power users in what follows.
But I think even power users should write clear and explicit codes and
with a style proper to a given language.

>  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.  

I disagree. The whole point of wrapping Fortran to Python is to make
coding easier. This includes that some technical details that can be
established from arguments in one-to-one fashion, can be hidden.
This is possible for Python types but not for Fortran or C types (think of
an array and its shape). So, coding is easier because: 1) less typing,
2) less possibility to make errors, 3) less bugs as there are some
internal checks.

So, I don't see why would one want to call a Fortran function from Python
if it would have, say, 20 arguments and consider all technical details
while the work would be the same as using Fortran directly. I would not
write Python code in this case, I would write it directly in Fortran then.

> Otherwise using pure Python will gets you reasonably close to the same
> speed.

Not true. If wrappers are used properly. The same holds for wrappers that
expose everything. Only with proper use one can have speed ups.
In the former case it takes less efford to learn proper use than in the
latter case.

> > 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.

I see your point but I find it impractical.
I don't see a point of experimenting with these low-level routines
from Python if I am later going to translate the algorithm to Fortran.
I would write the algorithm directly in Fortran if it is my goal to have
it in Fortran. Otherwise the work would be doubled. Note that developing
an algorithm in Fortran is as efficient as in Python if using tools like
f2py: all testing and I/O is done in Python while algorithm is exposed
to Python only with one command:
  f2py -c ...

I can imagine that such a developing style is not classical and many
developers hardly use it but it is acctually extremely efficient, both
the development and the resulting software.

> 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.

Exposing array shapes to Python has a minus because it increases the
wrappers code, which already is quite large. Even more, if it is not clear
whether they are useful or not.

I would suggest leaving them out and if it turns out that such arguments 
would be useful to have, then add them later. It is better to extend the
set of features rather than having useless features lying around.

>  I think these have been replaced in some sense by offx and offy which
> aren't even a part of the Fortran interface.

Actually, such things like calling 
  foo(x(10))
where x(10) refers to a sub-array starting at index 10 (and not only the
10th element), is part of Fortran language and widely used in Fortran
codes. The equivalent C codelet would be
  foo(x+10)
and Python one is
  foo(x,offx=10)

I don't know how else can you model Fortran statement `x(10)' other than
using offx keyword argument. Using
  foo(x[10:])
from Python would be quite inefficient compared to foo(x,offx=10): just 
think what is happening behind the scenes in those two cases.

>  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.

Well, see above, I have tried my best in assuring you. You seem to think
that I don't care about performance or generality. They are very much my
first priorities. Besides that I work hard to get the simplest wrapper
possible, even if it takes several iterations. May be that gives an
impression that they cannot be efficient or general. But I would hate to
produce complicated wrappers just that they would look sophisticated ;)

> >  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.

Yes, they are. See the foo(x,offx=10) example above.

> 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.

Yes, please do.

Pearu