On PackedArray, looking for advice for using them

Question

I have not used PackedArray before, but just started looking at using them from reading some discussion on them here today.

What I have is lots of large size 1D and 2D matrices of all reals, and no symbolic (it is a finite difference PDE solver), and so I thought that I should take advantage of using PackedArray.

I have an initialization function where I allocate all the data/grids needed. So I went and used ToPackedArray on them. It seems a bit faster, but I need to do more performance testing to better compare speed before and after and also compare RAM usage.

But while I was looking at this, I noticed that some operations in M automatically return lists in PackedArray already, and some do not.

For example, this does not return packed array

a = Table[RandomReal[], {5}, {5}];
Developer`PackedArrayQ[a]

But this does

a = RandomReal[1, {5, 5}];
Developer`PackedArrayQ[a]

and this does

a = Table[0, {5}, {5}];
b = ListConvolve[ {{0, 1, 0}, {1, 4, 1}, {0, 1, 1}}, a, 1];
Developer`PackedArrayQ[b]

and also matrix multiplication does return result in packed array

a = Table[0, {5}, {5}];
b = a.a;
Developer`PackedArrayQ[b]

But element wise multiplication does not

b = a*a;
Developer`PackedArrayQ[b]

My question : Is there a list somewhere which documents which M commands return PackedArray vs. not? (assuming data meets the requirements, such as Real, not mixed, no symbolic, etc..)

Also, a minor question, do you think it will be better to check first if a list/matrix created is already packed before calling calling ToPackedArray on it? I would think calling ToPackedArray on list already packed will not cost anything, as the call will return right away.

thanks,

update (1)

Just wanted to mention, that just found that PackedArray symbols not allowed in a demo CDF as I got an error uploading one with one. So, had to remove all my packing code out. Since I mainly write demos, now this topic is just of an academic interest for me. But wanted to thank everyone for time and good answers.

Answer 1

There isn't a comprehensive list. To point out a few things:

• Basic operations with packed arrays will tend to remain packed:

 
    In[66]:= a = RandomReal[1, {5, 5}];

    In[67]:= Developer`PackedArrayQ /@ {a, a.a, a*a}

    Out[67]= {True, True, True}

• Note above that that my version (8.0.4) doesn't unpack for element-wise multiplication.

• Whether a Table will result in a packed array depends on the number of elements:

 
    In[71]:= Developer`PackedArrayQ[Table[RandomReal[], {24}, {10}]]

    Out[71]= False

    In[72]:= Developer`PackedArrayQ[Table[RandomReal[], {24}, {11}]]

    Out[72]= True

    In[73]:= Developer`PackedArrayQ[Table[RandomReal[], {25}, {10}]]

    Out[73]= True

• On["Packing"] will turn on messages to let you know when things unpack:

 
    In[77]:= On["Packing"]

    In[78]:= a = RandomReal[1, 10];

    In[79]:= Developer`PackedArrayQ[a]

    Out[79]= True

    In[80]:= a[[1]] = 0 (* force unpacking due to type mismatch *)

       Developer`FromPackedArray::punpack1: Unpacking array with dimensions {10}. >>

    Out[80]= 0

• Operations that do per-element inspection will usually unpack the array,

    In[81]:= a = RandomReal[1, 10];

    In[82]:= Position[a, Max[a]]

       Developer`FromPackedArray::unpack: Unpacking array in call to Position. >>

    Out[82]= {{4}}

• There penalty for calling ToPackedArray on an already packed list is small enough that I wouldn't worry about it too much:

    In[90]:= a = RandomReal[1, 10^7];

    In[91]:= Timing[Do[Identity[a], {10^5}];]

    Out[91]= {0.028089, Null}

    In[92]:= Timing[Do[Developer`ToPackedArray[a], {10^5}];]

    Out[92]= {0.043788, Null}

• The frontend prefers packed to unpacked arrays, which can show up when dealing with Dynamic and Manipulate:

    In[97]:= Developer`PackedArrayQ[{1}]

    Out[97]= False

    In[98]:= Dynamic[Developer`PackedArrayQ[{1}]]

    Out[98]= True

• When looking into performance, focus on cases where large lists are getting unpacked, rather than the small ones. Unless the small ones are in big loops.

Answer 2

This is just an addendum to Brett's answer:

SystemOptions["CompileOptions"]

will give you the lengths being used for which a function will return a packed array. So if you did need to pack a small list, as an alternative to using Developer`ToPackedArray you could temporarily set a smaller number for one of the compile options. e.g.

SetSystemOptions["CompileOptions" -> {"TableCompileLength" -> 20}]

Note also some difference between functions which to me at least doesn't seem intuitive so I generally have to test these kind of things whenever I use them rather than instinctively knowing what will work best:

f = # + 1 &;
g[x_] := x + 1;
data = RandomReal[1, 10^6];

On["Packing"]
Timing[Developer`PackedArrayQ[f /@ data]]
{0.131565, True}


Timing[Developer`PackedArrayQ[g /@ data]]
Developer`FromPackedArray::punpack1: Unpacking array with dimensions {1000000}.
{1.95083, False}