%-------------------------------------------------------------------------%
%
% File: Lanczos(IM,rhs,tol,mmax,rbf,ep)
%
% Goal: perform the Lanczos procedure
%
% Inputs: IM:        the local interpolation matrix
%         rhs:       the function values
%         tol:       tolerance used during the Lanczos procedure
%         mmax:      maximum number of iterations for the Lanczos
%                    procedure
%         rbf:       the radial basis funcion
%         ep:        shape parameter
%
% Outputs: H:  the matrix change of basis
%          P:  the Lanczos basis
%
%-------------------------------------------------------------------------%
function [H,P] = Lanczos(IM,rhs,tol,mmax,rbf,ep)
% Initialize
n = size(IM, 1);
P = zeros(n, mmax + 2);
alpha = zeros(mmax, 1);
beta = zeros(mmax + 1, 1);
P(:, 2) = rhs / norm(rhs);
for j = 1 : mmax
    w = IM * P(:, j+1)- beta(j) * P(:, j); 
    alpha(j) = w' * P(:, j+1);
    w = w - alpha(j) * P(:, j+1);
    for i = j + 1 : -1 : 1
        w = w - (w' * P(:, i)) * P(:, i);
    end
    beta(j + 1) = norm(w);
    if beta(j+1)== 0
        P = P(:, 1 : j + 2);
        alpha = alpha(1 : j);
        beta = beta(1 : j + 1);
        break
    end    
    P(:, j + 2) = w / beta(j+1);
    if   abs(rbf(ep,0) - sum(alpha) / n) < tol 
        P = P(:, 1 : j + 2);
        alpha = alpha(1 : j);
        beta = beta(1 : j + 1);
        break
    end
end
if j == mmax
    w = IM * P(:, mmax + 1);
    alpha(mmax) = w' * P(:, mmax + 1);
end
P = P(:, 2 : end);
beta = beta(2 : end);
H = diag(alpha) + diag(beta(1 : end-1), 1) + diag(beta(1 : end-1), -1);
H(end +1, end) = beta(end);
end