How to Force Storage Order
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
On Tue, Mar 31, 2015, 12:50 AM Klemm, Michael <michael.klemm@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@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! -Ian Henriksen
On Di, 2015-03-31 at 07:11 +0000, Ian Henriksen wrote:
On Tue, Mar 31, 2015, 12:50 AM Klemm, Michael <michael.klemm@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@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@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Dear all,
-----Original Message----- From: numpy-discussion-bounces@scipy.org [mailto:numpy-discussion- bounces@scipy.org] On Behalf Of Sebastian Berg Sent: Tuesday, March 31, 2015 9:37 AM To: numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] How to Force Storage Order
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....
After Ian response I have decided to go for properly invoking dgemm with the right transpose parameters if the data is in transposed format. That's way more generic and works even if Numpy decided to switch again from C to F order or vice versa. All I assume now is that a temporary matrix is C order, but that I can control when constructing the matrix. Thanks for helping on this one! Kind regards, -michael 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
"Klemm, Michael" <michael.klemm@intel.com> wrote:
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.
NumPy's dot function will call cblas_dgemm in the most efficient way regardless of storage. It knows what to do with C and Fortran arrays. Sturla
participants (4)
-
Ian Henriksen -
Klemm, Michael -
Sebastian Berg -
Sturla Molden