clear all
ints=randi(20,3,1);     %   value of eigenvalues
sints=randi(10,3,1);     %   multiplicities
d=[ ints(1)*ones(sints(1),1);...
           ints(2)*ones(sints(2),1);...
                    ints(3)*ones(sints(3),1)];
d=diag(d);
xsol=randn(length(d),1);
[q,r]=qr(randn(length(d)));
A=q*d*q';
A=(A+A')/2;
b=A*xsol;
eig(A)
xi=randn(length(d),1);
for i=1:3,    % 3 distinct eigenvalues
	ri=b-A*xi;    % residual
	xi=xi+(1/ints(i))*ri;
end
residual=A*xi-b
