#ifdef _MSC_VER
#define _CRT_SECURE_NO_WARNINGS
#endif

#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>
#include <cmath>
#include <iostream>
#include <chrono>
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <cmath>
#define max(X,Y) ((X)>(Y)? (X):(Y))
#define min(X,Y) ((X)<(Y)? (X):(Y))
#define abs(X) ((X)>0? (X):-(X))

#define N 10000
#define M 5

using namespace std;

void multVM(double** A, double* x, double* y) {
    for (int i = 0; i < N; ++i) {
        y[i] = 0.0;
        for (int j = max(0, i - M); j <= min(N - 1, i + M); ++j) {
            y[i] += A[i][j] * x[j];
        }
    }
}

double scalarVV(double* x, double* y) {
    double z = 0;
    for (int i = 0; i < N; ++i) {
        z += x[i] * y[i];
    }
    return z;
}

void addVector(double *x, double* y, double* tmp, double modifier) {
    for (int i = 0; i < N; ++i) {
        tmp[i] = x[i] + (modifier*y[i]);
    }
}

void decVector(double *x, double* y, double* tmp, double modifier) {
    for (int i = 0; i < N; ++i) {
        tmp[i] = x[i] - (modifier*y[i]);
    }
}

int main() {
    time_t t1, t2;
    double **A = (double**)malloc(N * sizeof(double*));
    for (int i = 0; i < N; ++i) {
        A[i] = (double*)malloc(N * sizeof(double));
    }
    double x[N];
    double y[N];
    double b[N];
    double* v1 = new double[N];
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            if (abs(i - j) <= M)A[i][j] = 1 / (1.0 + abs(i - j));
            else A[i][j] = 0;
        }
        x[i] = 0.0;
        b[i] = i + 1;
    }
    //inicjalizacja
    multVM(A, x, v1);
    decVector(b, v1, v1, 1);

    double aK = 0;
    double bK = 0;
    double* r = new double[N];
    double* tmp = new double[N];
    for (int i = 0; i < N; ++i) {
        r[i] = v1[i];
    }
    double normr, normx;
    int iter = 0;
    //
    FILE* f = fopen("wyniki.dat", "w");

    std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();

    while (sqrt(scalarVV(r, r))>0.000001) {
        normr = sqrt(scalarVV(r, r));
        normx = sqrt(scalarVV(x, x));

        printf("x(k):%12lf    r(k):%12lf    ", normx, normr);
        if (f)fprintf(f, "x(k):%12lf    r(k):%12lf    ", normx, normr);
        double rrt = scalarVV(r, r);
        multVM(A, v1, tmp);
        aK = scalarVV(r, r) / scalarVV(v1, tmp);
        addVector(x, v1, x, aK);
        decVector(r, tmp, r, aK);
        bK = scalarVV(r, r) / rrt;
        printf("ak:%9lf    bk:%9lf    \n", aK, bK);
        if (f)fprintf(f, "ak:%9lf    bk:%9lf    \n", aK, bK);
        addVector(r, v1, v1, bK);
    }
    std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
    fprintf(f, "x(k) koncowe: %12lf  ", normx);

    std::cout << "norma x: " << normx <<std::endl;
    std::cout << "Time difference = " << std::chrono::duration_cast<std::chrono::microseconds>(end - begin).count() << std::endl;
    std::cout << "Time difference = " << std::chrono::duration_cast<std::chrono::nanoseconds> (end - begin).count() << std::endl;

    delete v1;
    delete r;
    delete tmp;
    for (int i = 0; i < N; ++i) {
        free(A[i]);
    }
    free(A);


    int signum;

    gsl_matrix *B = gsl_matrix_calloc(N, N);
    gsl_vector *x1 = gsl_vector_calloc(N);
    gsl_vector *b1 = gsl_vector_calloc(N);
    gsl_permutation *p = gsl_permutation_calloc(N);

    for (int i = 0; i < N; i++)
    {
        gsl_vector_set(b1, i, i + 1);
    }

    for (int i = 0; i<N; i++)
    {
        for (int j = 0; j<N; j++)
        {
            gsl_matrix_set(B, i, j, 0.0);
            if (abs(i - j) <= M)
                gsl_matrix_set(B, i, j, 1 / (1.0 + abs(i - j)));
            else
                gsl_matrix_set(B, i, j, 0.0);
        }
    }

    begin = std::chrono::steady_clock::now();

    gsl_linalg_LU_decomp(B, p, &signum);
    gsl_linalg_LU_solve(B, p, b1, x1);

    end = std::chrono::steady_clock::now();

    double x2[N];

    for (int i = 0; i < N; i++)
    {
        x2[i] = gsl_vector_get(x1, i);
    }

    double normx2 = sqrt(scalarVV(x2, x2));

    std::cout << std::endl << "norma x: " << std::endl;
    std::cout << "Time difference = " << std::chrono::duration_cast<std::chrono::microseconds>(end - begin).count() << std::endl;
    std::cout << "Time difference = " << std::chrono::duration_cast<std::chrono::nanoseconds> (end - begin).count() << std::endl;
    fprintf(f, "x(k)(LU) koncowe: %12lf\n", normx2);

    getchar();
    fclose(f);
    return 0;
}