function [x,G,c] = gsmp(A,b,n,z)
%
%   x = gsmp(A,b,n,z)
%    
%      Gauss-Seidel iteration on system A*x = b with printing
%         using matrix multiplication  (not optimal)
%      n -- number of iterations
%      z -- initial vector  (default 0)
%
%      x -- final iterate
%      G -- Gauss-Seidel matrix
%      c -- Gauss-Seidel vector
%

if nargin <=3, z=0*b; end

LD = tril(A);
G = -LD\triu(A,1);
c = LD\b;
x=z;
for i = 1:n
   x = G*x + c;
   fprintf(1,'%3d     ',i)
   fprintf(1,'%5.5f     ',x')
   fprintf(1,'\n')
end