function x = krylov(A,b,m,S);
%KRYLOV 


%         Reference: Y. Saad, Mathematics of Computation ...
 
%         Ellen McGrattan, 9-5-91
%         revised, ERM, 9-20-92 

N = length(b);
if nargin < 3; 
  m = round(.25*N); 
  S = 3;
end;
x = zeros(b);
r = b;
for s=1:S;
 v(:,1) = r/norm(r);
 for k = 1:m;
   sum1 = 0;
   for i = 1:k;
     h(i,k) = (A*v(:,k))'*v(:,i);
     sum1   = sum1+h(i,k)*v(:,i);
   end;
   w = A*v(:,k)-sum1;
   h(k+1,k) = norm(w);
   v(:,k+1) = w/h(k+1,k);
 end;
 y=norm(r)*(h(1:m,:)\(v(:,1:m)'*v(:,1)));
 z=v(:,1:m)*y;
 x=x+z;
 r=r-A*z;
end;
