search

JuliaLinearAlgebra / IterativeSolvers.jl

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

Some issues with non-determinancy of output #234

Closed pkofod closed 5 years ago

pkofod commented 5 years ago

Sorry if I'm not fully understanding these things, but is the behavior of gmres! and idrs! not meant to be deterministic? Please consider the following

julia> J = [-1.0 0.0; 24.0 10.0]
2×2 Array{Float64,2}:
 -1.0   0.0
 24.0  10.0

julia> v = [2.2, -4.4]
2-element Array{Float64,1}:
  2.2
 -4.4

julia> x = similar(v)
2-element Array{Float64,1}:
 6.9189794760214e-310
 6.9189606358082e-310

julia> idrs!(x, J, v)
2-element Array{Float64,1}:
 -2.2 
  4.84

julia> idrs!(x, J, v)
ERROR: UndefVarError: res not defined
Stacktrace:
 [1] #idrs_method!#39(::Bool, ::Bool, ::Function, ::ConvergenceHistory{false,Nothing}, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}, ::Int64, ::Float64, ::Int64) at /home/pkofod/.julia/packages/IterativeSolvers/QuajJ/src/idrs.jl:77
 [2] idrs_method!(::ConvergenceHistory{false,Nothing}, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}, ::Int64, ::Float64, ::Int64) at /home/pkofod/.julia/packages/IterativeSolvers/QuajJ/src/idrs.jl:76
 [3] #idrs!#38(::Int64, ::Float64, ::Int64, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}) at /home/pkofod/.julia/packages/IterativeSolvers/QuajJ/src/idrs.jl:50
 [4] idrs!(::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}) at /home/pkofod/.julia/packages/IterativeSolvers/QuajJ/src/idrs.jl:47
 [5] top-level scope at none:0

julia> gmres!(x, J', v)
2-element Array{Float64,1}:
 -12.759999999999987 
  -0.4399999999999977

julia> gmres!(x, J', v)
2-element Array{Float64,1}:
 -12.760000000000002  
  -0.44000000000000006

julia> gmres!(x, J, v)
2-element Array{Float64,1}:
 -2.200000000000003
  4.840000000000001

julia> gmres!(x, J, v)
2-element Array{Float64,1}:
 -2.2              
  4.840000000000001

related to #222

andreasnoack commented 5 years ago

This is expected since x is the initial guess. E.g.

julia> x = rand(length(v))
2-element Array{Float64,1}:
 0.0009348838551288541
 0.5852360604763667

julia> gmres!(copy(x), J', v)
2-element Array{Float64,1}:
 -12.760000000000023
  -0.4400000000000017

julia> gmres!(copy(x), J', v)
2-element Array{Float64,1}:
 -12.760000000000023
  -0.4400000000000017

julia> gmres!(copy(x), J', v)
2-element Array{Float64,1}:
 -12.760000000000023
  -0.4400000000000017
pkofod commented 5 years ago

Right, so I'm apparently too stupid to read a docstring :) I thought it was just a cache.