## Test what happens with Neumann series and singular matrices.
using LinearAlgebra, StableRNGs
B = randn(StableRNG(1), 6,4)
A = B*B'
A = A/10 # scale it so that it's largest eigenvalue is less than 10.
6×6 Matrix{Float64}: 0.623896 0.0521145 0.0356488 0.207312 -0.232289 0.14814 0.0521145 0.0964917 -0.026238 -0.0920069 -0.0980198 0.100822 0.0356488 -0.026238 0.0949782 0.0306672 -0.203343 -0.0672769 0.207312 -0.0920069 0.0306672 0.305068 0.0449488 -0.105023 -0.232289 -0.0980198 -0.203343 0.0449488 0.707332 -0.00139221 0.14814 0.100822 -0.0672769 -0.105023 -0.00139221 0.180155
##
b = randn(StableRNG(2), 6 )
6-element Vector{Float64}: -1.6747889454750067 0.2157655329992774 -0.696915102409161 -0.8839710394954117 -0.4990442125777455 -0.3251864489603265
##
M = I - A
6×6 Matrix{Float64}: 0.376104 -0.0521145 -0.0356488 -0.207312 0.232289 -0.14814 -0.0521145 0.903508 0.026238 0.0920069 0.0980198 -0.100822 -0.0356488 0.026238 0.905022 -0.0306672 0.203343 0.0672769 -0.207312 0.0920069 -0.0306672 0.694932 -0.0449488 0.105023 0.232289 0.0980198 0.203343 -0.0449488 0.292668 0.00139221 -0.14814 -0.100822 0.0672769 0.105023 0.00139221 0.819845
##
b + M*b + M^2*b + M^3*b + M^4*b
6-element Vector{Float64}: -2.6445822485320494 0.5620602606990042 -3.5079755150325713 -0.9672166444826965 -2.0946144344078608 -0.9774379920802301
## run a lot of terms of the Neumann series
x = similar(b)
fill!(x, 0)
for i=1:100
x += M*x + b
end
x
6-element Vector{Float64}: 4.63502821849714e28 -3.1092581158593055e29 -7.068666726792267e29 -7.468305836244001e28 -2.2749694456352612e29 -1.7447404548345065e29
## Let's use a right hand side that is consistent.
xtrue = randn(StableRNG(3), 6)
b = A*xtrue
x = similar(b)
fill!(x, 0)
for i=1:10
x += M*x + b
end
x
6-element Vector{Float64}: -19.734466426485692 -23.660925449252744 22.7721733952521 14.346103166963537 -17.113723770023633 -34.32715946509404