Метод Гаусса для решения СЛАУ

#include <iostream>
#include <math.h>
#include <iomanip>
using namespace std;
 
int main ()
{
    const int n = 5;
    double x[n], y[n], A[3][3], b[3], p[3], ai0, a21;
    int i, j, k;
 
    cout << "Vvedite pari x, y:\n";
    for (i = 0; i < n; i++)
        cin >> x[i] >> y[i];
 
    for (i = 0; i < 3; i++)
        for (j = 0; j < 3; j++) {
            A[i][j] = 0;
            for (k = 0; k < n; k++)
                A[i][j] += pow(x[k], i+j);
        }
 
    for (i = 0; i < 3; i++) {
        b[i] = 0;
        for (j = 0; j < n; j++)
            b[i] += y[j]*pow(x[j],i);
    }
 
    cout << "\nMatrica A:\n";
    for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++)
            cout << setw(16) << A[i][j];
        cout << setw(16) << b[i];
        cout << endl;
    }
 
    for (i = 1; i < 3; i++) {
        ai0 = A[i][0];
        for (j = 0; j < 3; j++)
            A[i][j] -= ai0*A[0][j]/A[0][0];
        b[i] -= ai0*b[0]/A[0][0];
    }
 
    a21 = A[2][1];
    A[2][1] -= a21*A[1][1]/A[1][1];
    A[2][2] -= a21*A[1][2]/A[1][1];
    b[2] -= a21*b[1]/A[1][1];
 
    p[2] = b[2]/A[2][2];
    p[1] = (b[1]-A[1][2]*x[2])/A[1][1];
    p[0] = (b[0]-A[0][2]*x[2]-A[0][1]*x[1])/A[0][0];
 
    cout << "\nY = " << p[2] << "*X^2  +  " << p[1] << "*X  +  " << p[0] << endl;
 
    return 0;
}

// файл data.h, шаблонный класс data, C++
#include <vector>
using namespace std;       // пространство имён std
typedef unsigned int uint; // беззнаковый целый
template<class Type>
class data {
public: // открытые свойства и методы
  data(uint ms) {   // конструктор
    try {
      msize = ms;
      vecB.assign(msize, (Type)0);
      vecX.assign(msize, (Type)0);
      matA.resize(msize);
      for (uint i=0; i<msize; i++) matA.at(i).assign(msize, (Type)0);
    } catch (...) {
      msize = 0;
      throw;
    }
  }
  ~data() {} // деструктор
  // методы для манипулирования данными
  uint size() {return (msize);}
  void savA(uint j, uint i, Type val) {matA.at(j).at(i) = val;}
  void savB(uint j, Type val) {vecB.at(j) = val;}
  Type readA(uint j, uint i) {return (matA.at(j).at(i));}
  Type readB(uint j) {return (vecB.at(j));}
  Type readX(uint j) {return (vecB.at(j));}
protected: // защищённые свойства и методы
  typedef vector<Type> vType; // тип вектор типа Type
  typedef vector<Type>::iterator pvType; // итератор для вектора типа Type
  vector<vType> matA; // матрица коэффициентов A
  vType vecB; // вектор правых частей B
  vType vecX; // вектор искомый X
  uint msize; // размер матриц и векторов
};

// листинг C++ кода (файл algorithm.h), шаблонный класс algorithm, C++
#define DEBUG true
#include <iostream>
#include <math.h>
#include "data.h"
// класс некоторых действий линейной алгебры
template<class Type> class algorithm : public data<Type> {
public: // открытые свойства и методы
  algorithm(uint ms) : data<Type>(ms) {} // конструктор
  bigelem() { // формирование системы с ведущим элементом
    uint k, i, j, kmax;
    Type maxval;
    for(k = 0; k < size(); k++) { 
      maxval = matA.at(k).at(k); // проход по диагонали
      kmax = k;
      for(j = k + 1; j < size(); j++) { // больший элемент в столбце
        if(abs(matA.at(j).at(k)) > abs(maxval)) {
          maxval = matA.at(j).at(k);
          kmax = j;
        }
      }
      if(kmax != k) { // если требуется перестановка
        for(i = 0; i < size(); i++) { // перестановка строк матрицы A
          maxval = matA.at(k).at(i);
          matA.at(k).at(i) = matA.at(kmax).at(i);
          matA.at(kmax).at(i) = maxval;
        }
        maxval = vecB.at(k); // перестановка элементов вектора B
        vecB.at(k) = vecB.at(kmax);
        vecB.at(kmax) = maxval;
      }
    }
  }
 
gauss() {  // метод исключения Гаусса
    int k, i, j;
    Type temp;
    for(k = 1; k < size(); k++) { // прямой ход
      for(i = k; i < size(); i++) {
        temp = matA.at(i).at(k - 1) / matA.at(k - 1).at(k - 1);
        for(j = 0; j < size(); j++) {
          matA.at(i).at(j) =
               matA.at(i).at(j) - temp * matA.at(k - 1).at(j);
        }
        vecB.at(i) = vecB.at(i) - temp * vecB.at(k - 1);
      }
#if (DEBUG == true)
      cout << "Forward stroke, cycle: " << k << endl;
      prnls();
#endif
    }
#if (DEBUG == true)
    cout<<"Reverse"<<endl;
#endif
    for(k = size() - 1; k >= 0; k--) { // обратный ход
      temp = vecB.at(k);
      for(j = size() - 1; j > k; j--) {
        temp -= matA.at(k).at(j) * vecX.at(j);
      }
      vecX.at(k) = temp / matA.at(k).at(k);
#if (DEBUG == true)
      cout<<"X["<<k<<"] = "<<vecX.at(k)<<endl;
#endif
    }
  }
  void prnls() { // вывод системы в поток экрана
    int j, i;
    for(j = 0; j < size(); j++) {
      for(i = 0; i < size(); i++) {
        if(matA.at(j).at(i) < 0)
          cout<<matA.at(j).at(i)<<"*X["<<i <<"]";
        else
          cout<<"+"<<matA.at(j).at(i)<<"*X[" <<i<<"]";
      }
      cout<<"="<<vecB.at(j)<<endl;
    }
    return;
  }
};

// листинг C++ кода проверки классов data и algorithm
#include <iostream>
#include <stdlib.h>
#include "algorithm.h"
using namespace std;
int main(int argc, char* argv[]) {
  cout << "Size of the system of equations:  ";
  int sz;
  float val;
  cin>>sz;  // порядок системы
  // объект класса шаблонного algorithm
  algorithm<float> *la = new algorithm<float>(sz);
  // ввод системы линейных уравнений
  cout<<"Enter lines of system of linear equations:  "<<endl;
  for (int j = 0; j < sz; j++) {
    cout<<"A["<<j<<",1]+...+A["<<j<<", "
        <<(sz-1)<<"]=B["<<j<<"]"<<endl;
    for (int i = 0; i < sz; i++) {
      cin>>val;
      la->savA(j, i, val);
    }
    cin>>val;
    la->savB(j, val);
  }
  cout<<"Linear system:"<<endl;
  la->prnls();
  la->bigelem();
  cout<<"Linear system:"<<endl;
  la->prnls();
  la->gauss();
  delete la;
  system("pause");
  return 0;
}