// gauss elimination
#include <fstream>
#include <iostream>
#include <iomanip>
#include <cmath>
enum { kSize = 3 };
typedef double SqMatT[][kSize];
void Dump(std::ostream& o, double const v[kSize])
{
for (int i = 0; i < kSize; ++i)
{
o << v[i] << '\t';
}
o << std::endl;
}
void Dump(std::ostream& o, SqMatT const m)
{
for (int i = 0; i < kSize; ++i)
{
Dump(o, m[i]);
}
}
void Eliminate(SqMatT a, double b[])
{
int n = kSize;
for (int i = 0; i < n -1; ++i)
{
for (int j = i + 1; j < n; ++j)
{
double c = a[j][i] / a[i][i];
for (int k = i; k < n; ++k)
{
a[j][k] -= c * a[i][k];
}
b[j] -= c * b[i];
}
}
}
void BackSubstitute(SqMatT U, double const c[], double x[])
{
for (int i = kSize - 1; 0 <= i; --i)
{
double s = 0.;
for (int j = i + 1; j < kSize; ++j)
{
s += U[i][j] * x[j];
}
x[i] = (c[i] - s)/U[i][i];
}
}
int main()
{
SqMatT m =
{
{ 2, 3, 0, },
{ 0, 1, 0, },
{ 0, 5, 1, },
};
double b[kSize] = { 2, 3, 4, };
Eliminate(m, b);
double x[kSize];
BackSubstitute(m, b, x);
//std::ofstream o("jjj.tsv");
Dump(std::cout, x);
}