% nsh.m ___________ Newton-Schultz-Hostelling
% @ K Chen (Nov 2002)
%__________________________________________________________________________ 
help nsh
if ~exist('A'), disp('New matrix A set for you');
  n =  7; A=rand(n); %%A=3*rand(n).*(rand(n)>.6); % sparse
          A=A+9*eye(n); 
end,       n=size(A,1), A0=A;
if ~exist('b'), disp('New RHS b set for you'); x=(1:n)'; b=A*x; clear x;end

 X=diag(1./diag(A));
for k=1:n                           %% Step k
 X=2*X-X*(A*X);
 x=X*b;
 r=A*x;
 r=b-r;
 r=norm(r,1)/norm(b,1);
 fprintf('nsh.m : step %2d  relative error = %e\n',k,r);
 if r<1.0e-6, iterations=k, break; end
end                                 %% k    ==>   Check solution x=b
  A=A0; clear A0
who A x X
 fprintf('\t nsh.m : check Ax=b and AX=eye(n)\n')
