eeompst.cpp

#include <omp.h>
#include <bits/stdc++.h>
// random function is still needed to get improvement.
using namespace std;

void print(int n, double** mat)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            cout << mat[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}

double** allocateMatrix(int n)
{
    double* data = (double*)malloc(n * n * sizeof(double));
    double** array = (double**)malloc(n * sizeof(double*));
    for (int i = 0; i < n; i++)
    {
        array[i] = &(data[n * i]);
    }
    return array;
}

void fillMatrix(int n, double**& mat, int num)
{
    srand(num);
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            // mat[i][j] = rand() % 5;
            int t = rand() % 2 + 1;
            mat[i][j] = pow(-1, t) * rand() / float(RAND_MAX);
        }
    }
}

void freeMatrix(int n, double** mat)
{
    free(mat[0]);
    free(mat);
}

double** naive(int n, double** mat1, double** mat2)
{
    double** prod = allocateMatrix(n);

    int i, j;

#pragma omp parallel for collapse(2)
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            prod[i][j] = 0;
            for (int k = 0; k < n; k++)
            {
                prod[i][j] += mat1[i][k] * mat2[k][j];
            }
        }
    }

    return prod;
}

double** getSlice(int n, double** mat, int offseti, int offsetj)
{
    int m = n / 2;
    double** slice = allocateMatrix(m);
    #pragma omp parallel for collapse(2)
    for (int i = 0; i < m; i++)
    {
        for (int j = 0; j < m; j++)
        {
            slice[i][j] = mat[offseti + i][offsetj + j];
        }
    }
    return slice;
}

double** addMatrices(int n, double** mat1, double** mat2, bool add)
{
    double** result = allocateMatrix(n);
    #pragma omp parallel for collapse(2)
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (add)
                result[i][j] = mat1[i][j] + mat2[i][j];
            else
                result[i][j] = mat1[i][j] - mat2[i][j];
        }
    }

    return result;
}

double** combineMatrices(int m, double** c11, double** c12, double** c21, double** c22)
{
    int n = 2 * m;
    double** result = allocateMatrix(n);
    #pragma omp parallel for collapse(2)
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (i < m && j < m)
                result[i][j] = c11[i][j];
            else if (i < m)
                result[i][j] = c12[i][j - m];
            else if (j < m)
                result[i][j] = c21[i - m][j];
            else
                result[i][j] = c22[i - m][j - m];
        }
    }

    return result;
}

