function T = trireducedemo(A)
%trireduce Orthogonal transformation of symmetric matrix A to a symmetric
%tridiagonal matrix.
%
%   T = trireduce(A) - assigns to T a symmetric tridiagonal matrix
%   orthogonally similar to A.
%
%   [P T] = trireduce(A) - assigns to T a symmetric tridiagonal
%   matrix orthogonally similar to A and an orthogonal matrix
%   P such that P'*A*P = T.

[m, n] = size(A);
if m ~= n
	error('The matrix must be square.');
end

if issym(A) == false
   error('The matrix is not symmetric.');
end

spy(A);
pause;

for i = 1:n-2
	[u, beta, sign] = houseparms(A(i+1:n,i));
	% efficiently compute Hi*A_nmi*Hi
	p = beta*A(i+1:n,i+1:n)*u;
	K = (beta*u'*p)/2;
	q = p - K*u;
	% both A(i+i,i) and A(i, i+1) have value plus or minus
	% the norm of [A(i+1,i) A(i+1,i) ... A(n,i)]'
	A(i+1,i) = sign*norm(A(i+1:n,i));
	A(i,i+1) = A(i+1,i);
	% now compute the orthogonal transformation
	A(i+1:n,i+1:n) = A(i+1:n,i+1:n) - u*q' - q*u';
   for j = i+2:n
      A(i,j) = 0;
      A(j,i) = 0;
   end
   A(abs(A)<1.e-7)=0;
   spy(A);
   pause;
end

T = diag(diag(A,-1),-1) + diag(diag(A)) + diag(diag(A,1),1);
