// =============================================================================
// Seidel.sci:  ガウス・ザイデル法で連立一次方程式 Ax=b を解く関数
// -----------------------------------------------------------------------------
// 入力  A:   係数行列A (n × n)
//       b:   右辺bの縦ベクトル (n × 1)
//       x0:  解の初期値の縦ベクトル (n × 1)
//       tol: ノルムの変化の2乗の値の許容範囲 (収束判定に利用)
//       itr: 最大繰返し回数
// 出力  x:   解xの縦ベクトル (n × 1)
// =============================================================================
// C: 式(5.41)の係数a_ij/a_ii の値を入れる行列 (右辺第1,2項目の係数÷対角成分)
// r: 式(5.41)の係数b_i/a_iiの値を入れる縦ベクトル (右辺第3項目の係数÷対角成分)
// S: 解の途中結果を保存するための行列

function x = Seidel(A, b, x0, tol, itr)
[n, m] = size(A);               // 行列の縦横サイズをチェック
C = -A;                         // C = -A (行列Aを右辺へ移項)
for i = 1:n       
  C(i,i) = 0;                   // 行列Cの対角成分を全て0に
  C(i,1:n) = C(i,1:n) / A(i,i); // 行列Cの各行の要素をAの対角成分で割る
  r(i) = b(i) / A(i,i);         // ベクトルbの要素をAの対角成分で割ったもの
end

S = x0;                         // 解の初期値を行列Sの一番左の列に保存
x = x0;                         // 解の初期値をxに代入
for k = 1:itr                   // 繰返し回数の最大値itr回まで繰り返す
  for i = 1:n                   // i=1,2,…の順にi番目の解x(i)を求める
    x(i) = C(i,:) * x + r(i);   // x(i)の結果が次の解x(i+1)に反映される
  end
  S(:,k+1) = x;                 // k回目繰返しにおける解を行列Sの右側に追加
  if norm(S(:,k) - x)^2 < tol   // 前の解との差のノルムの2乗を見て収束判定
    break;
  end
end
format(7);                      // 解の途中経過を5桁までの表示に変更
disp(S');                       // 途中経過を表示(Sを転置して表示)

endfunction