gher-uliege / DIVAnd.jl

DIVAnd performs an n-dimensional variational analysis of arbitrarily located observations
GNU General Public License v2.0
70 stars 11 forks source link

pcg demanding more memory than direct solver ? #138

Closed jmbeckers closed 9 months ago

jmbeckers commented 11 months ago
using SparseArrays
using LinearAlgebra
using DIVAnd
using Plots
using Statistics
using Random
using BenchmarkTools

Random.seed!(1234)
ND=10^3
NX=45
NY=46
NZ=47
len=0.1
# obs. error variance normalized by the background error variance
epsilon2 = 1.0

fun(x,y,z) = 2*(sin.(6x) * cos.(6y).* cos.(6z)) .+ (0 .- 1) .* x .* y .*z
x = 0.5.+0.25.*randn(ND);
y = 0.5  .+ 0.25 .* randn(ND);
z= 0.5  .+ 0.25 .* randn(ND);
f = fun.(x,y,z)+0.5*randn(ND)

x=[0.5]
y=[0.5]
z=[0.5]
f=[1.0]

# final grid
xi,yi,zi= ndgrid(range(0,stop=1,length=NX),range(0,stop=1,length=NY),range(0,stop=1,length=NZ));

# all points are valid points
mask = trues(size(xi));
#mask[:,:,1].=false
#mask[15:30,15:30,1:end].=false

pm = ones(size(xi)) / (xi[2,1,1]-xi[1,1,1]);
pn = ones(size(xi)) / (yi[1,2,1]-yi[1,1,1]);
po = ones(size(xi)) / (zi[1,1,2]-zi[1,1,1]);

pmn=(pm,pn,po)
xyi=(xi,yi,zi)
xy=(x,y,z)

fi,s= DIVAndrun(mask,pmn,xyi,xy,f,len,epsilon2);

#fi,s=DIVAndrun(mask,pmn,xyi,xy,f,len,epsilon2;inversion=:pcg,maxit=0)

If I replace the direct solver (which fits into memory), by an iterative solver with no iteration (and hence immediate warning of non-convergence), there seems to be, after the warning on non-convergence, a huge increase in memory demand and out of memory in my case:

Julia 1.8.5 and 1.7.1

[1] Array @ .\boot.jl:457 [inlined] [2] Array @ .\boot.jl:466 [inlined] [3] zeros @ .\array.jl:525 [inlined] [4] zeros @ .\array.jl:521 [inlined] [5] getindex(C::CovarIS{Float64, SparseMatrixCSC{Float64, Int64}}, i::Int64, j::Int64) @ DIVAnd C:\Users\jmbeckers.julia\packages\DIVAnd\MV3j9\src\special_matrices.jl:104 [6] isassigned(::CovarIS{Float64, SparseMatrixCSC{Float64, Int64}}, ::Int64, ::Int64) @ Base .\abstractarray.jl:553 [7] _show_nonempty(io::IOContext{IOBuffer}, X::AbstractMatrix, prefix::String, drop_brackets::Bool, axs::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}) @ Base .\arrayshow.jl:438 [8] _show_nonempty(io::IOContext{IOBuffer}, X::CovarIS{Float64, SparseMatrixCSC{Float64, Int64}}, prefix::String) @ Base .\arrayshow.jl:410 [9] show(io::IOContext{IOBuffer}, X::CovarIS{Float64, SparseMatrixCSC{Float64, Int64}}) @ Base .\arrayshow.jl:486 [10] _show_default(io::IOContext{IOBuffer}, x::Any) @ Base .\show.jl:413 [11] show_default @ .\show.jl:396 [inlined] [12] show @ .\show.jl:391 [inlined] [13] show_delim_array(io::IOContext{IOBuffer}, itr::Tuple{Array{Float64, 3}, DIVAnd.DIVAnd_struct{Float64, Int64, 3, SparseMatrixCSC{Float64, Int64}}}, op::Char, delim::Char, cl::Char, delim_one::Bool, i1::Int64, n::Int64) @ Base .\show.jl:1244 [14] show_delim_array @ .\show.jl:1229 [inlined] [15] show @ .\show.jl:1262 [inlined] [16] show @ .\multimedia.jl:47 [inlined] [17] limitstringmime(mime::MIME{Symbol("text/plain")}, x::Tuple{Array{Float64, 3}, DIVAnd.DIVAnd_struct{Float64, Int64, 3, SparseMatrixCSC{Float64, Int64}}}) @ IJulia C:\Users\jmbeckers.julia\packages\IJulia\e8kqU\src\inline.jl:43 [18] display_mimestring @ C:\Users\jmbeckers.julia\packages\IJulia\e8kqU\src\display.jl:71 [inlined] [19] display_dict(x::Tuple{Array{Float64, 3}, DIVAnd.DIVAnd_struct{Float64, Int64, 3, SparseMatrixCSC{Float64, Int64}}}) @ IJulia C:\Users\jmbeckers.julia\packages\IJulia\e8kqU\src\display.jl:102 [20] #invokelatest#2 @ .\essentials.jl:716 [inlined] [21] invokelatest @ .\essentials.jl:714 [inlined] [22] execute_request(socket::ZMQ.Socket, msg::IJulia.Msg) @ IJulia C:\Users\jmbeckers.julia\packages\IJulia\e8kqU\src\execute_request.jl:112 [23] #invokelatest#2 @ .\essentials.jl:716 [inlined] [24] invokelatest @ .\essentials.jl:714 [inlined] [25] eventloop(socket::ZMQ.Socket) @ IJulia C:\Users\jmbeckers.julia\packages\IJulia\e8kqU\src\eventloop.jl:8 [26] (::IJulia.var"#15#18")() @ IJulia .\task.jl:423

jmbeckers commented 11 months ago

Actually, if you suppress the output on screen by adding a ";" at the end, the problem disappears. It seems it is the screen showing that causes the problem.

ctroupin commented 11 months ago

Same here, indeed it seems the display saturates the memory in gnome-terminal