#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include <xmmintrin.h>
#include <time.h>
#include <omp.h>
#include <iostream>
#include <ctime>
#include <iostream>
#include <Windows.h>

#define Matrix_Size  1000
#define Cache_Line 8

double** Memory_Alloc();
void Free_Memory(double** matrix);
void Rand_Matrix(double** matrix, bool dst);
bool Check_Matrix(double** matrix_A, double** matrix_B);
int Get_Cache_Size(int cache_level);
int Get_Num_Elements(int cache_size);
void Division_Matrix(int cache_size, int &block_size, int &block_size_bord);
void Multiply_default(double** matrix_A, double** matrix_B, double** matrix_Res);
void Block_Matrix_Multiply(double **matrix_A, double **matrix_B, double **matrix_Res, int sizeBlock, int block_size_bord);
void Multiply_OpenMP(double** matrix_A, double** matrix_B, double** matrix_Res);

int main() {
    unsigned __int64 time;
    int block_size_L3 = 0;
    int block_size_bord_L3 = 0;
    int block_size_L2 = 0;
    int block_size_bord_L2 = 0;
    int block_size_L1 = 0;
    int block_size_bord_L1 = 0;

    double** matrix_A = Memory_Alloc();
    double** matrix_B = Memory_Alloc();
    double** result_matrix = Memory_Alloc();
    double** result_matrix2 = Memory_Alloc();
    double** result_matrix3 = Memory_Alloc();
    double** result_matrix4 = Memory_Alloc();
    double** result_matrix5 = Memory_Alloc();

    std::srand(std::time(0));
    Rand_Matrix(matrix_A, 0);
    Rand_Matrix(matrix_B, 0);
    Rand_Matrix(result_matrix, 1);
    Rand_Matrix(result_matrix2, 1);
    Rand_Matrix(result_matrix3, 1);
    Rand_Matrix(result_matrix4, 1);
    Rand_Matrix(result_matrix5, 1);

    std::cout << "Matrix Size = " << Matrix_Size << "\n";

    int cache_size_L3 = Get_Cache_Size(3);
    int cache_size_L2 = Get_Cache_Size(2);
    int cache_size_L1 = Get_Cache_Size(1);
    std::cout << "Cache L3 size = " << cache_size_L3 << "\n";
    std::cout << "Cache L2 size = " << cache_size_L2 << "\n";
    std::cout << "Cache L1 size = " << cache_size_L1 << "\n";

    int max_matrix_size_L3 = Get_Num_Elements(cache_size_L3);
    int max_matrix_size_L2 = Get_Num_Elements(cache_size_L2);
    int max_matrix_size_L1 = Get_Num_Elements(cache_size_L1);
    std::cout << "Max matrix size L3 = " << max_matrix_size_L3 << "\n";
    std::cout << "Max matrix size L2 = " << max_matrix_size_L2 << "\n";
    std::cout << "Max matrix size L1 = " << max_matrix_size_L1 << "\n";

    //Division_Matrix(max_matrix_size_L3, block_size_L3, block_size_bord_L3);
    //std::cout << "block size L3 = " << block_size_L3 << " block size bord L3 = " << block_size_bord_L3 << "\n";
    Division_Matrix(max_matrix_size_L2, block_size_L2, block_size_bord_L2);
    std::cout << "block size L2 = " << block_size_L2 << " block size bord L2 = " << block_size_bord_L2 << "\n";
    Division_Matrix(max_matrix_size_L1, block_size_L1, block_size_bord_L1);
    std::cout << "block size L1 = " << block_size_L1 << " block size bord L1 = " << block_size_bord_L1 << "\n";

    time = __rdtsc();
    Multiply_default(matrix_A, matrix_B, result_matrix);
    std::cout << "Time mux = \t\t" << __rdtsc() - time << " ticks.\n";

    time = __rdtsc();
    Multiply_OpenMP(matrix_A, matrix_B, result_matrix3);
    std::cout << "Time OpenMP = \t\t" << __rdtsc() - time << " ticks.\n";

    //time = __rdtsc();
    //Block_Matrix_Multiply(matrix_A, matrix_B, result_matrix2, block_size_L3, block_size_bord_L3);
    //std::cout << "Time mux L3 = \t\t" << __rdtsc() - time << " ticks.\n";

    time = __rdtsc();
    Block_Matrix_Multiply(matrix_A, matrix_B, result_matrix4, block_size_L2, block_size_bord_L2);
    std::cout << "Time mux L2 = \t\t" << __rdtsc() - time << " ticks.\n";

    time = __rdtsc();
    Block_Matrix_Multiply(matrix_A, matrix_B, result_matrix5, block_size_L1, block_size_bord_L1);
    std::cout << "Time mux L1 = \t\t" << __rdtsc() - time << " ticks.\n";

    std::cout << "Matrixes 'default' and 'Block' - ";
    Check_Matrix(result_matrix, result_matrix2);
    std::cout << "Matrixes 'Block_L3' and 'OpenMP' - ";
    Check_Matrix(result_matrix2, result_matrix3);
    std::cout << "Matrixes 'BlockL2' and 'OpenMP' - ";
    Check_Matrix(result_matrix3, result_matrix4);
    std::cout << "Matrixes 'BlockL1' and 'BlockL2' - ";
    Check_Matrix(result_matrix4, result_matrix5);

    Free_Memory(matrix_A);
    Free_Memory(matrix_B);
    Free_Memory(result_matrix);
    Free_Memory(result_matrix2);
    Free_Memory(result_matrix3);
    Free_Memory(result_matrix4);
    Free_Memory(result_matrix5);

    system("pause");
    return 0;
}

