5
$\begingroup$

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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
6
$\begingroup$

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

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ @Mr.Wizard passed ahead of me with a question of heads in SparseArray while I fought with MKL =) $\endgroup$
    – ybeltukov
    Commented 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$
    – Simon Woods
    Commented 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
    Commented Sep 6, 2014 at 18:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.