## We are going to look at CG and MINRES in terms of their properties. 
using LinearAlgebra
using Krylov

A = SymTridiagonal(2*ones(20), -ones(19))
#A = SymTridiagonal(2*rand(20), -rand(19))

function lanczos(A,b,k)
  n = size(A,1)
  V = zeros(n,k+1)
  T = Tridiagonal(zeros(k-1), zeros(k), zeros(k-1))
  rho = 0.0
  V[:,1] = b/norm(b)
  for j=1:k
    y = A*V[:,j]
    for i=max(j-1,1):j
      T[i,j] = V[:,i]'*y
      y -= T[i,j]*V[:,i]
    end
    rho = norm(y)
    V[:,j+1] = y/rho
    if j < k
      T[j+1,j] = rho
    end
  end
  return V,T,rho
end

function myminres(A, b, k) 
  n = size(A,1) 
  V,T,rho = lanczos(A,b,k+1)
  normb = norm(b)
  residuals = zeros(n, k)
  solutions = zeros(n, k)
  for i=1:k
    xk = V[:,1:i]*((T[1:i+1,1:i]\((begin; z = zeros(i+1); z[1] = 1; z; end)*normb)))
    rk = b - A*xk
    residuals[:,i] = rk
    solutions[:,i] = xk
  end 
  return residuals, solutions 
end
R, S = myminres(A, ones(size(A,1)), 8)

function mycg(A, b, k) 
  n = size(A,1) 
  V,T,rho = lanczos(A,b,k+1)
  normb = norm(b)
  residuals = zeros(n, k)
  solutions = zeros(n, k)
  for i=1:k
    xk = V[:,1:i]*((T[1:i,1:i]\((begin; z = zeros(i); z[1] = 1; z; end)*normb)))
    rk = b - A*xk
    residuals[:,i] = rk
    solutions[:,i] = xk
  end 
  return residuals, solutions 
end

Rcg, Scg = mycg(A, ones(size(A,1)), 8)

([-9.000000000000002 3.552713678800501e-15 … 1.0658141036401503e-14 -7.105427357601002e-15; 1.0 -8.000000000000004 … 3.552713678800501e-15 -3.552713678800501e-15; … ; 1.0 -8.000000000000004 … 7.105427357601002e-15 -8.881784197001252e-15; -9.000000000000002 3.552713678800501e-15 … 1.0658141036401503e-14 -3.552713678800501e-15], [10.000000000000002 10.0 … 10.00000000000001 10.00000000000003; 10.000000000000002 19.000000000000004 … 19.000000000000032 19.000000000000053; … ; 10.000000000000002 19.000000000000004 … 19.000000000000032 19.000000000000057; 10.000000000000002 10.0 … 10.00000000000001 10.00000000000003])

## Compare against Krylov 
xcg, cgstats = Krylov.cg(A, ones(size(A,1)); itmax=5)
Scg[:,5] ≈ xcg 

true

## Compare against Krylov 
xmr, mrstats = Krylov.minres(A, ones(size(A,1)); itmax=5)
S[:,5] ≈ xmr 

true