cg broken?

[madeleineudell]

On Julia 1.0, I get bounds errors from running cg. Here's a simple test:

using LinearAlgebra
import IterativeSolvers: cg, cg!

m,n = 10,5
A = randn(m,n)
x = randn(n)
b = A*x + .1*randn(m)
A\b                       # works

xhat = zeros(size(x))
cg(A, b)                  # bounds error
cg!(xhat, A, b)           # bounds error

the calls to cg give the following errors:

julia> cg(A, b)
ERROR: BoundsError
Stacktrace:
 [1] copyto! at ./array.jl:272 [inlined]
 [2] copyto! at ./array.jl:277 [inlined]
 [3] #cg_iterator!#21(::Float64, ::Int64, ::IterativeSolvers.CGStateVariables{Float64,Array{Float64,1}}, ::Bool, ::Function, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}, ::IterativeSolvers.Identity) at /Users/madeleine/.julia/packages/IterativeSolvers/O0pSc/src/cg.jl:126
 [4] #cg_iterator! at ./none:0 [inlined]
 [5] #cg!#23(::Float64, ::Int64, ::Bool, ::IterativeSolvers.CGStateVariables{Float64,Array{Float64,1}}, ::Bool, ::IterativeSolvers.Identity, ::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:initially_zero,),Tuple{Bool}}}, ::Function, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}) at /Users/madeleine/.julia/packages/IterativeSolvers/O0pSc/src/cg.jl:215
 [6] (::getfield(IterativeSolvers, Symbol("#kw##cg!")))(::NamedTuple{(:initially_zero,),Tuple{Bool}}, ::typeof(cg!), ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}) at ./none:0
 [7] cg(::Array{Float64,2}, ::Array{Float64,1}) at /Users/madeleine/.julia/packages/IterativeSolvers/O0pSc/src/cg.jl:161
 [8] top-level scope at none:0

julia> cg!(xhat, A, b)
ERROR: BoundsError
Stacktrace:
 [1] copyto! at ./array.jl:272 [inlined]
 [2] copyto! at ./array.jl:277 [inlined]
 [3] #cg_iterator!#21(::Float64, ::Int64, ::IterativeSolvers.CGStateVariables{Float64,Array{Float64,1}}, ::Bool, ::Function, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}, ::IterativeSolvers.Identity) at /Users/madeleine/.julia/packages/IterativeSolvers/O0pSc/src/cg.jl:126
 [4] #cg_iterator! at ./none:0 [inlined]
 [5] #cg!#23(::Float64, ::Int64, ::Bool, ::IterativeSolvers.CGStateVariables{Float64,Array{Float64,1}}, ::Bool, ::IterativeSolvers.Identity, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}) at /Users/madeleine/.julia/packages/IterativeSolvers/O0pSc/src/cg.jl:215
 [6] cg!(::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}) at /Users/madeleine/.julia/packages/IterativeSolvers/O0pSc/src/cg.jl:210
 [7] top-level scope at none:0

[madeleineudell]

Never mind! Of course A should be square (and positive definite).

[haampie]

Hm, yeah, the package is not really giving helpful feedback for incorrect input, maybe a size check is worth it!

[dkarrasch]

👍 , a quick LinearAlgebra.checksquare(A) would do the job. That doesn't cost much, and should preferably be included in all methods that require square matrices/linear maps.