I am using Preconditioners.jl to precondition the eigensolver lobpcg in IterativeSolvers.jl. But the preconditioners didn’t worked for more than one vectors. Any ideas?
This works when I need only one eigenpair.
using LinearAlgebra, IterativeSolvers, Preconditioners, SparseArrays
A = sprand(1000, 1000, 0.01)
A = A + A' + 30I
X0 = randn(1000,1)
tol = 1e-2
@time F = lobpcg(A, false, X0, tol=tol, maxiter=200, P=AMGPreconditioner{SmoothedAggregation}(A))
Output :
4.468253 seconds (34.34 M allocations: 1.563 GiB, 12.21% gc time, 99.73% compilation time)
However, when it came to more than one eigenpairs, it failed:
using LinearAlgebra, IterativeSolvers, Preconditioners, SparseArrays
A = sprand(1000, 1000, 0.01)
A = A + A' + 30I
X0 = randn(1000,2)
tol = 1e-2
@time F = lobpcg(A, false, X0, tol=tol, maxiter=200, P=AMGPreconditioner{SmoothedAggregation}(A))
Errors:
ERROR: LoadError: DimensionMismatch:
Stacktrace:
[1] mul!(C::Vector{Float64}, A::SparseMatrixCSC{Float64, Int64}, B::SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}, α::Bool, β::Bool)
@ SparseArrays /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:33
[2] mul!
@ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/LinearAlgebra/src/matmul.jl:276 [inlined]
[3] __solve!(x::SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}, ml::AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Adjoint{Float64, SparseMatrixCSC{Float64, Int64}}, AlgebraicMultigrid.MultiLevelWorkspace{Vector{Float64}, 1}}, cycle::AlgebraicMultigrid.V, b::SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}, lvl::Int64)
@ AlgebraicMultigrid ~/.julia/packages/AlgebraicMultigrid/ASpK7/src/multilevel.jl:216
[4] _solve!(x::SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}, ml::AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Adjoint{Float64, SparseMatrixCSC{Float64, Int64}}, AlgebraicMultigrid.MultiLevelWorkspace{Vector{Float64}, 1}}, b::SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}, cycle::AlgebraicMultigrid.V; maxiter::Int64, abstol::Float64, reltol::Float64, verbose::Bool, log::Bool, calculate_residual::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ AlgebraicMultigrid ~/.julia/packages/AlgebraicMultigrid/ASpK7/src/multilevel.jl:179
[5] ldiv!(x::SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}, p::AlgebraicMultigrid.Preconditioner{AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Adjoint{Float64, SparseMatrixCSC{Float64, Int64}}, AlgebraicMultigrid.MultiLevelWorkspace{Vector{Float64}, 1}}, AlgebraicMultigrid.V}, b::SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true})
@ AlgebraicMultigrid ~/.julia/packages/AlgebraicMultigrid/ASpK7/src/preconditioner.jl:18
[6] ldiv!
@ ~/.julia/packages/Preconditioners/CVS1Y/src/amg.jl:35 [inlined]
[7] (::IterativeSolvers.RPreconditioner{AMGPreconditioner{SmoothedAggregation, AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Adjoint{Float64, SparseMatrixCSC{Float64, Int64}}, AlgebraicMultigrid.MultiLevelWorkspace{Vector{Float64}, 1}}, AlgebraicMultigrid.V}, Float64, Matrix{Float64}})(X::SubArray{Float64, 2, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true})
@ IterativeSolvers ~/.julia/packages/IterativeSolvers/rhYBz/src/lobpcg.jl:238
[8] precond_constr!(block::Matrix{Float64}, temp_block::Matrix{Float64}, bs::Int64, precond!::IterativeSolvers.RPreconditioner{AMGPreconditioner{SmoothedAggregation, AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Adjoint{Float64, SparseMatrixCSC{Float64, Int64}}, AlgebraicMultigrid.MultiLevelWorkspace{Vector{Float64}, 1}}, AlgebraicMultigrid.V}, Float64, Matrix{Float64}}, constr!::IterativeSolvers.Constraint{Nothing, Matrix{Nothing}, Matrix{Nothing}, Nothing})
@ IterativeSolvers ~/.julia/packages/IterativeSolvers/rhYBz/src/lobpcg.jl:565
[9] (::LOBPCGIterator{false, Float64, SparseMatrixCSC{Float64, Int64}, Nothing, Vector{Float64}, Vector{Float64}, Vector{Int64}, Matrix{Float64}, IterativeSolvers.Blocks{false, Float64, Matrix{Float64}}, IterativeSolvers.CholQR{Matrix{Float64}}, IterativeSolvers.RPreconditioner{AMGPreconditioner{SmoothedAggregation, AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Adjoint{Float64, SparseMatrixCSC{Float64, Int64}}, AlgebraicMultigrid.MultiLevelWorkspace{Vector{Float64}, 1}}, AlgebraicMultigrid.V}, Float64, Matrix{Float64}}, IterativeSolvers.Constraint{Nothing, Matrix{Nothing}, Matrix{Nothing}, Nothing}, IterativeSolvers.BlockGram{false, Matrix{Float64}}, Vector{Bool}, Vector{IterativeSolvers.LOBPCGState{Vector{Float64}, Vector{Float64}}}})(residualTolerance::Float64, log::Bool)
@ IterativeSolvers ~/.julia/packages/IterativeSolvers/rhYBz/src/lobpcg.jl:709
[10] lobpcg!(iterator::LOBPCGIterator{false, Float64, SparseMatrixCSC{Float64, Int64}, Nothing, Vector{Float64}, Vector{Float64}, Vector{Int64}, Matrix{Float64}, IterativeSolvers.Blocks{false, Float64, Matrix{Float64}}, IterativeSolvers.CholQR{Matrix{Float64}}, IterativeSolvers.RPreconditioner{AMGPreconditioner{SmoothedAggregation, AlgebraicMultigrid.MultiLevel{AlgebraicMultigrid.Pinv{Float64}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, AlgebraicMultigrid.GaussSeidel{AlgebraicMultigrid.SymmetricSweep}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Adjoint{Float64, SparseMatrixCSC{Float64, Int64}}, AlgebraicMultigrid.MultiLevelWorkspace{Vector{Float64}, 1}}, AlgebraicMultigrid.V}, Float64, Matrix{Float64}}, IterativeSolvers.Constraint{Nothing, Matrix{Nothing}, Matrix{Nothing}, Nothing}, IterativeSolvers.BlockGram{false, Matrix{Float64}}, Vector{Bool}, Vector{IterativeSolvers.LOBPCGState{Vector{Float64}, Vector{Float64}}}}; log::Bool, maxiter::Int64, not_zeros::Bool, tol::Float64)
@ IterativeSolvers ~/.julia/packages/IterativeSolvers/rhYBz/src/lobpcg.jl:881
[11] #lobpcg#56
@ ~/.julia/packages/IterativeSolvers/rhYBz/src/lobpcg.jl:838 [inlined]
[12] #lobpcg#55
@ ~/.julia/packages/IterativeSolvers/rhYBz/src/lobpcg.jl:825 [inlined]
[13] top-level scope
@ ./timing.jl:262
in expression starting at /Users/neo/Desktop/work/oracle/gd/test_lobpcg_amg.jl:44