#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <cstdlib>
#include <cassert>
#include <cmath>

enum ErrorCode { None = 0, Singular, BadFormat, NotFound };
class Exception {
  enum ErrorCode ErrorCode;
public:
  Exception() : ErrorCode(None) {}
  Exception(enum ErrorCode e) : ErrorCode(e) {}
  enum ErrorCode getErrorCode() { return ErrorCode; }
};

class Matrix {
public:
  int n;
  std::vector<std::vector<double> > m;
  Matrix() : n(0) { m.clear(); };
  Matrix(int n, double diag);
  void readFromFile(char const *filename);

  friend std::ostream &operator<<(std::ostream &stream, Matrix m);
  friend Matrix operator*(Matrix &a, Matrix &b);
};

Matrix::Matrix(int n, double diag) {
  this->n = n;
  std::vector<double> v;
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++)
      if (i == j)
        v.push_back(diag);
      else
        v.push_back(0.0);
    m.push_back(v);
    v.clear();
  }
}

void Matrix::readFromFile(char const *filename) {
  std::ifstream fin(filename);
  if(fin) {
    std::vector<double> line;
    std::string buff;
    int index = 0;
    while (getline(fin, buff)) {
      double data;
      std::stringstream ss;
      ss << buff;
      int i;
      for (i = 0; ss >> data; i++) {
        line.push_back(data);
        ss.ignore();
      }
      m.push_back(line);
      line.clear();
      if (index == 0) {
        this->n = i;
      } else if (i != this->n) {
          std::cerr << "line: " << index << "'s colomn num is not " << this->n << " but " << i << std::endl;
          throw Exception();
      } else if (index >= this->n) {
          std::cerr << "line: " << index << " is too large: " << this->n << std::endl;
          throw Exception();
      }
      index++;
    }
    fin.close();
  } else {
    throw Exception(NotFound);
  }
}

std::ostream &operator<<(std::ostream &stream, Matrix m) {
  std::vector<std::vector<double> >::iterator p;
  std::vector<double>::iterator q;
  for (p = m.m.begin(); p != m.m.end(); p++) {
    for (q = (*p).begin(); q != (*p).end(); q++)
      stream << *q << ", ";
    stream << std::endl;
  }
  return stream;
}

Matrix operator*(Matrix &a, Matrix &b) {
  Matrix c(a.n, 0.0);
  for (int i = 0; i < a.n; i++)
    for (int j = 0; j < b.n; j++) {
      double s = 0.0;
      for (int k = 0; k < c.n; k++)
        s += a.m[i][k] * b.m[k][j];
      c.m[i][j] = s;
    }
  return c;
}

class GaussianEliminationPartial : public Matrix {
public:
  void operator()(Matrix &m, Matrix &e) {
    assert(m.n == e.n);
    /* forward */
    for (int i = 0; i < m.n; i++) {
      /* partial-pivotting */
      double max = std::fabs(m.m[i][i]);
      int k = i;
      for (int j = i + 1; j < m.n; j++) {
        double tmp;
        if (max < (tmp = std::fabs(m.m[j][i]))) {
          max = tmp;
          k = j;
        }
      }
      if (k != i)
        for (int j = 0; j < m.n; j++) {
          double r;
          r = m.m[k][j];
          m.m[k][j] = m.m[i][j];
          m.m[i][j] = r;
          r = e.m[k][j];
          e.m[k][j] = e.m[i][j];
          e.m[i][j] = r;
        }
      /* pivotting finished */
      if (std::fabs(max) < 1.0e-6) {
        throw Exception(Singular);
      }

      double r;
      r = m.m[i][i];
      for (int j = 0; j < m.n; j++) {
        m.m[i][j] /= r;
        e.m[i][j] /= r;
      }
      for (int j = i + 1; j < m.n; j++) {
        if ((r = m.m[j][i]) != 0) {
          assert(std::fabs(r) > 0);
          for (int k = 0; k < m.n; k++) {
            m.m[j][k] /= r;
            m.m[j][k] -= m.m[i][k];
            e.m[j][k] /= r;
            e.m[j][k] -= e.m[i][k];
          }
        }
      }
    }
    /* backward */
    for (int i = m.n - 1; i >= 0; --i) {
      for (int j = i - 1; j >= 0; --j) {
        e.m[j][0] -= e.m[i][0] * m.m[j][i];
        for (int k = 1; k < e.n; k++) {
          e.m[j][k] -= e.m[i][k] * m.m[j][i];
        }
        m.m[j][i] = 0.0;
      }
    }
  }
};

Matrix inverse(Matrix m) {
  Matrix e(m.n, 1.0);
  GaussianEliminationPartial ge;
  ge(/*ref*/m, /*ref*/e);
  std::cout << m << '\n' << e << std::endl;
  return e;
}

char const *fileName = "matrix.csv";
int main() {
  Matrix a, b;
  try {
    a.readFromFile(fileName);
    std::cout << "a = \n" << a << std::endl;
    b = inverse(a);
    std::cout << "b = \n" << b << std::endl;
    std::cout << "a * b = \n" << a * b << std::endl;
  } catch (Exception e) {
    switch (e.getErrorCode()) {
    case NotFound:
      std::cerr << "cannnot open the file." << std::endl;
      break;
    case BadFormat:
      std::cerr << "matrix: bad format." << std::endl;
      break;
    case Singular:
      std::cout << "not have an inverse matrix." << std::endl;    
      break;
    default:
      std::cout << "Internal Error." << std::endl;          
      ;
    }
  }
  return 0;
}
/* end */