# Write a smoother for AlgebraicMultigrid

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!

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

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?

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 `Array`s. 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.

Got it. Thanks for the help!