function [H,V] = reorth_arnoldi(A,u,degree)
%Input:
%   A:      Input matrix (best of sparse)
%   u:      Starting vector
%   degree: Number of steps
%
%Output:
%   H:      (degree+1) x degree upper Hessenberg matrix
%   V:      n x (degree+1) matrix that holds the Krylov basis V 
%
% Performs modified Gramm Schmidt Arnoldi on A for "degree" steps
% and then increases H by one line. A reorthogonalization scheme
% is applied to assure the orthogonality between v vectors.

[m,n] = size(A);
threshold = sqrt(eps)*normest(A);  %by Paige's theorem

H=zeros(degree+1,degree);
V=zeros(m,degree+1);
u = u/norm(u);
V(1:m,1) = u;
for j=1:degree,
   w = A * V(1:m,j);
   for i = 1:j,
      H(i,j) = w.'* V(1:m,i);
      w = w - H(i,j) * V(1:m,i);
   end;
   H(j+1,j) = norm(w);
   H(j+2:degree+1,j)=0;
   if ( H(j+1,j) <= eps)
      break;
      H(j+1,j)
   end;   
   V(:,j+1) = w/H(j+1,j);
   w=V(:,j+1);
   
   %check the orthogonality level between vectors
   for k=1:j
      dot=w'*V(:,k);
      if dot > (threshold)  
         w=w-dot*V(:,k);  %reorthogonalize w against V(:,k)         
      end   
   end   
   V(:,j+1)=w/norm(w);
   
end;