Write a smoother for AlgebraicMultigrid

New to Julia
,
1

I’m playing around with https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl and would like to write my own smoother, or even just use the Jacobi smoother which is in the package, but not the default. Unfortunately, I am still very new to Julia and I still don’t understand how or why they’ve set up the smoothers the way they have. Here’s how they have defined the Jacobi smoother

abstract type Smoother end
struct Jacobi{T,TX} <: Smoother
    ω::T
    temp::TX
    iter::Int
end
Jacobi(ω, x::TX; iter=1) where {T, TX<:AbstractArray{T}} = Jacobi{T,TX}(ω, similar(x), iter)
Jacobi(x::TX, ω=0.5; iter=1) where {T, TX<:AbstractArray{T}} = Jacobi{T,TX}(ω, similar(x), iter)

function (jacobi::Jacobi)(A, x, b)
...
end

What is the advantage of defining a Smoother type instead of allowing any function?

When I tried to use Jacobi as the smoother

ruge_stuben(A,
        presmoother = Jacobi(1, Array{Float64,1}, 1),
        postsmoother = Jacobi(1, Array{Float64,1}, 1))

I get an error, MethodError: no method matching setindex!(::Type{Array{Float64,1}}, ::Float64, ::Int64, ::Int64).

Some assistance and background explaining why this is a good way to set up functions in Julia would be helpful!

2

The second argument is expected to be an Array but you’ve given a type instead. x should be a Vector{Float64}(length).

3

Thanks! So is this the normal way to set up a temporary array for internal use in a function in Julia? Is there a way to size an array like this dynamically the first time a function is called so you don’t need to fix the length ahead of time?

4

First, I made a mistake that uninitialized array is allocated like Array{Float64}(undef, ...sizes in each dimension...).

It depends on who defines the callable. If you define a custom callable you can check the length of a vector and and use resize! to the needed length but it doesn’t work for higher dimensional Arrays. However, the aim for ahead of time allocation is to not check later for the size.

I think the x in Jacobi(ω, x::TX; iter=1) or Jacobi(x::TX, ω=0.5; iter=1) is expected to be the same x in the callable (jacobi::Jacobi)(A, x, b) so it allocates the compatible temporary array with similar(x) on the left side of the outer constructors. In this case, if A is square then size(x) == size(b) so Jacobi(1, b, iter=1) can be used if b is already known but x is not defined by you but the solver internally. If A is not square, the smoother is something like Jacobi(1, Array{Float64}(undef, size(A, 2), size(b, 1)), 1).

One more thing, Jacobi(1, vec, 1) differs from Jacobi(1, vec, iter=1). The latter uses the outer constructor with keyword argument iter but the former uses the default constructor with 3-argument that directly initializes fields with the corresponding arguments. Meaning, while in the former vec becomes the temporary array, the latter allocates another array with the shape of vec for the temporary.

5

Got it. Thanks for the help!