ArnoldiMethod: strange error

General Usage
1

For those who are familiar with ArnoldiMethod.jl:

This works

function (SI::ShiftAndInvert)(y, x)
    mul!(SI.temp, SI.B, x)
    y .= SI.A_factorization \ SI.temp

    # x -> PtL \ (B * (PtL' \ x)),
    # y .= SI.A_factorization.PtL' \ x
    # mul!(SI.temp, SI.B, y)
    # y .= SI.A_factorization.PtL \ SI.temp
end

And this does not

function (SI::ShiftAndInvert)(y, x)
    # mul!(SI.temp, SI.B, x)
    # y .= SI.A_factorization \ SI.temp

    # x -> PtL \ (B * (PtL' \ x)),
    y .= SI.A_factorization.PtL' \ x
    mul!(SI.temp, SI.B, y)
    y .= SI.A_factorization.PtL \ SI.temp
end

It gives this error:

ERROR: LoadError: CanonicalIndexError: getindex not defined for SparseArrays.CHOLMOD.FactorComponent{Float64, :UP, Int64}                                                                                        
Stacktrace:                                                                                                                                                                                                      
  [1] error_if_canonical_getindex(::IndexCartesian, ::SparseArrays.CHOLMOD.FactorComponent{Float64, :UP, Int64}, ::Int64, ::Int64)                                                                               
    @ Base .\abstractarray.jl:1307                                                                                                                                                                               
  [2] getindex                                                                                                                                                                                                   
    @ .\abstractarray.jl:1287 [inlined]                                                                                                                                                                          
  [3] getindex                                                                                                                                                                                                   
    @ .\subarray.jl:288 [inlined]                                                                                                                                                                                
  [4] iterate                                                                                                                                                                                                    
    @ .\abstractarray.jl:1214 [inlined]                                                                                                                                                                          
  [5] iterate                                                                                                                                                                                                    
    @ .\abstractarray.jl:1212 [inlined]                                                                                                                                                                          
  [6] _all                                                                                                                                                                                                       
    @ .\reduce.jl:1287 [inlined]                                                                                                                                                                                 
  [7] all                                                                                                                                                                                                        
    @ .\reducedim.jl:1023 [inlined]                                                                                                                                                                              
  [8] _istril                                                                                                                                                                                                    
    @ C:\Users\pkonl\SublimeJulia_1_10_0\assets\julia-1.10.0\share\julia\stdlib\v1.10\LinearAlgebra\src\generic.jl:1349 [inlined]                                                                                
  [9] istril(A::SparseArrays.CHOLMOD.FactorComponent{Float64, :UP, Int64}, k::Int64)                                                                                                                             
    @ LinearAlgebra C:\Users\pkonl\SublimeJulia_1_10_0\assets\julia-1.10.0\share\julia\stdlib\v1.10\LinearAlgebra\src\generic.jl:1342                                                                            
 [10] istril(A::SparseArrays.CHOLMOD.FactorComponent{Float64, :UP, Int64}, k::Int64)                                                                                                                             
    @ LinearAlgebra C:\Users\pkonl\SublimeJulia_1_10_0\assets\julia-1.10.0\share\julia\stdlib\v1.10\LinearAlgebra\src\generic.jl:1341 [inlined]                                                                  
 [11] \(A::SparseArrays.CHOLMOD.FactorComponent{Float64, :UP, Int64}, B::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})                                               
    @ LinearAlgebra C:\Users\pkonl\SublimeJulia_1_10_0\assets\julia-1.10.0\share\julia\stdlib\v1.10\LinearAlgebra\src\generic.jl:1114                                                                            
 [12] (::GEPHelpers.ShiftAndInvert{SparseArrays.CHOLMOD.Factor{Float64, Int64}, Symmetric{Float64, SparseMatrixCSC{Float64, Int64}}, Vector{Float64}})(y::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, x::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})                                                                                  
    @ GEPHelpers C:\Users\pkonl\Documents\00WIP\GEPHelpers.jl\src\gep_smallest.jl:82                                                                                                                             
 [13] _unsafe_mul!                                                                                                                                                                                               
    @ C:\Users\pkonl\SublimeJulia_1_10_0\assets\.julia-1.10.0-depot\packages\LinearMaps\6ej60\src\functionmap.jl:128 [inlined]                                                                                   
 [14] mul!                                                                                                                                                                                                       
    @ C:\Users\pkonl\SublimeJulia_1_10_0\assets\.julia-1.10.0-depot\packages\LinearMaps\6ej60\src\LinearMaps.jl:163 [inlined]                                                                                    
 [15] iterate_arnoldi!(A::LinearMaps.FunctionMap{Float64, GEPHelpers.ShiftAndInvert{SparseArrays.CHOLMOD.Factor{Float64, Int64}, Symmetric{Float64, SparseMatrixCSC{Float64, Int64}}, Vector{Float64}}, Nothing, 
true}, arnoldi::ArnoldiMethod.Arnoldi{Float64, Matrix{Float64}, Matrix{Float64}}, range::UnitRange{Int64})                                                                                                       
    @ ArnoldiMethod C:\Users\pkonl\SublimeJulia_1_10_0\assets\.julia-1.10.0-depot\packages\ArnoldiMethod\JdEiw\src\expansion.jl:117                                                                              
 [16] _partialschur(A::LinearMaps.FunctionMap{Float64, GEPHelpers.ShiftAndInvert{SparseArrays.CHOLMOD.Factor{Float64, Int64}, Symmetric{Float64, SparseMatrixCSC{Float64, Int64}}, Vector{Float64}}, Nothing, true}, ::Type{Float64}, mindim::Int64, maxdim::Int64, nev::Int64, tol::Float64, restarts::Int64, which::ArnoldiMethod.LM)                                                                                           
    @ ArnoldiMethod C:\Users\pkonl\SublimeJulia_1_10_0\assets\.julia-1.10.0-depot\packages\ArnoldiMethod\JdEiw\src\run.jl:185                                                                                    
 [17] partialschur(A::LinearMaps.FunctionMap{Float64, GEPHelpers.ShiftAndInvert{SparseArrays.CHOLMOD.Factor{Float64, Int64}, Symmetric{Float64, SparseMatrixCSC{Float64, Int64}}, Vector{Float64}}, Nothing, true}; nev::Int64, which::ArnoldiMethod.LM, tol::Float64, mindim::Int64, maxdim::Int64, restarts::Int64)                                                                                                             
    @ ArnoldiMethod C:\Users\pkonl\SublimeJulia_1_10_0\assets\.julia-1.10.0-depot\packages\ArnoldiMethod\JdEiw\src\run.jl:106                                                                                    
 [18] partialschur                                                                                                                                                                                               
    @ C:\Users\pkonl\SublimeJulia_1_10_0\assets\.julia-1.10.0-depot\packages\ArnoldiMethod\JdEiw\src\run.jl:94 [inlined]                                                                                         
 [19] __arnoldimethod_eigs(K::Symmetric{Float64, SparseMatrixCSC{Float64, Int64}}, M::Symmetric{Float64, SparseMatrixCSC{Float64, Int64}}; nev::Int64, ncv::Int64, tol::Float64, maxiter::Int64, v0::Vector{Float64}, check::Int64)                                                                                                                                                                                               
    @ GEPHelpers C:\Users\pkonl\Documents\00WIP\GEPHelpers.jl\src\gep_smallest.jl:103                                                                                                                            
 [20] __arnoldimethod_eigs                                                                                                                                                                                       
    @ C:\Users\pkonl\Documents\00WIP\GEPHelpers.jl\src\gep_smallest.jl:93 [inlined]

In both cases the linear map is created as

struct ShiftAndInvert{TA,TB,TT}
    A_factorization::TA
    B::TB
    temp::TT
end

function construct_linear_map(A, B)
    a = ShiftAndInvert(cholesky(A), B, Vector{eltype(A)}(undef, size(A, 1)))
    LinearMap{eltype(A)}(a, size(A, 1), ismutating = true)
end

I can’t understand why what is happening inside the function (SI::ShiftAndInvert)(y, x) should have this effect.

EDIT: In case you are wondering, I did try to construct the shift-and-invert object directly, and it works just fine:

    a = ShiftAndInvert(cholesky(K), M, Vector{eltype(K)}(undef, size(K, 1)))
    x = rand(size(K,1))
    y = rand(size(K,1))
    @show a(y, x)

produces the expected output.

2

You are trying to do F.PtL \ Av where Av is a view; the existing method is not designed for this (instead if falls back to a generic function which is incompatible with the sparse factorization). The method for F \ Av actually converts Av into a dense array object managed by the CHOLMOD library. You might file an issue against SparseArrays to see if the maintainers can do the same thing for F.PtL etc.

2 Likes
3

Thanks for the explanation, that makes sense. I did notice the view, but I did not expect the linear solved to have an issue with it!