#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;
}
Friday, August 15, 2014
Gaussian Elimination
Subscribe to:
Post Comments (Atom)

No comments:
Post a Comment