Решение СЛАУ с использованием метода релаксации

#include <iostream>
 
#define N 3
 
int main()
{
  const double omega = 0.5; // 0.0 < omega < 2.0
  const double eps = 0.000001;
  const int MAX_STEPS = 1000;
  double xk[N] = {0.0, 0.0, 0.0}; // начальная точка
 
  // Исходное СЛАУ
  double A[N][N] = {{3.0, -1.0, 1.0}, {2.0, -1.0, 3.0}, {1.0, -1.0, -0.5}};
  double b[N] = {10.0, 3.5, 0.75};
 
  // Новое СЛАУ
  double B[N][N];
  double c[N];
 
  double xk1[N]; // x^{(k+1)}
 
  int step = 0;
  double norm2;
 
  // Вычисление новой СЛАУ
  for (int i = 0; i < N; ++i) {
    c[i] = b[i] / A[i][i];
    for (int j = 0; j < N; ++j) {
      B[i][j] = -A[i][j] / A[i][i];
    }
  }
 
  do {
    // Следующая точка
    for (int i = 0; i < N; ++i) {
      xk1[i] = (1.0 - omega) * xk[i];
      double s = c[i];
      for (int j = 0; j < i; ++j) {
        s += B[i][j] * xk1[j];
      }
      for (int j = i + 1; j < N; ++j) {
        s += B[i][j] * xk[j];
      }
      xk1[i] += omega * s;
    }
 
    // вычисление нормы разности в квадрате
    norm2 = 0.0;
    for (int i = 0; i < N; ++i) {
      norm2 += (xk[i] - xk1[i]) * (xk[i] - xk1[i]);
    }
 
    // Сохраняем x^{k+1} в x^{k}
    for (int i = 0; i < N; ++i) {
      xk[i] = xk1[i];
    }
 
    // Печатаем информацию о текущем шаге
    ++step;
    std::cout << "step: " << step << std::endl;
    std::cout << "  norm^2: " << norm2 << std::endl;
    std::cout << "  x:";
    for (int i = 0; i < N; ++i) {
      std::cout << " " << xk[i];
    }
    std::cout << std::endl;
    std::cout << "  A*x: ";
    double r[N];
    for (int i = 0; i < N; ++i) {
      r[i] = 0.0;
      for (int j = 0; j < N; ++j) {
        r[i] += A[i][j] * xk[j];
      }
      std::cout << " " << r[i];
    }
    std::cout << std::endl;
  } while ((norm2 >= eps * eps) && (step < MAX_STEPS));
 
  return 0;
}