Sparse matrix assembly with bad timings

General Usage Performance
1

Hi,

I took the 2016 codes from Julia iFEM 2: Optimizing Stiffness Matrix Assembly and made some cosmetic changes to run them with Julia 1.7.

Chris’ old timings were:

tic; assemblingstandard(node,elem); toc;
tic; assemblingsparse(node,elem); toc;
tic; assembling(node,elem); toc;
 
Elapsed time is 0.874919 seconds.
Elapsed time is 0.698763 seconds.
Elapsed time is 0.020453 seconds.

mine:

  7.117206 seconds (1.17 M allocations: 880.707 MiB, 8.60% gc time, 18.66% compilation time)
  0.872631 seconds (1.09 M allocations: 90.582 MiB, 23.57% compilation 
time)
  1.388404 seconds (271.57 k allocations: 41.201 MiB, 93.56% compilation time)

Where did the speed of assembling() go? Where is the error?
Thx

UPDATE, using BenchmarkTools:

 4.168 s (1020002 allocations: 872.65 MiB)
 982.999 ms (1020033 allocations: 86.91 MiB)
 86.067 ms (259 allocations: 29.39 MiB)

using LinearAlgebra
using SparseArrays

# http://www.stochasticlifestyle.com/julia-ifem-1-mesh-generation-and-julias-type-system/
function squaremesh(square,h)
    #square = [x0 x1 y0 y1], h is stepsize
    x0 = square[1]; x1= square[2]
    y0 = square[3]; y1= square[4]
    x,y = meshgrid(x0:h:x1,y0:h:y1)
    node = [x[:] y[:]];
    ni = size(x,1); # number of rows
    N = size(node,1);
    t2nidxMap = 1:N-ni;
    topNode = ni:ni:N-ni;
    t2nidxMap = deleteat!(collect(t2nidxMap),collect(topNode))
    k = t2nidxMap;
    elem = [k.+ni k.+(ni+1) k ; k.+1 k k.+(ni+1)]
    return node,elem
  end

# http://www.stochasticlifestyle.com/julia-ifem2/
function localstiffness(p)
    At = Array{Float64}(undef, 3,3)
    B = [p[1,:]-p[3,:]; p[2,:]-p[3,:]];
    B = reshape(B,2,2)       # line added
    G = [1 0 ; 0 1 ; -1 -1]';
    area = 0.5*abs(det(B));
    for i = 1:3
      for j = 1:3
        At[i,j] = area*dot(B\G[:,i], B\G[:,j]); 
      end
    end
    return(At,area)
  end

function assemblingstandard(node,elem)
    N=size(node,1); NT=size(elem,1);
    A=zeros(N,N);
    for t=1:NT
      At,=localstiffness(node[vec(elem[t,:]),:]);
      for i=1:3
        @simd for j=1:3
          @inbounds A[elem[t,i],elem[t,j]]=A[elem[t,i],elem[t,j]]+At[i,j];
        end
      end
    end
    return(A)
  end

  function assemblingsparse(node,elem)
    N = size(node,1); 
    NT = size(elem,1);
    i = Array{Int64}(undef, NT,3,3);
    j = Array{Int64}(undef, NT,3,3);
    s = zeros(NT,3,3);
    index = 0;
    for t = 1:NT
      At, = localstiffness(node[vec(elem[t,:]),:]);
      for ti = 1:3, tj = 1:3
          i[t,ti,tj] = elem[t,ti];
          j[t,ti,tj] = elem[t,tj];
          s[t,ti,tj] = At[ti,tj];
      end
    end
    return(sparse(vec(i), vec(j), vec(s)));
  end

  function assembling(node,elem)
    N = size(node,1); NT = size(elem,1);
    ii = Vector{Int64}(undef, 9*NT);
    jj = Vector{Int64}(undef, 9*NT);
    sA = Vector{Float64}(undef, 9*NT);
    ve = Array{Float64}(undef, NT,2,3)
    ve[:,:,3] = node[vec(elem[:,2]),:]-node[vec(elem[:,1]),:];
    ve[:,:,1] = node[vec(elem[:,3]),:]-node[vec(elem[:,2]),:];
    ve[:,:,2] = node[vec(elem[:,1]),:]-node[vec(elem[:,3]),:];
    area = 0.5*abs.(-ve[:,1,3].*ve[:,2,2]+ve[:,2,3].*ve[:,1,2]);
    index = 0;
    for i = 1:3, j = 1:3
        @inbounds begin
            ii[index+1:index+NT] = elem[:,i]
            jj[index+1:index+NT] = elem[:,j]
            sA[index+1:index+NT] = sum(ve[:,:,i].* ve[:,:,j], dims=2) ./(4*area) # Replacing dot(ve[:,:,i],ve[:,:,j],2)
            index = index + NT;
        end
    end
    return(sparse(ii,jj,sA))
  end

node,elem = squaremesh([0 1 0 1], 0.01)
@time assemblingstandard(node,elem);
@time assemblingsparse(node,elem);
@time assembling(node,elem);
2

? You might try to run the timing twice, first for compilation and then for run time.

3

@goerch I thought I did re-run, but with BenchmarkTools you are on the safer side.
Thanks a lot!

1 Like