eempi4.cpp

#include <mpi.h>
#include <iomanip>
#include <iostream>
#include <cmath>
#include <fstream>
// NOT OK, for raspberrypi 4b, however, now is fine!!! Happy!!!
// Ready for testing

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++)
        {
            int t = rand() % 2 + 1;
            mat[i][j] = pow(-1, t) * rand() / float(RAND_MAX);
            // mat[i][j] = rand() % 10;
        }
    }
}

void freeMatrix(int n, double** mat) // 释放矩阵的内存使用
{
    free(mat[0]);
    free(mat);
}

double** naive(int n, double** mat1, double** mat2) // 传统矩阵乘法算法
{
    double** prod = allocateMatrix(n);

    for (int i = 0; i < n; i++)
    {
        for (int 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);
    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);
    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);

    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** moreZeros(int m, double** mat1) { // pad 0
    double** mat11 = allocateMatrix(m+1);
    if (true) {
        // cout << "Hello**" << endl;
        for (int i = 0; i < m + 1; i++) {
            for (int j = 0; j < m + 1; j++) {
                if (i < m && j < m) {
                    mat11[i][j] = mat1[i][j];
                    // cout << mat1[i][j] << endl;
                } else {
                    mat11[i][j] = 0;
                }
            }
        }
    }
    // print(m+1, mat11);

    return mat11;
}

double** strassen(int n, double** mat1, double** mat2) // Strassen 32以内回归的方式
{

    if (n <= 64 || n % 2 != 0) // based on testing, 64 is the best,
    {
        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** bds = allocateMatrix(m);
    double** gha = allocateMatrix(m);
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
            bds[i][j] = b[i][j] - d[i][j];
            gha[i][j] = g[i][j] + h[i][j];
        }
    }
    double** s1 = strassen(m, bds, gha);
    freeMatrix(m, bds);
    freeMatrix(m, gha);

    double** ada = allocateMatrix(m);
    double** eha = allocateMatrix(m);
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
            ada[i][j] = a[i][j] + d[i][j];
            eha[i][j] = e[i][j] + h[i][j];
        }
    }
    double** s2 = strassen(m, ada, eha);
    freeMatrix(m, ada);
    freeMatrix(m, eha);

    double** acs = allocateMatrix(m);
    double** efa = allocateMatrix(m);
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
            acs[i][j] = a[i][j] - c[i][j];
            efa[i][j] = e[i][j] + f[i][j];
        }
    }
    double** s3 = strassen(m, acs, efa);
    freeMatrix(m, acs);
    freeMatrix(m, efa);

    double** aba = addMatrices(m, a, b, true);
    double** s4 = strassen(m, aba, h);
    freeMatrix(m, aba);
    freeMatrix(m, b);

    double** fhs = addMatrices(m, f, h, false);
    double** s5 = strassen(m, a, fhs);
    freeMatrix(m, fhs);
    freeMatrix(m, a);
    freeMatrix(m, f);
    freeMatrix(m, h);

    double** ges = addMatrices(m, g, e, false);
    double** s6 = strassen(m, d, ges);
    freeMatrix(m, ges);
    freeMatrix(m, g);

    double** cda = addMatrices(m, c, d, true);
    double** s7 = strassen(m, cda, e);
    freeMatrix(m, cda);
    freeMatrix(m, c);
    freeMatrix(m, d);
    freeMatrix(m, e);

    double** s1s2a = addMatrices(m, s1, s2, true);
    double** s6s4s = addMatrices(m, s6, s4, false);
    double** c11 = addMatrices(m, s1s2a, s6s4s, true);
    freeMatrix(m, s1s2a);
    freeMatrix(m, s6s4s);
    freeMatrix(m, s1);

    double** c12 = addMatrices(m, s4, s5, true);
    freeMatrix(m, s4);

    double** c21 = addMatrices(m, s6, s7, true);
    freeMatrix(m, s6);

    double** s2s3s = addMatrices(m, s2, s3, false);
    double** s5s7s = addMatrices(m, s5, s7, false);
    double** c22 = addMatrices(m, s2s3s, s5s7s, true);
    freeMatrix(m, s2s3s);
    freeMatrix(m, s5s7s);
    freeMatrix(m, s2);
    freeMatrix(m, s3);
    freeMatrix(m, s5);
    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;
}

