Skip to content

Instantly share code, notes, and snippets.

@haampie
Created July 30, 2020 13:43
Show Gist options
  • Save haampie/3da8dd99f1433affc15739b4e597ece5 to your computer and use it in GitHub Desktop.
Save haampie/3da8dd99f1433affc15739b4e597ece5 to your computer and use it in GitHub Desktop.
expression templates
using LinearAlgebra, LoopVectorization
function error_norm(u::AbstractArray{T}, v::AbstractArray{T}) where T
result = T(0)
@avx for i = eachindex(u)
result += (u[i] - v[i])^2
end
return sqrt(result)
end
function finite_diff!(u::AbstractMatrix{T}, v::AbstractMatrix{T}) where T
copyto!(v, u)
@avx @views u[2:end-1,2:end-1] .= (
(v[1:end-2,2:end-1] .+ v[3:end,2:end-1] .+ v[2:end-1,1:end-2] .+ v[2:end-1,3:end]) .* T(4)
.+ v[1:end-2,1:end-2] .+ v[1:end-2,3:end] .+ v[3:end,1:end-2] .+ v[3:end,3:end]) ./ T(20)
return error_norm(u, v)
end
function run(num, ::Type{T} = Float64) where T
u = zeros(T, num, num)
v = similar(u)
x = sin.((0:num-1) .* T(π) ./ T(num-1))
u[:, 1] .= x
u[:, end] .= x .* exp(-T(π))
iter, err = 0, typemax(T)
while iter < 100000 && err > 1e-6
err = finite_diff!(u, v)
iter += 1
end
println("Relative error is: ", err)
println("Number of iterations: ", iter)
end
@haampie
Copy link
Author

haampie commented Jul 30, 2020

julia> @time run(100)
Relative error is: 9.999189386211492e-7
Number of iterations: 15231
  0.096602 seconds (26 allocations: 158.125 KiB)

julia> @time run(150)
Relative error is: 9.997364776246965e-7
Number of iterations: 32971
  0.461015 seconds (26 allocations: 353.953 KiB)

julia> @time run(200)
Relative error is: 9.999117654756488e-7
Number of iterations: 56877
  1.564764 seconds (26 allocations: 627.766 KiB)
$ ./eigen 100
 Relative error is:  9.99919e-07 
 Number of iterations:  15231 
 Elapsed time is: 0.104035 seconds 
$ ./eigen 150
 Relative error is:  9.99736e-07 
 Number of iterations:  32971 
 Elapsed time is: 0.519756 seconds 
$ ./eigen 200
 Relative error is:  9.99912e-07 
 Number of iterations:  56877 
 Elapsed time is: 1.56084 seconds 

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment