function [Q,T]=lanczos(A)
% LANCZOS: Outputs the factorization A = Q T Q^T, where Q is orthogonal and T is tridiagonal
%
% INPUTS:
%   A : nxn matrix (must be symmetric)
%
% OUTPUT:
%   Q : orthogonal matrix 
%   T : tridiagonal matrix
%       such that A = Q T Q^T
% 
% Written by:
% -- 
% John L. Weatherwax                2006-12-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

% get the dimension of A:
n = size(A,1); 

% initialized some space: 
a = zeros(n,1); 
b = zeros(n,1); 
R = zeros(n,n); 
Q = zeros(n,n); 

% this is the "i=0" step of the loop (indices are shifted by ONE left):
i=1; 
b_0 = 1;
rand('seed',0); R_0 = rand(n,1); R_0 = R_0/norm(R_0,2);
Q_0 = zeros(n,1); 
Q(:,i) = R_0/b_0;
  
for i=1:n-1,
  a(i) = (Q(:,i).')*A*Q(:,i);
  if( i>1 ) 
    R(:,i) = A*Q(:,i) - b(i-1)*Q(:,i-1) - a(i)*Q(:,i); 
  else
    R(:,i) = A*Q(:,i)                   - a(i)*Q(:,i); 
  end
  b(i) = norm(R(:,i),2); 
  
  Q(:,i+1) = R(:,i)/b(i);
end
i=n; a(i) = (Q(:,i).')*A*Q(:,i);

% assemble the matrix T: 
T = diag(b(1:end-1),-1) + diag(a) + diag(b(1:end-1),+1); 

% lets check our results (if desired): 
if( 0 ) 
  disp( 'unit testing...' ); 
  disp('Q^T Q is equal to' ); 
  (Q.') * Q 
  disp( 'A*Q is equal to' ); 
  A * Q 
  disp( 'Q*T is equal to' ); 
  Q * T 
end