Plot of the block-structure gives error

[chaozg]

Hi there,

It seems the plot(H) gives an error for even a small problem. To reproduce the error, I used the code from README.md and changed m to a very small integer, like 600.

using HMatrices, LinearAlgebra, StaticArrays
const Point3D = SVector{3,Float64}
# sample some points on a sphere
m = 600 #500 is fine but it's dense so nothing is interesting
X = Y = [Point3D(sin(θ)cos(ϕ),sin(θ)*sin(ϕ),cos(θ)) for (θ,ϕ) in zip(π*rand(m),2π*rand(m))]
function G(x,y) 
  d = norm(x-y) + 1e-8
  1/(4π*d)
end
K = KernelMatrix(G,X,Y)

H = assemble_hmatrix(K;atol=1e-6)

using Plots
plot(H)

And the error is as follows:

ERROR: method `getindex(::HMatrix,args...)` has been disabled to
avoid performance pitfalls. Unless you made an explicit call to `getindex`,
this error usually means that a linear algebra routine involving an
`HMatrix` has fallen back to a generic implementation.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] getindex(::HMatrix{ClusterTree{3, Float64}, Float64}, i::Int64, j::Int64)
    @ HMatrices ~/.julia/packages/HMatrices/OBct7/src/hmatrix.jl:87
  [3] _getindex
    @ ./abstractarray.jl:1341 [inlined]
  [4] getindex
    @ ./abstractarray.jl:1291 [inlined]
  [5] iterate
    @ ./abstractarray.jl:1217 [inlined]
  [6] iterate
    @ ./abstractarray.jl:1215 [inlined]
  [7] copyto_unaliased!
    @ ./abstractarray.jl:1093 [inlined]
  [8] copyto!(dest::Matrix{Float64}, src::HMatrix{ClusterTree{3, Float64}, Float64})
    @ Base ./abstractarray.jl:1068
  [9] copy_similar
    @ ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/LinearAlgebra.jl:416 [inlined]
 [10] eigencopy_oftype
    @ ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/eigen.jl:538 [inlined]
 [11] svdvals(A::HMatrix{ClusterTree{3, Float64}, Float64})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/svd.jl:242
 [12] rank(A::HMatrix{ClusterTree{3, Float64}, Float64}; atol::Float64, rtol::Float64)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/generic.jl:991
 [13] rank(A::HMatrix{ClusterTree{3, Float64}, Float64})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/generic.jl:989
 [14] macro expansion
    @ ~/.julia/packages/HMatrices/OBct7/src/hmatrix.jl:485 [inlined]
 [15] macro expansion
    @ ~/.julia/packages/RecipesBase/BRe07/src/RecipesBase.jl:340 [inlined]
 [16] macro expansion
    @ ~/.julia/packages/HMatrices/OBct7/src/hmatrix.jl:482 [inlined]
 [17] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, hmat::HMatrix)
    @ HMatrices ~/.julia/packages/RecipesBase/BRe07/src/RecipesBase.jl:300
 [18] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/user_recipe.jl:38
 [19] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/RecipesPipeline.jl:72
 [20] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
    @ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:223
 [21] #plot#188
    @ ~/.julia/packages/Plots/sxUvK/src/plot.jl:102 [inlined]
 [22] plot(args::Any)
    @ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:93
 [23] top-level scope
    @ REPL[21]:1

[maltezfaria]

Thanks for raising the issue! At some point I must have removed rank(hmatrix), which broke the plotting functionality. This should be fixed in main now. Let me know if that is not the case.

[chaozg]

Thanks, @maltezfaria. I just tested the code on main and everything is good now. Thanks!

[chaozg]

Just noticed that plot still throws an error for Hermitian HMatrix, for example with the following snippet:

using HMatrices, LinearAlgebra, StaticArrays, Random
using Plots
# copied from README.md
const Point3D = SVector{3,Float64}
m = 800
X = Y = [Point3D(sin(θ)cos(ϕ),sin(θ)*sin(ϕ),cos(θ)) for (θ,ϕ) in zip(π*rand(m),2π*rand(m))]
function G(x,y) 
  d = norm(x-y) + 1e-8
  exp(-d)
end
K = KernelMatrix(G,X,Y)

# copied from test/cholesky_test.jl
splitter = CardinalitySplitter(; nmax = 50)
Xclt = ClusterTree(X, splitter)
Yclt = ClusterTree(Y, splitter)
adm = StrongAdmissibilityStd(10)
comp = PartialACA(; atol = 1e-6)

H1 = assemble_hmatrix(K; atol=1e-6)
H2 = assemble_hmatrix(Hermitian(K), Xclt, Yclt; adm, comp, threads=false, distributed = false)
plot(H1) # this line runs fine
plot(H2) # this line throws an error

ERROR: method `getindex(::HMatrix,args...)` has been disabled to
avoid performance pitfalls. Unless you made an explicit call to `getindex`,
this error usually means that a linear algebra routine involving an
`HMatrix` has fallen back to a generic implementation.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] getindex(::Hermitian{Float64, HMatrix{ClusterTree{3, Float64}, Float64}}, i::Int64, j::Int64)
    @ HMatrices ~/.julia/packages/HMatrices/ZULqq/src/hmatrix.jl:87
  [3] macro expansion
    @ ./multidimensional.jl:917 [inlined]
  [4] macro expansion
    @ ./cartesian.jl:64 [inlined]
  [5] _unsafe_getindex!
    @ ./multidimensional.jl:912 [inlined]
  [6] _unsafe_getindex
    @ ./multidimensional.jl:903 [inlined]
  [7] _getindex
    @ ./multidimensional.jl:889 [inlined]
  [8] getindex
    @ ./abstractarray.jl:1291 [inlined]
  [9] #25
    @ ./none:0 [inlined]
 [10] iterate
    @ ./generator.jl:47 [inlined]
 [11] collect(itr::Base.Generator{Base.OneTo{Int64}, RecipesPipeline.var"#25#26"{Hermitian{Float64, HMatrix{…}}}})
    @ Base ./array.jl:834
 [12] _series_data_vector(v::Hermitian{Float64, HMatrix{ClusterTree{3, Float64}, Float64}}, plotattributes::Dict{Symbol, Any})
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/series.jl:60
 [13] macro expansion
    @ ~/.julia/packages/RecipesPipeline/BGM3l/src/series.jl:129 [inlined]
 [14] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, ::Type{RecipesPipeline.SliceIt}, x::Any, y::Any, z::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesBase/BRe07/src/RecipesBase.jl:300
 [15] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/user_recipe.jl:38
 [16] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/RecipesPipeline.jl:72
 [17] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
    @ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:223
 [18] #plot#188
    @ ~/.julia/packages/Plots/sxUvK/src/plot.jl:102 [inlined]
 [19] plot(args::Any)
    @ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:93
 [20] top-level scope
    @ REPL[31]:1
Some type information was truncated. Use `show(err)` to see complete types.

[maltezfaria]

Should be patched now with #68