function s=Gauss(a,b,x0,eps)
  % Gauss-Seidel迭代法解線性方程組
  % a為系數(shù)矩陣，b為方程組ax=b的右端項(xiàng)，x0為初值
  if nargin==3
      eps=1.0e-6;
  elseif nargin<3
      error
      return
  end
  D=diag(diag(a)); %求對(duì)角矩陣
  L=tril(a,-1); %求嚴(yán)格下三角矩陣
  U=triu(a,1); %求嚴(yán)格上三角矩陣
  C=inv(D+L);
  B=-C*U;
  f=C*b;
  s=B*x0+f;
  while norm(s-x0)>=eps
      x0=s;
      s=B*x0+f;
  end
  return