expression templates

julia.jl

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

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