Issues with interactions of LOBPCG and LinearMaps.LinearCombination

[tpolakovic]

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]

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]

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

[mohamed82008]

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]

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]

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

[dkarrasch]

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.