I really like the new iteration protocol, and I’m happy to see it being spread across packages. However, I usually have the problem, that I work with multidimensional arrays that is often not supported directly. In the MWE below, I use
gseidel! from IterativeSolvers.jl. The method accepts 1-and 2-d arrays but mine are 2-and 3d arrays, thus I need to copy and slice arrays to fit within a for loop.
Are there more idiomatic ways of doing the same?
using IterativeSolvers n = 10 m = 3 A = rand(n,n,m) b = rand(n,m) x1 = zeros(n) # 2d @time gauss_seidel!(x1, A[:,:,1], b[:,1]; maxiter=10) # 3d x2 = zeros(n,m) xt = zeros(n) @time for i = 1:m xt .= x2[:,i] gauss_seidel!(xt, A[:,:,i], b[:,i]; maxiter=10) x2[:,i] .= xt end x1 == x2[:,1]
Benchmarking gives 3 and 27 allocations, respectively.