#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <omp.h>

#define MATRIX_HEIGHT 15 * 1024
#define MATRIX_WIDTH  15 * 1024

#define EPS 0.0000001
#define TETTA 0.0001

long double countNorm(const double *vector) {
    long double result = 0;
#pragma omp parallel for reduction(+:result)
    for (int i = 0; i < MATRIX_WIDTH; i++) {
        result += (vector[i] * vector[i]);
    }
    result = sqrt(result);
    return result;
}

void printfMatrix(const double *m1) {
    for (int i = 0; i < MATRIX_WIDTH; i++) {
        printf("%lf\n",m1[i]);
    }
}

int main() {
    double *a = (double *)calloc(MATRIX_WIDTH * MATRIX_HEIGHT,sizeof(double));
    double *b = (double *)calloc(MATRIX_WIDTH,sizeof(double));

    for (int i = 0; i < MATRIX_HEIGHT; i++) {
        for (int j = 0; j < MATRIX_WIDTH; j++) {
            if (i == j) {
                a[i * MATRIX_WIDTH + j] = 2.0;
            }
            else {
                a[i * MATRIX_WIDTH + j] = 1.0;
            }
        }
        b[i] = MATRIX_WIDTH + 1;
    }

    long double b_norm = countNorm(b);

    double *x = (double *)calloc(MATRIX_WIDTH,sizeof(double));
    for (int i = 0; i < MATRIX_WIDTH; i++) {
        x[i] = 2.0;
    }

    long double result = 1.0;
    double *AxMinB = (double *)calloc(MATRIX_WIDTH,sizeof(double));
    double *AxMult = (double *)calloc(MATRIX_WIDTH,sizeof(double));
    int flag = 1;
    int counter = 2;
    while (flag) {
#pragma omp parallel for
        for (int i = 0; i < MATRIX_WIDTH; i++) {
            AxMult[i] = 0.0;
        }

#pragma omp parallel for
        for (int i = 0; i < MATRIX_HEIGHT; i++) {
            for (int j = 0; j < MATRIX_WIDTH; j++) {
                AxMult[i] += a[i*MATRIX_WIDTH + j]*x[j];
            }
            AxMinB[i] = AxMult[i] - b[i];
        }
        result = countNorm(AxMinB) / b_norm;
        if (result <= EPS) {
            flag = 0;
        }
#pragma omp parallel for
        for (int j = 0; j < MATRIX_WIDTH; j++) {
            x[j] -= AxMinB[j] * TETTA;
        }
    }

    printfMatrix(x);
    free(a);
    free(b);
    free(x);
    free(AxMinB);
    free(AxMult);

    return 0;
}