# How to speed up band matrix-matrix multiplication?

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


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.

# 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 = DeveloperPartitionMap[If[Length@#@"NonzeroValues" > 0, m1@Normal@#, 0] &,
a, {blockSize, blockSize}];
a2.u_ ^:= Flatten[DeveloperToPackedArray[
part.DeveloperPartitionMap[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

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.

• @Mr.Wizard passed ahead of me with a question of heads in SparseArray while I fought with MKL =) – ybeltukov Sep 6 '14 at 18:01
• 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`) – Simon Woods Sep 6 '14 at 18:37
• @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? – ybeltukov Sep 6 '14 at 18:54