I have a band matrix

$HistoryLength = 0;
n = 10000;
b = 300;
k = 300;
a = SparseArray[Flatten[#, 1] &@Table[{i, Mod[i + j, n, 1]}, {i, n}, {j, -b, b}] ->
     RandomComplex[{-1 - I, 1 + I}, n (2 b + 1)]];

and a dense matrix

u = RandomComplex[{-1 - I, 1 + I}, {n, k}];

I want to multiply them as fast as possible

v = a.u; // AbsoluteTiming

{6.748518, Null}

Visual representation of this multiplication:

draw = ArrayPlot[#[[;; ;; 30, ;; ;; 30]], ImageSize -> {Automatic, 200}] &;
Row@{draw[v], " = ", draw[a].draw[u]}

enter image description here

This problem usually comes up when you want to multiply a Hamiltonian by a set of wavefunctions (in a certain basis).

Why I expect a possibility of speeding up? When you multiply dense matrices you can use algorithms like Strassen algorithm and use the processor cache to operate with small blocks. The matrix a have a dense band. This knowledge can increase performance in contradiction to the sparse matrix of a general form.


1 Answer 1


Dense blocks

One can split matrices by blocks and use these blocks as a dense matrices

blockSize = 100;
m1[a_] m2[b_] ^:= a.b; 
part = Developer`PartitionMap[If[Length@#@"NonzeroValues" > 0, m1@Normal@#, 0] &, 
    a, {blockSize, blockSize}];
a2.u_ ^:= Flatten[Developer`ToPackedArray[
    part.Developer`PartitionMap[m2, u, blockSize]], 1]

I have to introduce the intermediate headers m1 and m2 since SparseArray doesn't allow List header. Block size should be several times smaller then the width of the band.

v2 = a2.u; // AbsoluteTiming

{2.581961, Null}

Almost 3x speed up from nothing!


MKL has a special function mkl_zdiamm for the band matrix - dense matrix multiplication.

However it is 5 times slower then build-in multiplication... But now I know how to use MKL from Mathematica. There was some problems and I'm going to discuss them in a separate question.

  • $\begingroup$ @Mr.Wizard passed ahead of me with a question of heads in SparseArray while I fought with MKL =) $\endgroup$
    – ybeltukov
    Sep 6, 2014 at 18:01
  • $\begingroup$ What version are you using? With version 10 on Windows I get 1.75 secs for a.u and 1.2 secs for a2.u (but with an added overhead of 2.4 secs to calculate part) $\endgroup$ Sep 6, 2014 at 18:37
  • $\begingroup$ @SimonWoods I use version 10 on Linux. I tested on my laptop with Core 2 Quad. Desktop with Core i7 shows even more speed up. Overhead doesn't matter in real calculations, when you multiply a fixed Hamiltonian by wavevectors over and over again. Did you test with the same parameters? $\endgroup$
    – ybeltukov
    Sep 6, 2014 at 18:54

This site is temporarily in read-only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .