[Numpy-discussion] Re: linalg.svd returns fortran arrays, causes problems

[Travis Oliphant]

Zachary Pincus wrote:
>> Can you give us some small, self-contained code that demonstrates the 
>> problem?
>
> u, s, vt = numpy.linalg.svd([[1,2],[3,4]])
> u.flags.fortran
>   True
>
> The problem, as per my original email, is that linalg.svd seems to be 
> now returning fortran-style arrays. I believe that it did not do this 
> in the past, but I am not certain.
I think it's always done this.  It is really very easy to generate 
fortran-order arrays (always has been --- the transpose operation has 
always returned a fortran-order array). 


a = numpy.rand(10,20).transpose()
a.flags.fortran


I think the issue that has raised this discussion is that the reshape 
function is now not running asarray if it doesn't have to --- this was 
done to preserve matrices...

Notice that currently the asarray function will convert Fortran-order 
arrays to C-order arrays (and thus make a copy).  Now that I'm aware 
it's doing this, I actually think that's a bug.   This is because 
Fortran order is really just a special case of striding, and other 
non-uniformly strided arrays are not automatically converted just 
because asarray is called.

So, I think as a side-issue that bug should be fixed.  Now, the real 
issue that brought this to the fore-front is what is meant by the 
"reshape" function.

It has been the case for several months that a fortran-order array has a 
different scan order for it's reshape *method* than any-other strided 
array.   The reshape method is one of those functions that you can't 
avoid thinking about layout of the array (so my previous statement about 
not concerning yourself with striding was too far-reaching). 

We are supporting more than one view about how layout should be done.  
But, we should be doing it in a backward-compatible fashion so that the 
standard is C-contiguous unless otherwise specified.  Thus, I think the 
bug is in the design of the reshape method in that it uses a different 
layout standard depending on whether or not the array is Fortran or 
not.  This bug has been there for awhile, I think, but only exposed by 
the re-working of the reshape function.

I propose to change the reshape method to always expect a C-contiguous 
perspective unless a fortran flag is set (defaults to False), in which 
case the perspective is Fortran-contiguous.  This is true regardless of 
the layout of the current array.

This will fix the reshape function as well but also make it possible for 
someone to use reshape who has a Fortran-order array-view.

>
>
> In any case, there are two issues at play:
> (1) Should any numpy functions be allowed to return fortran arrays 
> unless specifically requested?
> (2) Should all structural operations like reshape (etc.) behave the 
> same for fortran and c-style arrays?
>

The issue is that there are only a few shape-changing functions (ravel, 
reshape, etc.), that the Python user needs to have any perspective at 
all on how the data is layed out in memory.   Other routines should not 
depend on the striding (or if they do they need to ensure it's 
contiguous first). 

For those routines that do require memory-layout knowledge the default 
must continue to be C-contiguous, but we need to add a fortran keyword 
to allow that perspective as well should the user choose it.

This is already done for the ravel method.  And it should be available 
for reshape (and resize---though that one will take more effort and I 
will not do tonight---as well). 


Thanks for the bug-reports.

-Travis