// =============================================================================
// Gauss_pivot.sci: ガウスの消去法(ピボット選択を行う)で連立一次方程式を解く関数
// -----------------------------------------------------------------------------
//    入力 :  A : n行n列 の係数行列
//            b : n行1列 の定数ベクトル
//    出力:   x : n行1列 の連立一次方程式の解のベクトル
// =============================================================================
function x = Gauss_pivot(A, b)
[n, n1] = size(A);                 // 係数行列の縦横サイズ(n×n1)を調べる．
for i = 1:n-1                      // 係数行列の1〜n-1列目まで繰り返す．
   [m,k] = max(abs(A(i:n,i)));     // i列目のi行以降で係数の最大値を持つ行を探す．
                                   // k行目が最大値を持つ行 (実際は i+k-1行目)
   if k > 1                        // 行を入れ替える必要がある場合
      tmp1 = A(i,:);               // Aのi行目のベクトルをtmp1に一時的に退避する．
      tmp2 = b(i);                 // bのi行目の要素をtmp2に一時的に退避する．
      A(i,:) = A(i+k-1,:);         // Aのi行目を最大係数を持つi+k-1行目と入れ替える．
      b(i)   = b(i+k-1);           // bのi行目を最大係数を持つi+k-1行目と入れ替える．
      A(i+k-1,:) = tmp1;           // Aのi+k-1行目を退避してあったi行目に置き換える．
      b(i+k-1)   = tmp2;           // bのi+k-1行目を退避してあったi行目に置き換える．
   end

   for h = i+1:n                    // i+1行目以降に対して掃き出し法を適用する．
      r = A(h,i)/A(i,i);            // h行目係数のが対角要素の何倍かを求める．
      A(h,:) = A(h,:) - r * A(i,:); // 掃き出し法でAのh行目を更新する．(i列目を0にする)
      b(h)   = b(h)   - r * b(i);   // 掃き出し法でbのh行目を更新する．
   end;
end;                                // この時点で，掃き出し法の上三角行列ができ上がる．

x = zeros(n,1);                     // 解のベクトル x の要素を全て0で初期化する．
x(n) = b(n) / A(n,n);               // まず一番下の行(n行目)からn番目の解x(n)を計算する．

for i = n-1:-1:1                    // 上三角行列の下の方から順番に解を求めていく．
   x(i) = (b(i) - A(i,i+1:n) * x(i+1:n,1)) / A(i,i);
end                                 // ↑上三角行列i行目でi番目の解x(i)を計算する．
endfunction