Benchmark for BiCGSTB

As an example, we solve the curl-curl system arising in electromagnetics

$$ \nabla \times \nabla \times u + \imath w u = 0 $$

with natural boundary conditions. To generate the discrete PDE, we use a regular mesh and the mimetic finite volume discretization provided by jInv.Mesh

A good solver for these kinds of problems is BiCGSTAB and we test the implemenation from KrylovMethods for growing mesh sizes. If you have another implementation please add it here and make a PR.



In [ ]:

    
using jInv.Mesh
using KrylovMethods
using ParSpMatVec
using Base.Test

Ns        = 16*[1,2,3,4]  # number of grid points
nTrials   = 10            # number of repititions
nthreads  = 4             # number of threads used to accelerate sparse matvecs

# allocate space for results
timesKM    = zeros(length(Ns),nTrials)
timesK     = zeros(length(Ns),nTrials)
timesKMOpt = zeros(length(Ns),nTrials)
timesIS    = zeros(length(Ns),nTrials);



In [ ]:

    
relresKM = 0
iterKM = 0
iterKMOpt = 0
relresKMOpt = 0
for k=1:length(Ns)
    n = Ns[k]
    
    # build Curl Curl system
    M = getRegularMesh([0 1 0 1 0 1],[n,n,n])
    CURL = getCurlMatrix(M)
    w    = 1
    A    = CURL'*CURL - 1im*w*I(size(CURL,2))
    rhs  = randn(size(A,2)) + 1im *randn(size(A,2))
    
    if k==1 # warmup
        resKM = KrylovMethods.bicgstb(A,rhs,tol=1e-4,maxIter=2,out=-1);
    end
    for j=1:nTrials
        tic()
        resKM = KrylovMethods.bicgstb(A,rhs,tol=1e-12,maxIter=100,out=-1)
        timesKM[k,j] = toq();
        iterKM = resKM[4]
        relresKM = resKM[3]

        tic()
        yt = zeros(Complex128,size(A,1))
        Afun(x) = (yt[:]=Complex128(0.0); ParSpMatVec.Ac_mul_B!( Complex128(1.0), A, x, Complex128(0.0), yt, nthreads); return yt)
        resKMOpt = KrylovMethods.bicgstb(Afun,rhs,tol=1e-12,maxIter=100,out=-1)
        timesKMOpt[k,j] = toq();
        iterKMOpt   = resKM[4]
        relresKMOpt = resKMOpt[3]       
        @test iterKMOpt==iterKM
#         @test abs(relresKMOpt-relresKM)/relresKM < 1e-1
    end
println("n=$n, KM.jl: $(mean(timesKM[k,:]))  KM.jl+ParSpMat:$(mean(timesKMOpt[k,:]))")
end



In [ ]:

    
using PyPlot



In [ ]:

    
loglog(Ns,mean(timesKM,2),"--b",linewidth=3)
hold(true)
loglog(Ns,mean(timesKMOpt,2),"-b",linewidth=3)
xlabel("degrees of freedom")
ylabel("runtime")
title("Comparison of PCG runtimes for 3D CurlCurl")
legend(("KrylovMethods.jl","KrylovMethods+ParSpMatVec"),loc=4)



In [ ]: