/*
 * <日付を入力>
 * 例題6.6
 * 連立一次方程式の数値解法のクラス
 * Gaussの消去法，Gauss-Jordan法，逆行列計算を実装
 */


public class LinearEquationSolver3 {

  // フィールド
  int degree;  // 正方行列の次数
  double [][] coefficient_matrix;  // 係数行列
  double [] unknown_vector;  // 未知ベクトル
  double [] known_vector;  // 既知ベクトル

  // コンストラクタ
  LinearEquationSolver3(int n, double [][] matrix, double [] vector1, double [] vector2) {
    degree = n;  // 正方行列の次数を代入
    coefficient_matrix = matrix;  // 係数行列を代入
    unknown_vector = vector1;  // 解（未知）ベクトルを代入
    known_vector = vector2;  // 既知ベクトルを代入
  }

  // コンストラクタ
  LinearEquationSolver3(int n, double [][] matrix) {
    degree = n;  // 正方行列の次数を代入
    coefficient_matrix = matrix;  // 係数行列を代入
  }

  // Gaussの消去法を行うメソッド
  public double [] calculateGE() {  // GE: Gaussian Elimination

    double pivot;  // ピボット
    double solution_element;  // 未知ベクトルの各解要素
    double coefficient_element;  // 係数ベクトルの各係数要素

    // Gaussの消去法の計算
    // 前進消去
    for (int k=0; k<degree-1; ++k) {

      // k行目の計算
      // 第k行の式のピボットを代入．
      pivot = coefficient_matrix[k][k];

      // 係数行列k行目のj列以降を順にピボットで割る．
      for (int j=k+1; j<degree; ++j) {
        coefficient_matrix[k][j] /= pivot; 
      }

      // k行の既知ベクトルをピボットで割る．
      known_vector[k] /= pivot;


      // k+1行目以降の計算
      for (int i=k+1; i<degree; ++i) {
        // i行k列の係数要素をcoefficient_elementに代入
        coefficient_element = coefficient_matrix[i][k];

        for (int j=k+1; j<degree; ++j) {
          // 係数行列のi行j列からk行のcoefficient_matrix[i][k]倍を引く
          coefficient_matrix[i][j] -= coefficient_element * coefficient_matrix[k][j]; 
        }

        known_vector[i] -= coefficient_element * known_vector[k]; 
      } 
    }

    // 求める未知ベクトルの最終項が消去により求まる．
    unknown_vector[degree-1] = known_vector[degree-1] / coefficient_matrix[degree-1][degree-1];


    // 後退代入
    for (int k=degree-2; k>=0; --k) {
      // k行の定数ベクトルをsolution_elementに代入する．
      solution_element = known_vector[k];

      // 未知ベクトルのj行の値をk行の方程式に代入し，
      // 未知ベクトルのk行目を計算する．
      for (int j=k+1; j<degree; ++j) {
        solution_element -= coefficient_matrix[k][j] * unknown_vector[j]; 
      }

      // 計算結果を未知ベクトルのk行目に代入する．
      unknown_vector[k] = solution_element; 
    }
    

    return unknown_vector;  // 解行列を返す

  }  // calculateGE()


  // Gauss-Jordan法を行うメソッド
  public double [] calculateGJE() {  // GJE: Gaussian-Jordan Elimination

    double pivot;  // ピボット
    double solution_element;  // 未知ベクトルの各解要素
    double coefficient_element;  // 係数ベクトルの各係数要素

    // Gauss-Jordan法の計算
    for (int k=0; k<degree; ++k) {

      // k行目の計算
      // 第k行の式のピボットを代入．
      pivot = coefficient_matrix[k][k];

      // ピボットはピボット自身で割ると1になるので，計算を省く．
      coefficient_matrix[k][k] = 1.0;

      // 係数行列k行目のj列以降を順にピボットで割る．
      for (int j=k+1; j<degree; ++j) {
        coefficient_matrix[k][j] /= pivot; 
      }

      // k行の既知ベクトルをピボットで割る．
      known_vector[k] /= pivot;


      // k行目以外の計算
      // 最終的な定数ベクトルknown_vector[]が計算で上書きされて解となる．
      for (int i=0; i<degree; ++i) {
        if (i != k) {  // k行目以外の行のx項の消去．

          // i行k列の係数要素をcoefficient_elementに代入
          coefficient_element = coefficient_matrix[i][k];
          coefficient_matrix[i][k] = 0.0;

          for (int j=k+1; j<degree; ++j) {
            // 係数行列のi行j列からk行のcoefficient_matrix[i][k]倍を引く
            coefficient_matrix[i][j] -= coefficient_element * coefficient_matrix[k][j]; 
          }

          // 定数ベクトルのk行目以外の行からk行目のcoefficient_matrix[k+1][k]倍を引く．
          known_vector[i] -= coefficient_element * known_vector[k];
        }
      } 
    }

    return known_vector;  // 解ベクトルを返す

  }  // calculateGJE()


  /* 
   * ガウス・ジョルダン法による逆行列を求める関数 
   * および演習の関数を変更
   * 既知ベクトルknown_vector[]は使用しないので削除
   */
  public double [][] calculateInverseMatrix() {
    double pivot;  // ピボット
    double coefficient_element;  // 係数ベクトルの各係数要素

    // Gauss-Jordan法の計算
    for (int k=0; k<degree; ++k) {

      // k行目の計算
      // 第k行の式のピボットを代入．
      pivot = coefficient_matrix[k][k];

      // ピボットはピボット自身で割ると1になるので，計算を省く．
      coefficient_matrix[k][k] = 1.0;

      // 係数行列j行を順にピボットで割る．
      for (int j=0; j<degree; ++j) {
        coefficient_matrix[k][j] /= pivot; 
      }

      // k行目以外の計算
      // 元の行列に逆行列が上書きされる．
      for (int i=0; i<degree; ++i) {
        if (i != k) {  // k行目以外の行のx項の消去．

          // i行k列の係数要素をcoefficient_elementに代入
          coefficient_element = coefficient_matrix[i][k];
          coefficient_matrix[i][k] = 0.0;


          for (int j=0; j<degree; ++j) {
            // 係数行列のi行j列からk行のcoefficient_matrix[i][k]倍を引く
            coefficient_matrix[i][j] -= coefficient_element * coefficient_matrix[k][j]; 
          }
        }
      } 
    }

    return coefficient_matrix;  // 係数行列を返す

  }  // calculateInverseMatrix()


}  // class LinearEquationSolver