<br>Hi Christopher,<br><br>thanks for taking the time to reply at length. I do understand the concept of striding in general but was not familiar with the  Numpy way of accessing that information. So thanks for pointing me to .flag and .stride.<br>
<br>That said, BLAS/LAPACK do have apis that take the stride length into account. But for sparse arrays I think its a hopeless situation. That is a bummer, because sparse is what I need. Oh well, I will probably do it in C++<br>
<br>-- srean<br><br>p.s. I hope top posting is not frowned upon here. If so, I will keep that in mind in my future posts.<br><br><div class="gmail_quote">On Sat, Mar 26, 2011 at 1:31 PM, Christopher Barker <span dir="ltr"><<a href="mailto:Chris.Barker@noaa.gov">Chris.Barker@noaa.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><div class="im">
<br>
</div>Probably not -- the trick is that when an array is a view of a slice of<br>
another array, it may not be laid out in memory in a way that other libs<br>
(like LAPACK, BLAS, etc) require, so the data needs to be copied to call<br>
those routines.<br>
<br>
To understand all this, you'll need to study up a bit on how numpy<br>
arrays lay out and access the memory that they use: they use a concept<br>
of "strided" memory. It's very powerful and flexible, but most other<br>
numeric libs can't use those same data structures. I"m not sure what a<br>
good doc is to read to learn about this -- I learned it from messing<br>
with the C API. TAke a look at any docs that talk about "strides", and<br>
maybe playing with the "stride tricks" tools will help.<br>
<br>
A simple example:<br>
<br>
In [3]: a = np.ones((3,4))<br>
<br>
In [4]: a<br>
Out[4]:<br>
array([[ 1.,  1.,  1.,  1.],<br>
        [ 1.,  1.,  1.,  1.],<br>
        [ 1.,  1.,  1.,  1.]])<br>
<br>
In [5]: a.flags<br>
Out[5]:<br>
   C_CONTIGUOUS : True<br>
   F_CONTIGUOUS : False<br>
   OWNDATA : True<br>
   WRITEABLE : True<br>
   ALIGNED : True<br>
   UPDATEIFCOPY : False<br>
<br>
So a is a (3,4) array, stored in C_contiguous fashion, jsut like a<br>
"regular old C array". A lib expecting data in this fashion could use<br>
the data pointer just like regular C code.<br>
<br>
In [6]: a.strides<br>
Out[6]: (32, 8)<br>
<br>
this means is is 32 bytes from the start of one row to the next, and 8<br>
bytes from the start of one element to the next -- which makes sense for<br>
a 64bit double.<br>
<br>
<br>
In [7]: b = a[:,1]<br>
<br>
In [10]: b<br>
Out[10]: array([ 1.,  1.,  1.])<br>
<br>
so b is a 1-d array with three elements.<br>
<br>
In [8]: b.flags<br>
Out[8]:<br>
   C_CONTIGUOUS : False<br>
   F_CONTIGUOUS : False<br>
   OWNDATA : False<br>
   WRITEABLE : True<br>
   ALIGNED : True<br>
   UPDATEIFCOPY : False<br>
<br>
but it is NOT C_Contiguous - the data is laid out differently that a<br>
standard C array.<br>
<br>
In [9]: b.strides<br>
Out[9]: (32,)<br>
<br>
so this means that it is 32 bytes from one element to the next -- for a<br>
8 byte data type. This is because the elements are each one element in a<br>
row of the a array -- they are not all next to each other. A regular C<br>
library generally won't be able to work with data laid out like this.<br>
<br></blockquote></div>