I have the following code

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];

will lead to error message

"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?


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?


even weirder bug behaviour is this.


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??!!!

  • $\begingroup$ 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? $\endgroup$ – J. M. is away Apr 6 '16 at 11:40
  • $\begingroup$ @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 $\endgroup$ – matheorem Apr 6 '16 at 11:47
  • $\begingroup$ Okay, but I was asking about the algorithm itself; what are you trying to do with your matrix? $\endgroup$ – J. M. is away Apr 6 '16 at 11:48
  • $\begingroup$ @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. $\endgroup$ – matheorem Apr 6 '16 at 12:01
  • $\begingroup$ 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 = Internal`Bag[{1., hh[[i, i]]}]; Do[{v0, v1} = {v1, 2 hh.v1 - v0}; Internal`StuffBag[ub, v1[[i]]], {NN - 2}]; Internal`BagPart[ub, All]] $\endgroup$ – J. M. is away Apr 6 '16 at 12:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.