% Mfun = pc_gs(A)
%
% Generates a symmetric Gauss-Seidel preconditioner.
%
% Input:
%   A - Original matrix 
%
% Output:
%   Mfun - function to solve with the preconditioner M
%
% Usage:
%   Mfun = pc_2grid(A);
%   x = pcg(A,b,tol,maxit,Mfun);

function Mfun = pc_gs(A)

  L = tril(A);
  U = triu(A);
  Mfun = @(r) apply_pc(L,U,A,r);
  
function x = apply_pc(L,U,A,r)

  x = L\r;
  x = x + U\(r-A*x);

