using MatrixMarket, SuiteSparseMatrixCollection
using Krylov, LinearOperators
using LinearAlgebra, Printf

function residuals(A, b, shifts, x)
  nshifts = length(shifts)
  r = [ A' * (A * x[i] - b) + shifts[i] * x[i] for i = 1 : nshifts ]
  return r
end
ssmc = ssmc_db(verbose=false)
matrix = ssmc_matrices(ssmc, "HB", "well1033")
path = fetch_ssmc(matrix, format="MM")

A = MatrixMarket.mmread(joinpath(path[1], "$(matrix.name[1]).mtx"))
b = MatrixMarket.mmread(joinpath(path[1], "$(matrix.name[1])_b.mtx"))[:]
(m, n) = size(A)
@printf("System size: %d rows and %d columns\n", m, n)

# Define regularization parameters.
shifts = [1.0, 2.0, 3.0, 4.0]

(x, stats) = cgls_lanczos_shift(A, b, shifts)
show(stats)
r = residuals(A, b, shifts, x)

resids = map(norm, r) / norm(b)
@printf("CGLS: Relative residuals with shifts:\n")
for resid in resids
  @printf(" %8.1e", resid)
end
@printf("\n")

System size: 1033 rows and 320 columns
LanczosShiftStats
 niter: 15
 solved: true
 residuals: [Float64[], Float64[], Float64[], Float64[]]
 indefinite: Bool[0, 0, 0, 0]
 ‖A‖F: NaN
 κ₂(A): NaN
 timer: 40.10ms
 status: solution good enough given atol and rtol
CGLS: Relative residuals with shifts:
  1.3e-08  6.5e-09  1.0e-08  1.1e-08