# C++ my Gauss–Seidel method doesn’t work wery well

I want to solve system of equations (3×3) in C++ useing Gauss–Seidel method. I read some of pages about this: Wiki Page2 and Page3. By pattern oneself on those pages I wrote this code:

```#include <iostream>
#include <cmath>
#define epsilon 0.001

int main()
{
//declaration and initialization of variables
long double a[3][3], b[3], x[3], xo[3];

a[0][0] = 3;   a[0][1] = -6.5; a[0][2] = 3.33; b[0] = 7;
a[1][0] = 4.7; a[1][1] = 8;    a[1][2] = -17;  b[1] = 5;
a[2][0] = 6;   a[2][1] = 2;    a[2][2] = 8;    b[2] = -5.2;
//By calculator: x=0,5063353908 y=-1,215147148 z=-0,725964756

int n = 3;

//initialization of x array
for (int i = 0; i < n; i++)
{
x[i] = b[i] / a[i][i];
}

//start Gauss-Seidel method
do
{
for (int i = 0; i < n; i++)
{
//xo is useing to chceck accuracy
xo[i] = x[i];
for (int j = 0; j < n; j++)
{
if (i != j) x[i] -= (a[i][j] * x[j]) / a[i][i];
}
}
//chceck the accuracy
} while (fabs(x[0] - xo[0]) < epsilon &&
fabs(x[1] - xo[1]) < epsilon &&
fabs(x[2] - xo[2]) < epsilon);

std::cout << "x = " << x[0] << std::endl;
std::cout << "y = " << x[1] << std::endl;
std::cout << "z = " << x[2] << std::endl;

return 0;
}
```

This code doesn’t work very well – the result of the program isn’t the same like solving by calculator. I really don’t have an idea what is going on.

There are several issues in your code:

1. The accuracy checking is not correct
2. In calculation of the solution, a term `b[i]` is missing
3. The matrix is not diagonal dominant, far from it: the algorithm cannot converge

In addition, to avoid getting `nan` when the algorithm doesn’t converge, I have added a condition to limit the number of iterations.

I have also added a final check, i.e. calculation of `ax`, to see if a good result is effectively obtained.

```#include <iostream>
#include <cmath>
#define epsilon 0.001

int main()
{
//declaration and initialization of variables
long double a[3][3], b[3], x[3], xo[3];

a[0][0] = 3;   a[0][1] = -6.5; a[0][2] = 3.33; b[0] = 7;
a[1][0] = 4.7; a[1][1] = 8;    a[1][2] = -17;  b[1] = 5;
a[2][0] = 6;   a[2][1] = 2;    a[2][2] = 8;    b[2] = -5.2;
//By calculator: x=0,5063353908 y=-1,215147148 z=-0,725964756

a[0][0] = 12;       // to make matrix diadonal dominant
a[1][2] = -1;

int n = 3;

//initialization of x array
for (int i = 0; i < n; i++) {
x[i] = b[i] / a[i][i];
}

//start Gauss-Seidel method
int iter = 0;
int n_iter_max = 20;
do {
for (int i = 0; i < n; i++) {
//xo is used to check accuracy
xo[i] = x[i];
x[i] = b[i];
for (int j = 0; j < n; j++) {
if (i != j) x[i] -= (a[i][j] * x[j]);
}
x[i] /= a[i][i];
}
iter++;
//check the accuracy
} while ((std::abs(x[0] - xo[0]) > epsilon ||
std::abs(x[1] - xo[1]) > epsilon ||
std::abs(x[2] - xo[2]) > epsilon) && (iter <= n_iter_max));

std::cout << "x = " << x[0] << std::endl;
std::cout << "y = " << x[1] << std::endl;
std::cout << "z = " << x[2] << std::endl;

//  Final check

double b_cal[3];
for (int i = 0; i < n; ++i) {
b_cal[i] = 0.0;
for (int j = 0; j < n; j++) {
b_cal[i] += a[i][j] * x[j];
}
}
std::cout << "A*x = ";
for (int i = 0; i < n; ++i) {
std::cout << b_cal[i] << " (" << b_cal[i] - b[i] << ")  ";
}
std::cout << "n";

return 0;
}
```