[SciPy-User] How to efficiently do dot(dot( A.T, diag(d) ), A ) ?
Hugh Perkins
hughperkins@gmail....
Tue Sep 11 23:56:55 CDT 2012
(pressed send by accident too soon, sorry :-(
The good news is: if I use larger matrices, which is actually closer
to my target scenario, then it seems the scipy with openblas is *very*
good on a
dual-core machine:
matlab:
>> n = 200000;
k = 100;
a = spdiags(rand(n,1),0,n,n);
c = rand(k,n);
tic, d = c*a*c'; toc
tic, d = c*a*c'; toc
Elapsed time is 1.322931 seconds.
Elapsed time is 1.308237 seconds.
python:
python testmult.py
Elapsed time: 0.408543109894
Elapsed time: 1.15053105354
However, on a 12-core machine, it seems that scipy/openblas is rather
less competitive:
matlab:
>> n = 200000;
>> k = 100;
>> a = spdiags(rand(n,1),0,n,n);
>> c = rand(k,n);
>> tic, d = c*a; toc
Elapsed time is 0.106401 seconds.
>> tic, d = d*c'; toc
Elapsed time is 0.098882 seconds.
python:
$ python testmult2.py
Elapsed time: 0.168606996536
Elapsed time: 0.584333896637
python code:
from __future__ import division
from scipy import *
import scipy.sparse as sparse
import time
start = 0
def tic():
global start
start = time.time()
def toc():
global start
print "Elapsed time: " + str( time.time() - start )
start = time.time()
n = 200000;
k = 100;
a = sparse.spdiags( rand(n), 0, n, n )
c = rand(n,k);
tic(); d = c.T * a; toc()
tic(); e = dot( d, c ); toc()
Note that openblas seems to use multithreading by default:
$ export OMP_NUM_THREADS=1
$ python testmult2.py
Elapsed time: 0.145273923874
Elapsed time: 0.750689983368
$ unset OMP_NUM_THREADS
$ python testmult2.py
Elapsed time: 0.168558120728
Elapsed time: 0.578655004501
$ export OMP_NUM_THREADS=12
$ python testmult2.py
Elapsed time: 0.168849945068
Elapsed time: 0.591191053391
... but this is apparently not enough to get quite as good as matlab :-(
Looking at the ratio of the times for 12 threads and 1 thread, it
looks like the openblas performance doesn't change much with the
number of threads.
... but perhaps this is now out of the scope of scipy, and is
something I should check with openblas?
More information about the SciPy-User
mailing list