Instantly share code, notes, and snippets.

# haampie/julia.jl

Created July 30, 2020 13:43
Show Gist options
• Save haampie/3da8dd99f1433affc15739b4e597ece5 to your computer and use it in GitHub Desktop.
expression templates
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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 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
``````