Element stiffness benchmark

keTest.jl

function keTestVec(N::Int64, gradN::Array{Float64,2}, kappa::Array{Float64,2}, keDim::Int64, sdim::Int64)
	ke = zeros(Float64,keDim,keDim)
	detJ = 1.0
	w = 1.0/2.0
	for n = 1:N
		fill!(ke,0.0)
		ke += (w*detJ).*(gradN*kappa*gradN')
	end
	return ke
end

function keTestDeVec(N::Int64, gradN::Array{Float64,2}, kappa::Array{Float64,2}, keDim::Int64, sdim::Int64)
	ke = zeros(Float64,keDim,keDim)
	detJ = 1.0
	w = 1.0/2.0
	preFac = 0.0
	for n = 1:N
		fill!(ke,0.0)
		preFac = w*detJ
		for i = 1:keDim
			for l = 1:keDim
				for j = 1:sdim
					for k = 1:sdim
						ke[i,l] = ke[i,l] + gradN[i,j]*kappa[j,k]*gradN[l,k] 
					end
				end
				ke[i,l] = preFac*ke[i,l]
			end
		end
	end
	return ke
end

function keTestDeVecInbounds(N::Int64, gradN::Array{Float64,2}, kappa::Array{Float64,2}, keDim::Int64, sdim::Int64)
	ke = zeros(Float64,keDim,keDim)
	detJ = 1.0
	w = 1.0/2.0
	preFac = 0.0
	for n = 1:N
		fill!(ke,0.0)
		preFac = w*detJ
		for i = 1:keDim
			for l = 1:keDim
				for j = 1:sdim
					for k = 1:sdim
						@inbounds ke[i,l] = ke[i,l] + gradN[i,j]*kappa[j,k]*gradN[l,k] 
					end
				end
				ke[i,l] = preFac*ke[i,l]
			end
		end
	end
	return ke
end

function keTest(N::Int64)
	keDim = 3
	sdim = 2
	kappa = rand(Float64,sdim,sdim)
	kappa = kappa+kappa'
	gradN = rand(Float64,keDim,sdim)
	println("Vectorized:")
	@time ke1 = keTestVec(N, gradN, kappa, keDim, sdim)
	println("DeVectorized:")
	@time ke2 = keTestDeVec(N, gradN, kappa, keDim, sdim)
	println("DeVectorized InBounds:")
	@time ke3 = keTestDeVecInbounds(N, gradN, kappa, keDim, sdim)
	println("Error norm deVec: $(1-norm(ke2./ke1)/keDim)")
	println("Error norm inBnd: $(1-norm(ke3./ke1)/keDim)")
end