# How to compile code that involve sparse array?

I have the following code

Clear[test];
test[i_, hh_, NN_] := Module[{v0, v1, vtmp},
dim = Length[hh];
u = ConstantArray[0., NN];
v0 = SparseArray[i->1.,{dim}];
v1 = hh.v0;
u[[1]] = 1.;
u[[2]] = v1[[i]];
Do[vtmp = 2 hh.v1 - v0; u[[k]] = vtmp[[i]]; v0 = v1;
v1 = vtmp;, {k, 3, NN}]; u]


Now hh is a sparse matrix

hh = SparseArray[{Band[{1, 2}] -> 0.1, Band[{2, 1}] -> 0.1}, {1000, 1000}]/0.201;


and this is timing by varying NN

Table[test[1, hh, 10^i]; // AbsoluteTiming, {i, 1, 5}]
(*{{0.000357558, Null}, {0.00196176, Null}, {0.0288201,
Null}, {0.344731, Null}, {3.32809, Null}}*)


Kind of a slow when NN=100000. So I need to speed it up.

But Mathematica's Compile seems not supporting sparse array, since simple

ff = Compile[{{x, _Real, 2}}, x*x];
ff[hh]


"Argument SparseArray[Automatic,{1000,1000},0.,{<<1>>}] at position 1 should be a rank 2 tensor of "machine-size real number"

So is it possible to compile or any other way to speed up this code?

update

I just found if I replace SparseArray[i->1.,{dim}] with UnitVector[dim,i], then the timing becomes

In[8]:= test[1, hh, 10^5]; // AbsoluteTiming

Out[8]= {1.96876, Null}


much faster but quite odd, why faster?

update2

even weirder bug behaviour is this.

define

comtest =
Compile[{{n, _Integer}, {NN, _Integer}}, Sin[(\[Pi] n)/(NN + 1)],
RuntimeOptions -> "Speed", RuntimeAttributes -> {Listable}]


and then run this once

comtest[Range[0, 10^5 - 1], 10^5] test[1, hh, 10^5];


Now the timing becomes

In[10]:= test[1, hh, 10^5]; // AbsoluteTiming

Out[10]= {5.80249, Null}


This is ... ridiculous. I tried many times, always this result. I am using mathematica 10.4.

How could two irrelavent evaluation affect each other??!!!

• Only packed arrays are admissible in Compile[]. But, why not discuss what you're trying to implement here, and we can see if there are better ways to get what you want? – J. M. is away Apr 6 '16 at 11:40
• @J.M. Hi, J.M. The most timing consuming part as you can see is the do loop. v0 and v1 are vectors, we can get new vector from v0 and v1 by 2 hh.v1 - v0 and extract element at position i of it, stored it in List u. and then the v0 is replaced by v1, and v1 by the new vector, and it goes on and on, until we get a whole list of u, so approximately NN iterations – matheorem Apr 6 '16 at 11:47
• Okay, but I was asking about the algorithm itself; what are you trying to do with your matrix? – J. M. is away Apr 6 '16 at 11:48
• @J.M. I am implementing this algorithm arxiv.org/abs/cond-mat/0504627 The most time consuming part is iteration of Multiplication of sparse matrix and vector. Because I know do loop is slow in MMA, so I want to compile it and expecting it should be a few times faster. – matheorem Apr 6 '16 at 12:01
• Hmm... try this uncompiled version first: test[i_, hh_, NN_] := Module[{dim = Length[hh], ub, v0, v1}, v0 = UnitVector[dim, i]; v1 = hh.v0; ub = InternalBag[{1., hh[[i, i]]}]; Do[{v0, v1} = {v1, 2 hh.v1 - v0}; InternalStuffBag[ub, v1[[i]]], {NN - 2}]; InternalBagPart[ub, All]]` – J. M. is away Apr 6 '16 at 12:04