void strassen(int n, double** mat1, double** mat2, double**& prod, int rank) // Strassen 不回归的方式
{

    if (n == 1)
    {
        prod = allocateMatrix(1);
        prod[0][0] = mat1[0][0] * mat2[0][0];
    }

    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** s1 = allocateMatrix(m);
    double** s2 = allocateMatrix(m);
    double** s3 = allocateMatrix(m);
    double** s4 = allocateMatrix(m);
    double** s5 = allocateMatrix(m);
    double** s6 = allocateMatrix(m);
    double** s7 = allocateMatrix(m);

    if (rank == 0)
    {
        MPI_Recv(&(s1[0][0]), m * m, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        MPI_Recv(&(s2[0][0]), m * m, MPI_DOUBLE, 2, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        MPI_Recv(&(s3[0][0]), m * m, MPI_DOUBLE, 3, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        MPI_Recv(&(s4[0][0]), m * m, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        MPI_Recv(&(s5[0][0]), m * m, MPI_DOUBLE, 2, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        MPI_Recv(&(s6[0][0]), m * m, MPI_DOUBLE, 3, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
        MPI_Recv(&(s7[0][0]), m * m, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    }

    if (rank == 1)
    {
        double** bds = allocateMatrix(m);
        double** gha = allocateMatrix(m);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < m; j++) {
                bds[i][j] = b[i][j] - d[i][j];
                gha[i][j] = g[i][j] + h[i][j];
            }
        }
        s1 = strassen(m, bds, gha);
        freeMatrix(m, bds);
        freeMatrix(m, gha);
        MPI_Send(&(s1[0][0]), m * m, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    }

    if (rank == 2)
    {
        double** ada = allocateMatrix(m);
        double** eha = allocateMatrix(m);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < m; j++) {
                ada[i][j] = a[i][j] + d[i][j];
                eha[i][j] = e[i][j] + h[i][j];
            }
        }
        s2 = strassen(m, ada, eha);
        freeMatrix(m, ada);
        freeMatrix(m, eha);
        MPI_Send(&(s2[0][0]), m * m, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    }

    if (rank == 3)
    {
        double** acs = allocateMatrix(m);
        double** efa = allocateMatrix(m);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < m; j++) {
                acs[i][j] = a[i][j] - c[i][j];
                efa[i][j] = e[i][j] + f[i][j];
            }
        }
        s3 = strassen(m, acs, efa);
        freeMatrix(m, acs);
        freeMatrix(m, efa);
        MPI_Send(&(s3[0][0]), m * m, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    }

    if (rank == 1)
    {
        double** aba = addMatrices(m, a, b, true);
        s4 = strassen(m, aba, h);
        freeMatrix(m, aba);
        MPI_Send(&(s4[0][0]), m * m, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    }

    if (rank == 2)
    {
        double** fhs = addMatrices(m, f, h, false);
        s5 = strassen(m, a, fhs);
        freeMatrix(m, fhs);
        MPI_Send(&(s5[0][0]), m * m, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    }

    if (rank == 3)
    {
        double** ges = addMatrices(m, g, e, false);
        s6 = strassen(m, d, ges);
        freeMatrix(m, ges);
        MPI_Send(&(s6[0][0]), m * m, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    }


    if (rank == 1)
    {
        double** cda = addMatrices(m, c, d, true);
        s7 = strassen(m, cda, e);
        freeMatrix(m, cda);
        MPI_Send(&(s7[0][0]), m * m, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    }
    freeMatrix(m, b);
    freeMatrix(m, a);
    freeMatrix(m, f);
    freeMatrix(m, h);
    freeMatrix(m, g);
    freeMatrix(m, c);
    freeMatrix(m, d);
    freeMatrix(m, e);

    MPI_Barrier(MPI_COMM_WORLD);

    if (rank == 0)
    {
        double** s1s2a = addMatrices(m, s1, s2, true);
        double** s6s4s = addMatrices(m, s6, s4, false);
        double** c11 = addMatrices(m, s1s2a, s6s4s, true);
        freeMatrix(m, s1s2a);
        freeMatrix(m, s6s4s);

        double** c12 = addMatrices(m, s4, s5, true);

        double** c21 = addMatrices(m, s6, s7, true);

        double** s2s3s = addMatrices(m, s2, s3, false);
        double** s5s7s = addMatrices(m, s5, s7, false);
        double** c22 = addMatrices(m, s2s3s, s5s7s, true);
        freeMatrix(m, s2s3s);
        freeMatrix(m, s5s7s);

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

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

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

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 argc, char* argv[])
{
    int p_rank;
    int num_process;
	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};
	int randomNum[5] = {1001, 8808, 9909, 1242, 1017};
	int randomNum1[5] = {1243, 7838, 1893, 1283, 1893};
    fstream ff;
    ff.open("eempi4.txt", ios::out);
    ff << "MPI -- ST" << endl;
    ff.close();
    ff.open("eempi4.txt", ios::out|ios::app);

    if (MPI_Init(&argc, &argv) != MPI_SUCCESS)
    {
        printf("MPI-INIT Failed\n");
        return 0;
    }

    MPI_Comm_rank(MPI_COMM_WORLD, &p_rank);
    MPI_Comm_size(MPI_COMM_WORLD, &num_process);


    for (int i = 0; i < 57; i++) {
        if (num[i] % 2 == 0) {
            MPI_Barrier(MPI_COMM_WORLD);
            MPI_Bcast(&num[i], 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

            double** mat1 = allocateMatrix(num[i]);
            double** mat2 = allocateMatrix(num[i]);
            double** prod = allocateMatrix(num[i]);
            // cout << p_rank  << endl;

            // for (int j = 0; j < 5; i++) {
            if (p_rank == 0)
            {
                fillMatrix(num[i], mat1, 1017);
                fillMatrix(num[i], mat2, 1000);
            }

            MPI_Bcast(&(mat1[0][0]), num[i] * num[i], MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Bcast(&(mat2[0][0]), num[i] * num[i], MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Bcast(&(prod[0][0]), num[i] * num[i], MPI_DOUBLE, 0, MPI_COMM_WORLD);

            double startTime = MPI_Wtime();

            strassen(num[i], mat1, mat2, prod, p_rank);

            double endTime = MPI_Wtime();

            if (p_rank == 0)
            {
                // cout << "\nParallel Strassen Runtime (MPI): ";
                cout << (endTime - startTime)*1000 << endl;
                ff << (endTime - startTime)*1000 << endl;
                // cout << endl;
                // print(n, mat1);
                // print(n, mat2);
                // print(n, prod);
            }
            // }
            freeMatrix(num[i], mat1);
            freeMatrix(num[i], mat2);
            freeMatrix(num[i], prod);
        } else {
            int n = num[i] + 1;
            MPI_Barrier(MPI_COMM_WORLD);
            MPI_Bcast(&n, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

            double** mat1 = allocateMatrix(n-1);
            double** mat2 = allocateMatrix(n-1);
            double** mat11 = allocateMatrix(n);
            double** mat22 = allocateMatrix(n);
            double** prod = allocateMatrix(n);
            // cout << p_rank  << endl;

            if (p_rank == 0)
            {
                fillMatrix(n-1, mat1, 1001);
                mat11 = moreZeros(n-1, mat1);
                fillMatrix(n-1, mat2, 1017);
                mat22 = moreZeros(n-1, mat2);
            }

            freeMatrix(n-1, mat1);
            freeMatrix(n-1, mat2);

            MPI_Bcast(&(mat11[0][0]), n * n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Bcast(&(mat22[0][0]), n * n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            MPI_Bcast(&(prod[0][0]), n * n, MPI_DOUBLE, 0, MPI_COMM_WORLD);

            double startTime = MPI_Wtime();

            strassen(n, mat11, mat22, prod, p_rank);

            double endTime = MPI_Wtime();

            if (p_rank == 0)
            {
                // cout << "\nParallel Strassen Runtime (MPI): ";
                cout << (endTime - startTime)*1000 << endl;
                ff << (endTime - startTime)*1000 << endl;
                // cout << endl;
                // print(n, mat11);
                // print(n, mat22);
                // print(n, prod);
            }

            freeMatrix(n, mat11);
            freeMatrix(n, mat22);
            freeMatrix(n, prod);
        }
    }

    ff.close();
    MPI_Finalize();
    return 0;
}