#include "pch.h"
#include <iostream>
#include <omp.h>
#include <cstdlib>
#include <ctime>
#include <time.h>
#include <stdio.h>

using namespace std;

#define N 600
#define R 100
#define NUM_OF_THREADS 4

FILE *result_file;

float matrix_a[N][N];
float matrix_b[N][N];
float matrix_r[N][N];
auto threeForTempMatrix = new float[10][N][N];

void initialize_matrices()
{
    // zdefiniowanie zawarosci poczatkowej macierzy
    //#pragma omp parallel for
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            matrix_a[i][j] = (float)rand() / RAND_MAX;
            matrix_b[i][j] = (float)rand() / RAND_MAX;
            matrix_r[i][j] = 0.0;
        }
    }
}

void initialize_matricesZ()
{
    // zdefiniowanie zawarosci poczatkowej macierzy
#pragma omp parallel for
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            matrix_r[i][j] = 0.0;
        }
    }
}

void initialize_temp()
{
    // zdefiniowanie zawarosci poczatkowej macierzy
#pragma omp parallel for
    for (int k = 0; k < 6; k++) {
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                threeForTempMatrix[k][i][j] = 0.0;
            }
        }
    }
}

void threeLoops() {
#pragma omp parallel
    {
        int num = omp_get_thread_num();
#pragma omp for nowait
        for (int k = 0; k < N; k++) {
            for (int i = 0; i < N; i++) {
                for (int j = 0; j < N; j++) {
                    threeForTempMatrix[num][i][j] += matrix_a[i][k] * matrix_b[k][j];
                }
            }
        }

        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
#pragma omp atomic
                matrix_r[i][j] += threeForTempMatrix[num][i][j];
            }
        }
    }

}

void sixLoops() {
    for (int i = 0; i < N; i += R)
#pragma omp parallel for
        for (int j = 0; j < N; j += R)
            for (int k = 0; k < N; k += R)
                for (int ii = i; ii < i + R; ii++)
                    for (int kk = k; kk < k + R; kk++)
                        for (int jj = j; jj < j + R; jj++)
                            matrix_r[ii][jj] += matrix_a[ii][kk] * matrix_b[kk][jj];
}

void sequenceLoopsSmall() {
    for (int i = 0; i < N; i++)
        for (int k = 0; k < N; k++)
            for (int j = 0; j < N; j++)
                matrix_r[i][j] += matrix_a[i][k] * matrix_b[k][j];
}

void sequenceLoopsBig() {
    for (int i = 0; i < N; i += R)
        for (int j = 0; j < N; j += R)
            for (int k = 0; k < N; k += R)
                for (int ii = i; ii < i + R; ii++)
                    for (int kk = k; kk < k + R; kk++)
                        for (int jj = j; jj < j + R; jj++)
                            matrix_r[ii][jj] += matrix_a[ii][kk] * matrix_b[kk][jj];
}

int main() {
    srand(time(NULL));

    if (fopen_s(&result_file, "results.txt", "a") < 0) {
        fprintf(stderr, "nie mozna otworzyc pliku wyniku \n");
        perror("results");
        return(EXIT_FAILURE);
    }

    omp_set_num_threads(NUM_OF_THREADS);

    fprintf(result_file, "N for 3-LOOP: %d\n", N);
    fprintf(result_file, "N for 6-LOOP: %d\n", N);
    fprintf(result_file, "R for 6-LOOP: %d\n", R);
    fprintf(result_file, "NUM_OF_THREADS : %d\n", NUM_OF_THREADS);

    initialize_matrices();
    cout << "Matrices are prepared..." << endl;

    initialize_matricesZ();
    double start = (double)clock() / CLK_TCK;
    sequenceLoopsSmall();
    std::cout << "\n" << matrix_r[0][0] << "\n";
    double end = (double)clock() / CLK_TCK;
    printf("3-loops sequence processing: %8.4f sec \n", end - start);
    fprintf(result_file,
        "3-loops sequence processing : %8.4f sec (%6.4f sec rozdzielczosc pomiaru)\n",
        end - start, 1.0 / CLK_TCK);

    initialize_matricesZ();
    start = (double)clock() / CLK_TCK;
    sequenceLoopsBig();
    std::cout << "\n" << matrix_r[0][0] << "\n";
    end = (double)clock() / CLK_TCK;
    printf("6-loops sequence processing: %8.4f sec \n", end - start);
    fprintf(result_file,
        "6-loops sequence processing : %8.4f sec (%6.4f sec rozdzielczosc pomiaru)\n",
        end - start, 1.0 / CLK_TCK);

    initialize_matricesZ();
    initialize_temp();
    start = (double)clock() / CLK_TCK;
    threeLoops();
    std::cout << "\n" << matrix_r[0][0] << "\n";
    end = (double)clock() / CLK_TCK;
    printf("3-loops parallel processing: %8.4f sec \n", end - start);
    fprintf(result_file,
        "3-loops parallel processing : %8.4f sec (%6.4f sec rozdzielczosc pomiaru)\n",
        end - start, 1.0 / CLK_TCK);

    initialize_matricesZ();
    start = (double)clock() / CLK_TCK;
    sixLoops();
    std::cout << "\n" << matrix_r[0][0] << "\n";
    end = (double)clock() / CLK_TCK;
    printf("6-loops parallel processing: %8.4f sec \n", end - start);
    fprintf(result_file,
        "6-loops parallel processing : %8.4f sec (%6.4f sec rozdzielczosc pomiaru)\n",
        end - start, 1.0 / CLK_TCK);

    fprintf(result_file, "===========================================\n");
    char c;
    std::cin >> c;
    fclose(result_file);
    return 0;
}