# Write a smoother for AlgebraicMultigrid

**URL:** <https://discourse.julialang.org/t/write-a-smoother-for-algebraicmultigrid/41790>\
**Category:** New to Julia\
**Tags:** question, package\
**Created:** [June 20, 2020, 6:54pm UTC](https://discourse.julialang.org/t/write-a-smoother-for-algebraicmultigrid/41790 "2020-06-20T18:54:41Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [June 20, 2020, 6:54pm UTC](https://discourse.julialang.org/t/write-a-smoother-for-algebraicmultigrid/41790/1 "2020-06-20T18:54:41Z")

</div>

I’m playing around with [https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl](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

```nohighlight
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

```nohighlight
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!

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [June 20, 2020, 7:58pm UTC](https://discourse.julialang.org/t/write-a-smoother-for-algebraicmultigrid/41790/2 "2020-06-20T19:58:06Z")

</div>

> [@weymouth](#):
>
> `Jacobi(1, Array{Float64,1}, 1)`

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

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [June 20, 2020, 10:37pm UTC](https://discourse.julialang.org/t/write-a-smoother-for-algebraicmultigrid/41790/3 "2020-06-20T22:37:35Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [June 21, 2020, 8:19am UTC](https://discourse.julialang.org/t/write-a-smoother-for-algebraicmultigrid/41790/4 "2020-06-21T08:19:24Z")

</div>

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

> [@weymouth](#):
>
> 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?

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.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [June 21, 2020, 9:07am UTC](https://discourse.julialang.org/t/write-a-smoother-for-algebraicmultigrid/41790/5 "2020-06-21T09:07:32Z")

</div>

Got it. Thanks for the help!