search

JuliaLinearAlgebra / IterativeSolvers.jl

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

support matrix RHSs #248

Open chriscoey opened 5 years ago

chriscoey commented 5 years ago

I need to solve linear systems A * x = b where x and b have 2 or 3 columns (the application is for interior point optimization routines: specifically finding prediction and correction directions simultaneously). A is a BlockArray with structured and sparse blocks that I do not want to materialize and factorize due to memory constraints. Currently if I try to use different iterative solvers from this package, I get various errors such as BoundsError when b has more than 1 column.

mohamed82008 commented 5 years ago

This works but the allocating version errors, for a trivially fixable reason.

using LinearAlgebra, IterativeSolvers

cg!(rand(2,2), Diagonal(rand(2)), rand(2,2))

In https://github.com/JuliaMath/IterativeSolvers.jl/blob/master/src/common.jl#L20, size(A, 2) should be size(b).

mohamed82008 commented 5 years ago

Or actually the size is completely unnecessary for cg. Not sure if there is another solver that uses this function differently.

chriscoey commented 5 years ago

OK the in-place cg! works, but gmres! (and lsmr!, lsqr! and others) errors:

julia> gmres!(x, B, b)
ERROR: BoundsError: attempt to access 4-element view(::Array{Float64,2}, :, 1) with eltype Float64 at index [1, 2, 3, 4, 5, 6, 7, 8]
Stacktrace:
 [1] copyto!(::IndexLinear, ::SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true}, ::IndexLinear, ::Array{Float64,2}) at ./abstractarray.jl:804
 [2] copyto! at ./abstractarray.jl:799 [inlined]
 [3] #init!#33(::Bool, ::typeof(IterativeSolvers.init!), ::IterativeSolvers.ArnoldiDecomp{Float64,Array{Float64,2}}, ::Array{Float64,2}, ::Array{Float64,2}, ::Identity, ::Array{Float64,2}) at /home/coey/.julia/packages/IterativeSolvers/QuajJ/src/gmres.jl:220
 [4] #init! at ./none:0 [inlined]
 [5] #gmres_iterable!#30(::Identity, ::Identity, ::Float64, ::Int64, ::Int64, ::Bool, ::typeof(IterativeSolvers.gmres_iterable!), ::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}) at /home/coey/.julia/packages/IterativeSolvers/QuajJ/src/gmres.jl:119
 [6] #gmres_iterable! at ./none:0 [inlined]
 [7] #gmres!#32(::Identity, ::Identity, ::Float64, ::Int64, ::Int64, ::Bool, ::Bool, ::Bool, ::typeof(gmres!), ::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}) at /home/coey/.julia/packages/IterativeSolvers/QuajJ/src/gmres.jl:186
 [8] gmres!(::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}) at /home/coey/.julia/packages/IterativeSolvers/QuajJ/src/gmres.jl:182
 [9] top-level scope at REPL[137]:1