QR ver 2 ЧМ 1 лаба

#include <iostream>

using namespace std;

const double PI = 3.1415926536;

bool isSimmetrial(double ** coefficients, int numberOfEquation) {
  bool result = true;
  int i, j;
  for (i = 0; i < numberOfEquation; i++) {
    for (j = i + 1; j < numberOfEquation; j++) {
      if (coefficients[i][j] != coefficients[j][i]) {
        result = false;
        break;
      }
    }
    if (!result) {
      break;
    }
  }
  return result;
}
int wrachenie(double ** coefficients, int numberOfEquation,
  double ** solution, double precision) {
  int result = 1;
  int i, j, k;
  int maxI, maxJ;
  double max, fi;
  double ** matricaPoworota;
  matricaPoworota = new double * [numberOfEquation];
  for (i = 0; i < numberOfEquation; i++) {
    matricaPoworota[i] = new double[numberOfEquation];
  }
  double ** temp;
  temp = new double * [numberOfEquation];
  for (i = 0; i < numberOfEquation; i++) {
    temp[i] = new double[numberOfEquation];
  }
  double fault = 0.0;
  for (i = 0; i < numberOfEquation; i++) {
    for (j = i + 1; j < numberOfEquation; j++) {
      fault = fault + coefficients[i][j] * coefficients[i][j];
    }
  }
  fault = sqrt(2 * fault);
  while (fault > precision) {
    max = 0.0;
    for (i = 0; i < numberOfEquation; i++) {
      for (j = i + 1; j < numberOfEquation; j++) {
        if (coefficients[i][j] > 0 && coefficients[i][j] > max) {
          max = coefficients[i][j];
          maxI = i;
          maxJ = j;
        } else if (coefficients[i][j] < 0 && -coefficients[i][j] > max) {
          max = -coefficients[i][j];
          maxI = i;
          maxJ = j;
        }
      }
    }
    for (i = 0; i < numberOfEquation; i++) {
      for (j = 0; j < numberOfEquation; j++) {
        matricaPoworota[i][j] = 0;
      }
      matricaPoworota[i][i] = 1;
    }
    if (coefficients[maxI][maxI] == coefficients[maxJ][maxJ]) {
      matricaPoworota[maxI][maxI] = matricaPoworota[maxJ][maxJ] =
        matricaPoworota[maxJ][maxI] = sqrt(2.0) / 2.0;
      matricaPoworota[maxI][maxJ] = -sqrt(2.0) / 2.0;
    } else {
      fi = 0.5 * atan((2.0 * coefficients[maxI][maxJ]) /
        (coefficients[maxI][maxI] - coefficients[maxJ][maxJ]));
      matricaPoworota[maxI][maxI] = matricaPoworota[maxJ][maxJ] = cos(fi);
      matricaPoworota[maxI][maxJ] = -sin(fi);
      matricaPoworota[maxJ][maxI] = sin(fi);
    }
    for (i = 0; i < numberOfEquation; i++) {
      for (j = 0; j < numberOfEquation; j++) {
        temp[i][j] = 0.0;
      }
    }
    for (i = 0; i < numberOfEquation; i++) {
      for (j = 0; j < numberOfEquation; j++) {
        for (k = 0; k < numberOfEquation; k++) {
          temp[i][j] = temp[i][j] + matricaPoworota[k][i] * coefficients[k][j];
        }
      }
    }
    for (i = 0; i < numberOfEquation; i++) {
      for (j = 0; j < numberOfEquation; j++) {
        coefficients[i][j] = 0.0;
      }
    }
    for (i = 0; i < numberOfEquation; i++) {
      for (j = 0; j < numberOfEquation; j++) {
        for (k = 0; k < numberOfEquation; k++) {
          coefficients[i][j] = coefficients[i][j] +
            temp[i][k] * matricaPoworota[k][j];
        }
      }
    }
    fault = 0.0;
    for (i = 0; i < numberOfEquation; i++) {
      for (j = i + 1; j < numberOfEquation; j++) {
        fault = fault + coefficients[i][j] * coefficients[i][j];
      }
    }
    fault = sqrt(2 * fault);
    for (i = 0; i < numberOfEquation; i++) {
      for (j = 0; j < numberOfEquation; j++) {
        temp[i][j] = 0.0;
      }
    }
    for (i = 0; i < numberOfEquation; i++) {
      for (j = 0; j < numberOfEquation; j++) {
        for (k = 0; k < numberOfEquation; k++) {
          temp[i][j] = temp[i][j] + solution[i][k] * matricaPoworota[k][j];
        }
      }
    }
    for (i = 0; i < numberOfEquation; i++) {
      for (j = 0; j < numberOfEquation; j++) {
        solution[i][j] = temp[i][j];
      }
    }
    result++;
  }
  return result;
}

using namespace std;

int main() {
  setlocale(LC_ALL, "rus");
  int i, j;
  int size = 21; // N
  double ** coefficients, ** solution, precision;
  coefficients = new double * [size];
  solution = new double * [size];
  for (i = 0; i < size; i++) {
    coefficients[i] = new double[size];
    solution[i] = new double[size];
  }
  for (i = 0; i < size; i++) {
    for (j = 0; j < size; j++) {
      solution[i][j] = 0;
    }
    solution[i][i] = 1;
  }
  for (int i = 0; i < size; i++) {
    for (int j = 0; j < size; j++) {
      if (i == j) {
        coefficients[i][j] = 100;
      } else {
        coefficients[i][j] = 1 + i - j;
      }
    }
  }

  cout << endl;

  for (int i = 0; i < size; i++) {
    for (int j = 0; j < size; j++) {
      cout << coefficients[i][j] << " ";

    }
    cout << endl;
  }

  cout << endl;

  cout << "Введите точность расчета: ";
  cin >> precision;
  if (!isSimmetrial(coefficients, size)) {
    cout << "Матрица не симметричная";
  } else {
    int steps = wrachenie(coefficients, size, solution, precision);
    cout << "Reshenie:\n";
    for (i = 0; i < size; i++) {
      cout << "Собственный вектор k " <<
        i + 1 << ":\n";
      for (j = 0; j < size; j++) {
        cout << solution[j][i] << // полагаю это Х
          "\n";
      }
    }
    cout << "Собственные значения:\n";
    for (i = 0; i < size; i++) {
      cout << coefficients[i][i] <<
        "\n";
    }
    cout << "Общее число шагов: " <<
      steps;
  }
  cout << "\nPress \"Enter\" to continue..." <<
    endl;

  return 0;
}