double** strassen(int n, double** mat1, double** mat2)
{

    if (n <= 32  || n % 2 != 0)
    {
        return naive(n, mat1, mat2);
    }

    int m = n / 2;

    double** a = getSlice(n, mat1, 0, 0);
    double** b = getSlice(n, mat1, 0, m);
    double** c = getSlice(n, mat1, m, 0);
    double** d = getSlice(n, mat1, m, m);
    double** e = getSlice(n, mat2, 0, 0);
    double** f = getSlice(n, mat2, 0, m);
    double** g = getSlice(n, mat2, m, 0);
    double** h = getSlice(n, mat2, m, m);

    // double** gha;
    // #pragma omp task shared(gha)
    // {
    //     gha = addMatrices(m, g, h, true);
    // }

    double** s1;
    #pragma omp task shared(s1)
    {
        double** bds = addMatrices(m, b, d, false);
        double** gha = addMatrices(m, g, h, true);
        s1 = strassen(m, bds, gha);
        freeMatrix(m, bds);
        freeMatrix(m, gha);
    }

    double** s2;
    #pragma omp task shared(s2)
    {
        double** ada = addMatrices(m, a, d, true);
        double** eha = addMatrices(m, e, h, true);
        s2 = strassen(m, ada, eha);
        freeMatrix(m, ada);
        freeMatrix(m, eha);
    }

    double** s3;
    #pragma omp task shared(s3)
    {
        double** acs = addMatrices(m, a, c, false);
        double** efa = addMatrices(m, e, f, true);
        s3 = strassen(m, acs, efa);
        freeMatrix(m, acs);
        freeMatrix(m, efa);
    }

    double** s4;
    #pragma omp task shared(s4)
    {
        double** aba = addMatrices(m, a, b, true);
        s4 = strassen(m, aba, h);
        freeMatrix(m, aba);
    }

    double** s5;
    #pragma omp task shared(s5)
    {
        double** fhs = addMatrices(m, f, h, false);
        s5 = strassen(m, a, fhs);
        freeMatrix(m, fhs);
    }

    double** s6;
    #pragma omp task shared(s6)
    {
        double** ges = addMatrices(m, g, e, false);
        s6 = strassen(m, d, ges);
        freeMatrix(m, ges);
    }

    double** s7;
    #pragma omp task shared(s7)
    {
        double** cda = addMatrices(m, c, d, true);
        s7 = strassen(m, cda, e);
        freeMatrix(m, cda);
    }

    #pragma omp taskwait

    freeMatrix(m, a);
    freeMatrix(m, b);
    freeMatrix(m, c);
    freeMatrix(m, d);
    freeMatrix(m, e);
    freeMatrix(m, f);
    freeMatrix(m, g);
    freeMatrix(m, h);

    double** c11;
    #pragma omp task shared(c11)
    {
        double** s1s2a = addMatrices(m, s1, s2, true);
        double** s6s4s = addMatrices(m, s6, s4, false);
        c11 = addMatrices(m, s1s2a, s6s4s, true);
        freeMatrix(m, s1s2a);
        freeMatrix(m, s6s4s);
    }

    double** c12;
    #pragma omp task shared(c12)
    {
        c12 = addMatrices(m, s4, s5, true);
    }

    double** c21;
    #pragma omp task shared(c21)
    {
        c21 = addMatrices(m, s6, s7, true);
    }

    double** c22;
    #pragma omp task shared(c22)
    {
        double** s2s3s = addMatrices(m, s2, s3, false);
        double** s5s7s = addMatrices(m, s5, s7, false);
        c22 = addMatrices(m, s2s3s, s5s7s, true);
        freeMatrix(m, s2s3s);
        freeMatrix(m, s5s7s);
    }

    #pragma omp taskwait

    freeMatrix(m, s1);
    freeMatrix(m, s2);
    freeMatrix(m, s3);
    freeMatrix(m, s4);
    freeMatrix(m, s5);
    freeMatrix(m, s6);
    freeMatrix(m, s7);

    double** prod = combineMatrices(m, c11, c12, c21, c22);

    freeMatrix(m, c11);
    freeMatrix(m, c12);
    freeMatrix(m, c21);
    freeMatrix(m, c22);

    return prod;
}

bool check(int n, double** prod1, double** prod2)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (prod1[i][j] != prod2[i][j])
                return false;
        }
    }
    return true;
}

int main()
{
	int num[57] = {32, 64, 128, 256, 512, 1024, 2048, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000, 3100, 3200, 3300, 3400, 3500, 3600, 3700, 3800, 3900, 4000, 4100, 4200, 4300, 4400, 4500, 4600, 4700, 4800, 4900, 5000};
    fstream ff;
    // int n;
    // cout << "\nEnter matrix Dimension: ";
    // cin >> n;
    ff.open("eeompst.txt", ios::out);
    ff << "OMP -- ST" << endl;
    ff.close();
    ff.open("eeompst.txt", ios::out|ios::app);

    for (int i = 0; i < 57; i++) {
        double** mat1 = allocateMatrix(num[i]);
        fillMatrix(num[i], mat1, 1017);

        double** mat2 = allocateMatrix(num[i]);
        fillMatrix(num[i], mat2, 1001);

        double startParStrassen = omp_get_wtime();
        double** prod;

        omp_set_num_threads(4);

        #pragma omp parallel
        {
        #pragma omp single
            {
                prod = strassen(num[i], mat1, mat2);
            }
        }
        double endParStrassen = omp_get_wtime();
        // cout << "\nParallel Strassen Runtime (OMP): ";
        cout << (endParStrassen - startParStrassen)*1000 << endl;
        ff << (endParStrassen - startParStrassen)*1000 << endl;
        // cout << endl;
    }

    return 0;
}