Inf on preconditioner \ b

[ranjanan]

using LinearAlgebra: Tridiagonal
using SparseArrays: sparse
using AlgebraicMultigrid
n = 10_000
A = Tridiagonal(rand(n-1), rand(n), rand(n-1))
A = A + A'
ml = ruge_stuben(sparse(A))
M = aspreconditioner(ml)
b = rand(n)
M \ b

[ranjanan]

This is failing at the gauss seidel iteration, which is guaranteed to converge only for strongly diagonally dominant matrices. IterativeSolvers gauss seidel also seems to fail:

julia> IterativeSolvers.gauss_seidel(A, b)
10000-element Vector{Float64}:
   1.5800257460438275e8
  -1.3554143007232122e9
   3.678727322552991e9
  -2.9398100121906166e9
   9.609920630162582e8
  -9.265191029407964e8
   1.4179182468803146e9
  -1.953846198176463e9
   2.8147302728149977e9
  -2.2543662309779e9
   5.296700245710633e9
  -3.518278556801988e10
   7.519165127084735e11
  -9.174574863473485e11
   8.46692988019672e11
  -4.6588672065469116e11
   3.672492872009494e11
  -3.434314920660273e11
   2.2295164294441867e11
  -2.2530896628448767e11
   1.8095117462242758e11
  -1.0104953470842676e11
   6.047258133417195e11
  -2.651110781815847e12
   9.80595210606752e12
  -7.792056661898067e13
   2.329786496840364e13
  -1.1476908757642984e13
   4.62634925450458e12
  -5.074849546586282e12
   ⋮
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf
  Inf
 -Inf

[ranjanan]

using LinearAlgebra: Tridiagonal
using SparseArrays: sparse
using AlgebraicMultigrid
n = 10_000
A = Tridiagonal(rand(n-1), 100*rand(n), rand(n-1))
A = A + A'
ml = ruge_stuben(sparse(A))
M = aspreconditioner(ml)
b = rand(n)
M \ b

works.