#include <iostream>
#include <vector>
#include <cstdlib>
#include <iomanip>

const int N1 = 3;
const int N2 = 4;

template <class T>
class LUDecompose {
  typedef std::vector<std::vector<T> > VV;
  typedef std::vector<T> V;
  VV a, l, u, x, y;
  V b;
  int n;

  void MakeSquareMatrix(VV& a) {
    a.resize(n);
    for (int i = 0; i < n; i++)
      a[i].resize(n);
  }
  void MakeVector(V& b) {
    b.resize(n);
  }
  // LU分解(ドゥーリトル法)
  void LUDecomp(VV& a, VV& l, VV& u) {
    int i, j, k;

    for (i = 0; i < n; i++)
      l[i][i] = 1.0;
    for (i = 0; i < n; i++)
      u[0][i] = a[0][i];

    for (i = 1; i < n; i++) {
      if (u[i - 1][i - 1] == 0.0) {
        std::cerr << "行列Uは正則ではありません" << std::endl;
        std::exit(1);
      }
      for (j = i; j < n; j++) {
        // L成分の算出
        l[j][i - 1] = a[j][i - 1];
        for (k = 0; k < i - 1; k++)
          l[j][i - 1] -= l[j][k] * u[k][i - 1];
        l[j][i - 1] /= u[i - 1][i - 1];
        // U成分の算出
        u[i][j] = a[i][j];
        for (k = 0; k < i; k++)
          u[i][j] -= l[i][k] * u[k][j];
      }
    }
  }
  void PrintMatrix(VV& a, const char* s) {
    std::cout << "*** " << s << "***\n";

    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) {
        std::cout << std::setw(11) << std::setprecision(7) << std::fixed
          << std::showpoint << std::right << a[i][j];
        std::cout << "  ";
      }
      std::cout << '\n';
    }
  }
  // LU分解を利用して逆行列を求める
  void LUInv(VV& l, VV& u, VV& x) {
    int i, j, k;
    double r;
    VV lu;

    // 上三角行列と下三角行列を合成する
    MakeSquareMatrix(lu);

    for (i = 0; i < n; i++)
      for (j = 0; j < n; j++)
        lu[i][j] = (i <= j) ? u[i][j] : l[i][j];

    // 逆行列を求める
    for (k = 0; k < n; k++) {
      for (i = 0; i < n; i++) {
        r = (i == k) ? 1.0 : 0.0;
        for (j = 0; j < n; j++)
          r -= lu[i][j] * x[j][k];
        x[i][k] = r;
      }
      for (i = n - 1; i >= 0; i--) {
        r = x[i][k];
        for (j = i + 1; j < n; j++)
          r -= lu[i][j] * x[j][k];
        x[i][k] = r / lu[i][i];
      }
    }
  }
  // 行列の積
  void MatMul(VV& z, VV& x, VV& y) {
    int i, j, k;
    double s;

    for (i = 0; i < n; i++)
      for (j = 0; j < n; j++) {
        s = 0.0;
      for (k = 0; k < n; k++)
        s += x[i][k] * y[k][j];
      z[i][j] = s;
    }
  }
  /* 三角行列の行列式は対角成分の積 */
  double DetTriangularMatrix(VV& a) {
    double det = 1.0;
    
    for (int i = 0; i < n; i++)
      det *= a[i][i];
  
    return det;
  }

public:
  LUDecompose(int nn, T* init) : n(nn) {
    MakeSquareMatrix(a);
    for (int i = 0; i < n; i++)
      for (int j = 0; j < n; j++)
        a[i][j] = init[i * n + j];
    MakeSquareMatrix(l);
    MakeSquareMatrix(u);
    MakeSquareMatrix(x);
    MakeSquareMatrix(y);
  }
  void Calculate() {
    LUDecomp(a, l, u);
    LUInv(l, u, x);
    MatMul(y, a, x);
  }
  void PrintMatrix() {
    PrintMatrix(a, "行列 A");
    PrintMatrix(l, "行列 L");
    PrintMatrix(u, "行列 U");
    PrintMatrix(x, "逆行列 X(A^(-1))");
    PrintMatrix(y, "行列 Y(A * X)");
    std::cout << "行列式 |A| = " << DetTriangularMatrix(l) * DetTriangularMatrix(u) << std::endl;
  }
};

int main(void)
{
  double a1[] = {2.0, 3.0, -1.0, 4.0, 4.0, -3.0, 2.0, -3.0, 1.0};
  double a2[] = {1.0, 1.0, 0.0, 3.0, 2.0, 1.0, -1.0, 1.0, 3.0, -1.0, -1.0, 2.0, -1.0, 2.0, 3.0, -2.0};

  LUDecompose<double> lu1(N1, a1);
  lu1.Calculate();
  lu1.PrintMatrix();

  std::cout << '\n';

  LUDecompose<double> lu2(N2, a2);
  lu2.Calculate();
  lu2.PrintMatrix();
}