Hi all,
New to Julia, but with experience in HPC. I am looking for a canonical way to have module variables that can be shared between the different functions of the module without having to pass these variables as arguments to every function.
The use case is the following. I am solving (many times) sparse linear systems of equations Mx = b. For this I have different solvers: Conjugate Gradient (CG), BiCGstab,
, some problem dependent preconditioning, etc… Each solver needs some some additional allocated memory to actually solve the system of equations (for example, CG needs three adittional vectors).
Usually in fortran or c, one would define these vectors as module private variables, allocate them once, and re-use them in the routine that actually performs the CG. This avoids the allocation/deallocation of memory, that in large sparse systems would kill any performance.
But surprisingly there is no simple way to do something like this in Julia. One can define global variables in a file, but this is produces terrible performance. I always receive the advice to pass this memory as a struct to the function, but this is just terrible programming: There is a long chain of function calls, and one would need to pass the solver workspace to each of them in order to avoid the multiple allocation.
I managed to do something like this using a let construct:
let Ap = similar(Vector{Float64}, 0), r = similar(Vector{Float64}, 0), p = similar(Vector{Float64}, 0)
global function alloc_ws(n::Int64)
Ap = zeros(n)
r = zeros(n)
p = zeros(n)
end
global function CG!(x::Vector{Float64}, func::Function, b::Vector{Float64}, msq::Float64)::Int64
func(Ap, x, msq)
r .= b .- Ap
res_old::Float64 = dot(r,r)
p .= r
niter::Int64 = 0
res_new::Float64 = 0.0
alpha::Float64 = 0.0
for i in 1:length(b)
niter += 1
func(Ap,p,msq)
alpha = res_old / dot(p,Ap)
x .= x .+ alpha .* p
r .= r .- alpha .* Ap
res_new = dot(r,r)
if (res_new < 1.0E-10)
break
end
p .= r .+ (res_new/res_old) .* p
res_old = res_new
end
return niter
end
end
This looks a reasonable solution, since the temporary arrays Ap, r, p are alllocated once. Also the code is type stable (so I do not see any reason for this not to be optimized).
Still I find it strange that such a common situation in practical numerical situations (i.e. a routine that requires some workspace and you want to avoid multiple allocations by using a type variable valid in the scope of the module), has to be done in such a strange way…
Since I am very new to Julia, am I missing something here?