Richardson iterative using Chebyshev polynomials

#include <iostream>
#include <fstream>
#include <iomanip>

#define sign(x) (x < 0 ? (-1) : (1))
#define eps 0.00000001

using namespace std;
ofstream fout;
ifstream fin;

double dotProduct(double **x, double **y, int n) {
    double product = 0;
    for (int i = 0; i < n; i++)
        product += x[i][0] * y[i][0];
    return product;
}

double **matrixMultL(double **A, double **B, int row, int column, int n) {
    int i, j;
    double **product = new double*[row]();
    for (i = 0; i < row; i++) {
        product[i] = new double[column]();
        for (j = 0; j < column; j++) {
            for (int inner = 0; inner < n; inner++) {
                product[i][j] += A[i][inner] * B[inner][j];
            }
        }
    }
    //for (i = 0; i < row; i++) for (j = 0; j < column; j++) if (abs(product[i][j]) <= eps) product[i][j] = 0;
    return product;
}

double **matrixTranspose(double **M, int row, int column) {
    double **product = new double*[column];
    for (int i = 0; i < column; i++) product[i] = new double[row];
    for (int i = 0; i < row; i++) for (int j = 0; j < column; j++) product[j][i] = M[i][j];
    return product;
}

double **squareMatrixSub(double **A, double **B, int n) {
    double **product = new double*[n]();
    for (int i = 0; i < n; i++) {
        product[i] = new double[n];
        for (int j = 0; j < n; j++) product[i][j] = A[i][j] - B[i][j];
    }
    return product;
}

double **vectorSub(double **x, double **y, int n) {
    double ** product = new double*[n];
    for (int i = 0; i < n; i++) {
        product[i] = new double[1];
        product[i][0] = x[i][0] - y[i][0];
    }
    return product;
}

double **multByScalar(double **M, int row, int column, double a) {
    for (int i = 0; i < row; i++)
    for (int j = 0; j < column; ++j)
        M[i][j] *= a;
    return M;
}

void rowSwap(double **M, int row1, int row2, int n) {
    int i, j, t;
    for (j = 0; j < n; j++) {
        t = M[row1][j];
        M[row1][j] = M[row2][j];
        M[row2][j] = t;
    }
}

int check(double **M, int row, int column, int n) {
    if (M[row][column] != 0) return 1;
    else return 0;
}

double **richardsonIterative(double **M, double **V, double **X0, double **X1, double w1, double w2, int n, double a, double b, int *cnt) {
    double **X = new double*[n], **Y = new double*[n], **product = new double*[n], u;
    for (int i = 0; i < n; i++) {
        X[i] = new double[1];
        Y[i] = new double[1];
    }
    u = 1 / (2 * (1 / w1) - w2);
    X = vectorSub(vectorSub(X1, multByScalar(vectorSub(X1, X0, n), n, 1, -u * w2), n), multByScalar(vectorSub(matrixMultL(M, X1, n, 1, n), V, n), n, 1, (1 + u * w2) * 2 / (b + a)), n);
    product = vectorSub(X, X1, n);
    if (sqrt(dotProduct(product, product, n)) <= eps) return X1;
    else {
        (*cnt)++;
        X = richardsonIterative(M, V, X1, X, w1, u, n, a, b, cnt);
    }
    return X;
}


int main() {
    fin.open("input.txt");
    fout.open("output.txt");
    int n, cntTmp = 0;
    int *cnt = &cntTmp;
    double a, b, u1, u2, w1, w2;
    fin >> a >> b;
    fin >> n;
    double **X = new double*[n](), **M = new double*[n], **V = new double*[n], **X0 = new double*[n], **Y = new double*[n];
    for (int i = 0; i < n; i++) {
        M[i] = new double[n + 1];
        V[i] = new double[1];
        X0[i] = new double[1]();
        X[i] = new double[1];
        for (int j = 0; j <= n; j++) {
            X0[i][0] = 1;
            if (j == n) fin >> V[i][0];
            else fin >> M[i][j];
        }
    }
    w1 = -(b - a) / (b + a);
    w2 = 1 / (2 * (1 / w1) - w1);
    X = vectorSub(X0, multByScalar((vectorSub(matrixMultL(M, X0, n, 1, n), V, n)), n, 1, 2 / (a + b)), n);
    Y = vectorSub(vectorSub(X, multByScalar(vectorSub(X, X0, n), n, 1, w1 * w2), n), multByScalar(vectorSub(matrixMultL(M, X, n, 1, n), V, n), n, 1, (1 + w1 * w2) * 2 / (b + a)), n);
    X = richardsonIterative(M, V, X, Y, w1, w2, n, a, b, cnt);
    fout << "Ax = b solution set:" << endl;
    for (int i = 0; i < n; i++) fout << "x[" << i << "] = " << setprecision(10) << X[i][0] << endl;
    fout << "Number of iterations:" << endl << *cnt;
}