Friday, August 15, 2014

Gaussian Elimination


#include <vector>
#include <algorithm>
#include <iostream>
#include <cmath>
using namespace std;
const double eps = 1e-9;
enum sol {
 NOSOL, UNIQUE, INF
};
sol GaussianElimination(vector<vector<double> >& lhsMat, vector<double>& rhs,
  vector<double>&ans) {
 // the final answer will be in the ans
 int eqnCnt = (int) lhsMat.size(), varCnt = (int) lhsMat[0].size();
 sol ret = UNIQUE;
 for (int curCol = 0, curRow = 0; curCol < min(eqnCnt, varCnt);
   curCol++, curRow++) {
  for (int j = curRow; j < eqnCnt; j++) { // get a row where mat[i][i] is not zero
   if (fabs(lhsMat[j][curCol]) > eps) {
    swap(lhsMat[curRow], lhsMat[j]); // swapping vectors O(1)
    swap(rhs[curRow], rhs[j]);
    break;
   }
  }
  if (fabs(lhsMat[curRow][curCol]) < eps) { // if the selected cell is still 0, stay in the same row.
   curRow--;
   continue;
  }

  // make all cells below curCell 0
  for (int j = curRow + 1; j < eqnCnt; j++) {
   if (fabs(lhsMat[j][curCol]) < eps)
    continue;
   double mulByThis = -lhsMat[j][curCol] / lhsMat[curRow][curCol];
   for (int k = curCol; k < varCnt; k++) {
    lhsMat[j][k] += mulByThis * lhsMat[curRow][k];
   }
   rhs[j] += mulByThis * rhs[curRow];
  }
 }
 int calculatedValues = 0;
 ans = vector<double>(varCnt, 0.0);
 // go from bottom to top and get the answers
 for (int i = eqnCnt - 1; i >= 0; --i) {
  int k = 0;
  for (k = 0; k < varCnt; k++) {
   if (fabs(lhsMat[i][k]) > eps) {
    break;
   }
  }
  if (k == varCnt) {
   if (fabs(rhs[i]) > eps)
    return NOSOL;
   continue;
  }
  double val = rhs[i];
  for (int j = k + 1; j < varCnt; j++) {
   val -= ans[j] * lhsMat[i][j];
  }
  val /= lhsMat[i][k];

  ans[k] = val;
  calculatedValues++;
 }
 if (calculatedValues != varCnt)
  ret = INF;
 return ret;
}
int main(){
 int n, m;
 cin >> n >> m;
 vector<vector<double> > lhsMat(n, vector<double>(m));
 vector<double> rhs(n), ans;
 for(int i = 0; i < n; i++){
  for(int k = 0; k < m; k++){
   cin >> lhsMat[i][k];
  }
  cin >> rhs[i];
 }
 GaussianElimination(lhsMat,rhs,ans);
 for(int i = 0; i < (int) ans.size(); i++){
  cout << ans[i] << endl;
 }
 return 0;
}

No comments:

Post a Comment