Comparing ParSpMatVec to Julia's sparse matvecs

In [56]:

    

"""
L,Lin,Lib = getLaplacian(n::Int)

returns discrete 2D Laplacian, interior part, and boundary part, on regular mesh with n^2 cells.
"""
function getLaplacian(n::Int)
    h  = 1/n
    dx = spdiagm((-ones(n,1),ones(n,1)),0:1,n,n+1)/h
    d2x = dx'*dx
    L  = kron(speye(n+1),d2x) + kron(d2x,speye(n+1))
   return L
end


    

    




WARNING: Method definition getLaplacian(Int64) in module Main at In[54]:7 overwritten at In[56]:7.
WARNING: replacing docs for 'getLaplacian :: Tuple{Int64}' in module 'Main'.






    
Out[56]:





getLaplacian

In [57]:

    

using ParSpMatVec
using Base.Test

# modify these parameters to test different cases
testRandom    = false    # sprandn or 2D Laplacian
testTranspose = true     # test A'*x or A*x
nvec          = 1        # number of vectors to multiply 
numProcs      = 1        # number of processors (for ParSpMatVec)
nrun          = 20       # number of runs

if testTranspose
    f1 = Base.LinAlg.At_mul_B!
    f2 = ParSpMatVec.Ac_mul_B!
else
    f1 = Base.LinAlg.A_mul_B!
    f2 = ParSpMatVec.A_mul_B!
end

if testRandom
    getA = n-> sprandn(n,n,5/n)
else
    getA = n-> getLaplacian(round(Int64,sqrt(n)))
end

N     = logspace(4,7,21)
time1 = zeros(length(N),nrun)
time2 = zeros(length(N),nrun)

for k=1:length(N)
    n = round(Int64,N[k])
    
    A = getA(n)

    x = rand(size(A,2),nvec);  x = x*10 - 5;
    y = rand(size(A,2),nvec);  y = y*10 - 5;
    
    if nvec==1
        x = vec(x)
        y = vec(y)
    end

    alpha = 123.56
    beta = 543.21
    y2 = copy(y)
    
    # use julia's matvec
    for j=1:nrun
        tic(); 
            y2 = f1(alpha,A,x,beta,y2)
        time1[k,j] = toq()
    end

    # use Fortran code
    y3 = copy(y)
    tic(); 
    for j=1:nrun
        tic()
        f2( alpha, A, x, beta, y3, numProcs ); 
        time2[k,j] = toq()
    end
    
    # print results
    @printf "n=%d\t Base=%1.4f\t ParSpMatVec=%1.4f\t speedup=%1.4f\n" n mean(time1[k,:]) mean(time2[k,:]) mean(time1[k,:])/mean(time2[k,:])
    @test norm(y3-y2) / norm(y) < 1.e-12
end


    

    




n=10000	 Base=0.0001	 ParSpMatVec=0.0001	 speedup=1.3283
n=14125	 Base=0.0001	 ParSpMatVec=0.0001	 speedup=1.6935
n=19953	 Base=0.0002	 ParSpMatVec=0.0001	 speedup=1.2694
n=28184	 Base=0.0002	 ParSpMatVec=0.0003	 speedup=0.6741
n=39811	 Base=0.0003	 ParSpMatVec=0.0002	 speedup=1.5173
n=56234	 Base=0.0006	 ParSpMatVec=0.0003	 speedup=1.7774
n=79433	 Base=0.0006	 ParSpMatVec=0.0005	 speedup=1.3409
n=112202	 Base=0.0011	 ParSpMatVec=0.0007	 speedup=1.6192
n=158489	 Base=0.0013	 ParSpMatVec=0.0010	 speedup=1.2669
n=223872	 Base=0.0018	 ParSpMatVec=0.0017	 speedup=1.0566
n=316228	 Base=0.0026	 ParSpMatVec=0.0020	 speedup=1.2822
n=446684	 Base=0.0039	 ParSpMatVec=0.0032	 speedup=1.2020
n=630957	 Base=0.0054	 ParSpMatVec=0.0049	 speedup=1.0955
n=891251	 Base=0.0075	 ParSpMatVec=0.0065	 speedup=1.1540
n=1258925	 Base=0.0126	 ParSpMatVec=0.0092	 speedup=1.3721
n=1778279	 Base=0.0156	 ParSpMatVec=0.0149	 speedup=1.0464
n=2511886	 Base=0.0235	 ParSpMatVec=0.0228	 speedup=1.0286
n=3548134	 Base=0.0357	 ParSpMatVec=0.0300	 speedup=1.1888
n=5011872	 Base=0.0604	 ParSpMatVec=0.0438	 speedup=1.3790
n=7079458	 Base=0.0775	 ParSpMatVec=0.0628	 speedup=1.2336
n=10000000	 Base=0.1044	 ParSpMatVec=0.0924	 speedup=1.1298

In [ ]: