function DD=qralg(N,fflag)
%
% naive implementation of the QR algorithm -- for 
% demo/illustration purposes only.
% If N is a scalar, use it for A=rand(N,N)
% If A is a matrix, use it directly (symmetrize if needed)
%    (typical use: sort some numbers!!!!!!!)
% fflag = 1 is the default, always symmetrize A=A'*A
% fflag = 0 allows experiments w/ nonsymmetric A
%
% qralg(diag(floor(20*rand(10,1)-10)),0);
%       nothing happens
% rand perturb
% S=rand(10,10);qralg(S*diag(floor(20*rand(10,1)-10))/S,0);
% All rights reserved. June 2006. Matthias Kawski
% http://math.asu.edu/~kawski

if nargin < 1
    %size
    N=8;
    A=rand(N,N);
else
   if max(size(N)) > 1
       A=N;
       N=size(A,1);
   else
       A=rand(N,N);
   end
end

if nargin < 2
   fflag = 1;
end;
       
% skip this for experiments with non-symmetric matrices
if 1-fflag
   if max(max(abs(A-A')))>0, A=A'*A; end
end

A
D=eig(A);
% number of iterations
n=10;
DD=zeros(N,n+1);
DD(:,1)=diag(A);
for k=[1:n]
    [Q,R]=qr(A);
    A=R*Q
    DD(:,k+1)=diag(A);
end
clf
plot([1;n+1],[1;1]*eig(A)','c-')
hold on
plot(DD')
pause
axis([0 n+1 -3 3])