[Numpy-discussion] How to Force Storage Order

Sebastian Berg sebastian at sipsolutions.net
Tue Mar 31 03:36:47 EDT 2015


On Di, 2015-03-31 at 07:11 +0000, Ian Henriksen wrote:
> On Tue, Mar 31, 2015, 12:50 AM Klemm, Michael
> <michael.klemm at intel.com> wrote:
> 
> 
>         Dear all,
>         
>         I have found a bug in one of my codes and the way it passes a
>         Numpy matrix to MKL's dgemm routine.  Up to now I was assuming
>         that the matrixes are using C order.  I guess I have to
>         correct this assumption :-).
>         
>         I have found that the numpy.linalg.svd algorithm creates the
>         resulting U, sigma, and V matrixes with Fortran storage.  Is
>         there any way to force these kind of algorithms to not change
>         the storage order?  That would make passing the matrixes to
>         the native dgemm operation much easier.
>         
>         Cheers,
>                 -michael
>         
>         Dr.-Ing. Michael Klemm
>         Senior Application Engineer
>         Software and Services Group
>         Developer Relations Division
>         Phone   +49 89 9914 2340
>         Cell    +49 174 2417583
>         
>         
>         Intel GmbH
>         Dornacher Strasse 1
>         85622 Feldkirchen/Muenchen, Deutschland
>         Sitz der Gesellschaft: Feldkirchen bei Muenchen
>         Geschaeftsfuehrer: Christian Lamprechter, Hannes Schwaderer,
>         Douglas Lusk
>         Registergericht: Muenchen HRB 47456
>         Ust.-IdNr./VAT Registration No.: DE129385895
>         Citibank Frankfurt a.M. (BLZ 502 109 00) 600119052
>         
>         _______________________________________________
>         NumPy-Discussion mailing list
>         NumPy-Discussion at scipy.org
>         http://mail.scipy.org/mailman/listinfo/numpy-discussion
> 
> 
> 
> Why not just call the algorithm on the transpose of the original
> array? That will transpose and reverse the order of the SVD, but
> taking the transposes once the algorithm is finished will ensure they
> are C ordered. You could also use np.ascontiguousarray on the output
> arrays, though that results in unnecessary copies that change the
> memory layout.
> Best of luck!
> 

Frankly, I would suggest to call some function like that in any case (or
at least assert the memory order or have a test suit), just to make sure
that if some implementation detail changes you are not relying on "this
function will give me back a fortran ordered array".
Of course you can do tricks to make this extra function call a no-op,
since it will do nothing if the array is already in the correct memory
order. But the quick check that everything is in order should not matter
speed wise. To be honest, after doing an SVD, even the copy likely
doesn't matter speed wise....

- Sebastian


> -Ian Henriksen 
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 819 bytes
Desc: This is a digitally signed message part
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20150331/2c7e3079/attachment.sig>


More information about the NumPy-Discussion mailing list