function [x,p] = dc(A,C,b,n,z)
%
%   [x,res] = dc(A,C,b,n,z)
%    
%      Defect correction iteration on system A*x = b with printing
%      n -- number of iterations
%      z -- initial vector  (default 0)
%      b -- right hand side of system
%
%      x -- final iterate
%      A -- system matrix
%      C -- Approximate inverse of A according to the scheme
%           x^{k+1} = x^k - C*(A*x^k - b)
%

if nargin <=4, z=0*b; end
p=zeros(n,2);
m = size(A,1); 
x=z;
for k = 1:n
   x = x - C*(A*x - b);
   defect = A*x -b;
   res = sqrt(defect'*defect);
   fprintf(1,'%3d     ',k)
   fprintf(1,'%5.4f     ',x')
   fprintf(1,'res=%5.4f     ',res)
   fprintf(1,'\n')
   p(k,1)= k;
   p(k,2)= res;
end
plot(p(:,1),p(:,2));
%semilogy(p(:,1),p(:,2));