double** Memory_Alloc() {
    double** matrix;

    try {
        matrix = new double*[Matrix_Size];
        for (int i = 0; i < Matrix_Size; i++) {
            matrix[i] = new double[Matrix_Size];
        }
    }
    catch (std::bad_alloc excp) {
        std::cout << "Can't allocate memory. " << std::endl;
        matrix = NULL;
    }

    return matrix;
}

void Free_Memory(double** matrix) {

    for (int i = 0; i < Matrix_Size; i++) {
        delete[] matrix[i];
    }
    delete[] matrix;
}

void Rand_Matrix(double** matrix, bool dst) {

    std::srand(std::time(0));

    for (int i = 0; i < Matrix_Size; i++) {
        for (int j = 0; j < Matrix_Size; j++) {
            if (dst) {
                matrix[i][j] = 0;
            }
            else {
                matrix[i][j] = std::rand() % 9;
            }
        }
    }
}

bool Check_Matrix(double** matrix_A, double** matrix_B) {

    for (int i = 0; i < Matrix_Size; i++) {
        for (int j = 0; j < Matrix_Size; j++) {
            if (matrix_A[i][j] != matrix_B[i][j]) {
                std::cout << "Not equal\n";
                return false;
            }
        }
    }
    std::cout << "Equal\n";
    return true;
}

int Get_Cache_Size(int cache_level) {

    int cache_size = 0;
    DWORD buffer_size = 0;
    SYSTEM_LOGICAL_PROCESSOR_INFORMATION* buffer = NULL;

    GetLogicalProcessorInformation(0, &buffer_size);
    try {
        buffer = new SYSTEM_LOGICAL_PROCESSOR_INFORMATION[buffer_size];
    }
    catch (std::bad_alloc excp) {
        std::cout << "Can't allocate memory for L3buffer. " << std::endl;
        buffer = NULL;
        return 0;
    }

    GetLogicalProcessorInformation(&buffer[0], &buffer_size);

    for (DWORD i = 0; i < buffer_size; i++) {
        if (buffer[i].Cache.Level == 3 && cache_level == 3) {
            cache_size = buffer[i].Cache.Size;
            break;
        }
        if (buffer[i].Cache.Level == 2 && cache_level == 2) {
            cache_size = buffer[i].Cache.Size;
            break;
        }
        if (buffer[i].Cache.Level == 1 && cache_level == 1) {
            cache_size = buffer[i].Cache.Size;
            break;
        }
    }

    delete[] buffer;
    return cache_size;
}

int Get_Num_Elements(int cache_size) {

    int memory = (cache_size * 0.8 / 3);        // размер кэша в байтах, процент заполнения кэша матрицай, количество равных частей кэша (3 матрицы в нашем случае)
    int double_num = memory / sizeof(double);   // количество элементов типа double
    double_num = sqrt(double_num);              // количество строк/столбцов в матрице

    while (double_num % Cache_Line != 0) {      // кратность 64 байтам
        double_num--;
    }

    return double_num;
}

void Division_Matrix(int cache_size, int &block_sizee, int &block_size_bord) {

    for (int i = cache_size; i > Cache_Line; i--) {
        if ((i % Cache_Line == 0) && (Matrix_Size % i == 0)) {
            block_sizee = i;
            block_size_bord = i;
            return;
        }
    }

    block_sizee = cache_size;
    block_size_bord = Matrix_Size;

    while (block_size_bord > block_sizee) {
        block_size_bord -= block_sizee;
    }

    return;
}

void Multiply_default(double** matrix_A, double** matrix_B, double** matrix_Res) {

    for (int i = 0; i < Matrix_Size; i++) {
        for (int j = 0; j < Matrix_Size; j++) {
            for (int k = 0; k < Matrix_Size; k++) {
                matrix_Res[i][j] += matrix_A[i][k] * matrix_B[k][j];
            }
        }
    }
}

void Block_Matrix_Multiply(double **matrix_A, double **matrix_B, double **matrix_Res, int sizeBlock, int block_size_bord) {

    int count = Matrix_Size / sizeBlock;
    if (Matrix_Size % sizeBlock != 0) {
        count++;
    }

    for (int i = 0, sizeIBlock = sizeBlock; i < count; i++) {
        if (i == count - 1)
            sizeIBlock = block_size_bord;
        for (int j = 0, sizeJBlock = sizeBlock; j < count; j++) {
            if (j == count - 1)
                sizeJBlock = block_size_bord;
            for (int k = 0, sizeKBlock = sizeBlock; k < count; k++) {
                if (k == count - 1) {
                    sizeKBlock = block_size_bord;
                }

                for (int l = 0; l < sizeIBlock; l++) {
                    for (int n = 0; n < sizeKBlock; n++) {
                        for (int m = 0; m < sizeJBlock; m++) {
                            matrix_Res[i * sizeBlock + l][j * sizeBlock + m] +=
                                matrix_A[i * sizeBlock + l][k * sizeBlock + n] *
                                matrix_B[k * sizeBlock + n][j * sizeBlock + m];
                        }
                    }
                }
            }
        }
    }
}

void Multiply_OpenMP(double** matrix_A, double** matrix_B, double** matrix_Res) {

#pragma omp parallel for
    for (int i = 0; i < Matrix_Size; i++) {
        for (int k = 0; k < Matrix_Size; k++) {
            for (int j = 0; j < Matrix_Size; j++) {
                matrix_Res[i][j] += matrix_A[i][k] * matrix_B[k][j];
            }
        }
    }
}