using Krylov
using LinearAlgebra, Printf

m = 5
n = 8
λ = 1.0e-3
A = rand(m, n)
b = A * ones(n)
xy_exact = [A  λ*I] \ b # In Julia, this is the min-norm solution!

(x, y, stats) = craig(A, b, λ=λ, atol=0.0, rtol=1.0e-20, verbose=1)
show(stats)

# Check that we have a minimum-norm solution.
# When λ > 0 we solve min ‖(x,s)‖  s.t. Ax + λs = b, and we get s = λy.
@printf("Primal feasibility: %7.1e\n", norm(b - A * x - λ^2 * y) / norm(b))
@printf("Dual   feasibility: %7.1e\n", norm(x - A' * y) / norm(x))
@printf("Error in x: %7.1e\n", norm(x - xy_exact[1:n]) / norm(xy_exact[1:n]))
if λ > 0.0
  @printf("Error in y: %7.1e\n", norm(λ * y - xy_exact[n+1:n+m]) / norm(xy_exact[n+1:n+m]))
end

CRAIG: system of 5 equations in 8 variables
    k       ‖r‖       ‖x‖       ‖A‖      κ(A)         α        β  timer
    0  8.21e+00  0.00e+00  0.00e+00  0.00e+00   ✗ ✗ ✗ ✗  ✗ ✗ ✗ ✗  0.00s
    1  1.23e-01  2.78e+00  2.95e+00  2.95e+00   2.9e+00  4.4e-02  0.00s
    2  1.03e-01  2.79e+00  3.11e+00  4.40e+00   7.6e-01  6.3e-01  0.00s
    3  1.18e-01  2.80e+00  3.20e+00  5.82e+00   5.0e-01  5.7e-01  0.00s
    4  1.68e-03  2.80e+00  3.32e+00  7.50e+00   8.8e-01  1.3e-02  0.00s
    5  1.13e-11  2.80e+00  3.41e+00  8.43e+00   7.7e-01  5.2e-09  0.00s

SimpleStats
 niter: 5
 solved: true
 inconsistent: false
 indefinite: false
 npcCount: 0
 residuals: []
 Aresiduals: []
 qvals: []
 κ₂(A): []
 allocation timer: 2.66μs
 timer: 90.27μs
 status: solution good enough for the tolerances given
Primal feasibility: 1.4e-12
Dual   feasibility: 2.3e-16
Error in x: 1.4e-12
Error in y: 9.7e-13