search

JuliaLinearAlgebra / IterativeSolvers.jl

Iterative algorithms for solving linear systems, eigensystems, and singular value problems
MIT License
399 stars 105 forks source link

Issues with interactions of LOBPCG and LinearMaps.LinearCombination #259

Closed tpolakovic closed 4 years ago

tpolakovic commented 4 years ago

I'm not entirely sure whether this issue should be reported at LinearMaps or here, but here goes the MWE (on Julia 1.3):

> a = LinearMap(identity, 10)
> b = a + a
> lobpcg(b, true, 3)

ERROR: MethodError: no method matching A_mul_B!(::Array{Float64,2}, ::LinearMaps.FunctionMap{Float64,typeof(identity),Nothing}, ::Array{Float64,2})

Strangely enough, if I wrap the map, it works:

> c = LinearMap(b)
> lobpcg(c, true, 3)

Results of LOBPCG Algorithm
 * Algorithm: LOBPCG - CholQR
 * λ: [2.000000000000001,2.0000000000000004, ...]
 * Residual norm(s): [8.455847260422304e-16,4.335559509131367e-16, ...]
 * Convergence
   * Iterations: 1
   * Converged: true
   * Iterations limit: 200

The offender seems to be the Array{Float64,2}, where it tries to do multiplication with a row vector instead a column vector Array{Float64,1} but I'm not sure why that is the case.

dkarrasch commented 4 years ago

Can you show the full stacktrace? What is calling A_mul_B!? Because mul!tiplying a FunctionMap with a matrix works:

julia> using LinearMaps
[ Info: Precompiling LinearMaps [7a12625a-238d-50fd-b39a-03d52299707e]

julia> a = LinearMap(identity, 10)
LinearMaps.FunctionMap{Float64}(identity, 10, 10; ismutating=false, issymmetric=false, ishermitian=false, isposdef=false)

julia> b = a + a
LinearMaps.LinearCombination{Float64,Tuple{LinearMaps.FunctionMap{Float64,typeof(identity),Nothing},LinearMaps.FunctionMap{Float64,typeof(identity),Nothing}}}((LinearMaps.FunctionMap{Float64}(identity, 10, 10; ismutating=false, issymmetric=false, ishermitian=false, isposdef=false), LinearMaps.FunctionMap{Float64}(identity, 10, 10; ismutating=false, issymmetric=false, ishermitian=false, isposdef=false)))

julia> using LinearAlgebra

julia> mul!(zeros(10,3), b, rand(10, 3))
10×3 Array{Float64,2}:
 1.77102   1.00707   1.0886  
 0.718971  0.853393  0.156217
 0.49425   0.908157  1.91741 
 0.928176  0.230759  0.580476
 0.818615  0.511816  0.195229
 1.52616   1.5837    0.964583
 1.16619   1.55489   1.22776 
 1.17506   0.316413  0.779997
 0.28235   0.558206  1.7167  
 0.200653  0.873268  0.634083
mohamed82008 commented 4 years ago

Please report the versions of the packages using ]st LinearMaps and ]st IterativeSolvers.

mohamed82008 commented 4 years ago

I suspect the new Pkg upper bounds got you an ancient version of IterativeSolvers that only supports Julia 0.6 (but claims to support everything).

tpolakovic commented 4 years ago

I'm using Julia 1.3 (2019-11-26) LinearMaps is pulled from the most recent commit because I need the kron functionality.

> st LinearMaps
    Status `~/.julia/environments/v1.3/Project.toml`
  [7a12625a] LinearMaps v2.5.2 #d9ca214 (https://github.com/Jutho/LinearMaps.jl.git)
> st IterativeSolvers
    Status `~/.julia/environments/v1.3/Project.toml`
  [42fd0dbc] IterativeSolvers v0.8.1

Here's the stack trace:

ERROR: MethodError: no method matching A_mul_B!(::Array{Float64,2}, ::LinearMaps.FunctionMap{Float64,typeof(identity),Nothing}, ::Array{Float64,2})
Closest candidates are:
  A_mul_B!(::AbstractArray{T,1} where T, ::LinearMaps.FunctionMap, ::AbstractArray{T,1} where T) at /Users/bubu/.julia/packages/LinearMaps/njgRB/src/functionmap.jl:53
Stacktrace:
 [1] _mul! at /Users/bubu/.julia/packages/LinearMaps/njgRB/src/linearcombination.jl:99 [inlined]
 [2] mul! at /Users/bubu/.julia/packages/LinearMaps/njgRB/src/linearcombination.jl:79 [inlined]
 [3] mul! at /Users/bubu/.julia/packages/LinearMaps/njgRB/src/linearcombination.jl:70 [inlined]
 [4] A_mul_X!(::IterativeSolvers.Blocks{false,Float64,Array{Float64,2}}, ::LinearMaps.LinearCombination{Float64,Tuple{LinearMaps.FunctionMap{Float64,typeof(identity),Nothing},LinearMaps.FunctionMap{Float64,typeof(identity),Nothing}}}) at /Users/bubu/.julia/packages/IterativeSolvers/QuajJ/src/lobpcg.jl:125
 [5] ortho_AB_mul_X!(::IterativeSolvers.Blocks{false,Float64,Array{Float64,2}}, ::IterativeSolvers.CholQR{Array{Float64,2}}, ::LinearMaps.LinearCombination{Float64,Tuple{LinearMaps.FunctionMap{Float64,typeof(identity),Nothing},LinearMaps.FunctionMap{Float64,typeof(identity),Nothing}}}, ::Nothing, ::Int64) at /Users/bubu/.julia/packages/IterativeSolvers/QuajJ/src/lobpcg.jl:530
 [6] ortho_AB_mul_X!(::IterativeSolvers.Blocks{false,Float64,Array{Float64,2}}, ::IterativeSolvers.CholQR{Array{Float64,2}}, ::LinearMaps.LinearCombination{Float64,Tuple{LinearMaps.FunctionMap{Float64,typeof(identity),Nothing},LinearMaps.FunctionMap{Float64,typeof(identity),Nothing}}}, ::Nothing) at /Users/bubu/.julia/packages/IterativeSolvers/QuajJ/src/lobpcg.jl:526
 [7] (::LOBPCGIterator{false,Float64,LinearMaps.LinearCombination{Float64,Tuple{LinearMaps.FunctionMap{Float64,typeof(identity),Nothing},LinearMaps.FunctionMap{Float64,typeof(identity),Nothing}}},Nothing,Array{Float64,1},Array{Float64,1},Array{Int64,1},Array{Float64,2},IterativeSolvers.Blocks{false,Float64,Array{Float64,2}},IterativeSolvers.CholQR{Array{Float64,2}},IterativeSolvers.RPreconditioner{Nothing,Float64,Array{Float64,2}},IterativeSolvers.Constraint{Nothing,Array{Nothing,2},Array{Nothing,2},Nothing},IterativeSolvers.BlockGram{false,Array{Float64,2}},Array{Bool,1},Array{IterativeSolvers.LOBPCGState{Array{Float64,1},Array{Float64,1}},1}})(::Float64, ::Bool) at /Users/bubu/.julia/packages/IterativeSolvers/QuajJ/src/lobpcg.jl:696
 [8] #lobpcg!#53(::Bool, ::Int64, ::Bool, ::Float64, ::typeof(lobpcg!), ::LOBPCGIterator{false,Float64,LinearMaps.LinearCombination{Float64,Tuple{LinearMaps.FunctionMap{Float64,typeof(identity),Nothing},LinearMaps.FunctionMap{Float64,typeof(identity),Nothing}}},Nothing,Array{Float64,1},Array{Float64,1},Array{Int64,1},Array{Float64,2},IterativeSolvers.Blocks{false,Float64,Array{Float64,2}},IterativeSolvers.CholQR{Array{Float64,2}},IterativeSolvers.RPreconditioner{Nothing,Float64,Array{Float64,2}},IterativeSolvers.Constraint{Nothing,Array{Nothing,2},Array{Nothing,2},Nothing},IterativeSolvers.BlockGram{false,Array{Float64,2}},Array{Bool,1},Array{IterativeSolvers.LOBPCGState{Array{Float64,1},Array{Float64,1}},1}}) at /Users/bubu/.julia/packages/IterativeSolvers/QuajJ/src/lobpcg.jl:880
 [9] #lobpcg! at ./none:0 [inlined]
 [10] #lobpcg#52 at /Users/bubu/.julia/packages/IterativeSolvers/QuajJ/src/lobpcg.jl:837 [inlined]
 [11] #lobpcg at ./none:0 [inlined]
 [12] #lobpcg#50(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(lobpcg), ::LinearMaps.LinearCombination{Float64,Tuple{LinearMaps.FunctionMap{Float64,typeof(identity),Nothing},LinearMaps.FunctionMap{Float64,typeof(identity),Nothing}}}, ::Nothing, ::Bool, ::Int64) at /Users/bubu/.julia/packages/IterativeSolvers/QuajJ/src/lobpcg.jl:791
 [13] lobpcg at /Users/bubu/.julia/packages/IterativeSolvers/QuajJ/src/lobpcg.jl:791 [inlined]
 [14] #lobpcg#49 at /Users/bubu/.julia/packages/IterativeSolvers/QuajJ/src/lobpcg.jl:788 [inlined]
 [15] lobpcg(::LinearMaps.LinearCombination{Float64,Tuple{LinearMaps.FunctionMap{Float64,typeof(identity),Nothing},LinearMaps.FunctionMap{Float64,typeof(identity),Nothing}}}, ::Bool, ::Int64) at /Users/bubu/.julia/packages/IterativeSolvers/QuajJ/src/lobpcg.jl:788
 [16] top-level scope at REPL[5]:1
tpolakovic commented 4 years ago

I tried using LinearMaps from the registry (2.5.2) and it works with there. So there's an issue on that side.

dkarrasch commented 4 years ago

I'm glad somebody is interested in/using the kron functionality. I can confirm that this is on LinearMaps.jl side. It is fixed very soon.