mpi中的矩阵乘法和全局归约操作

Question

例如，如果我想将 4 的 2 乘以 2 矩阵相乘，为简单起见，假设它们都相同并且具有条目 1。

然后我想知道我应该如何使用全局减少 mpi 来并行化这个？假设大小为 4。

你能给我一个这样做的想法吗？谢谢！

# include <stdio.h>
# include <mpi.h>
# define N 4

//Create the 2 times 2 matrix type
typedef double Matrix[2][2];

void printMatrix(Matrix m);
void unitMatrix(Matrix m);
void randomMatrix(Matrix m);
void multMatrix(Matrix r, Matrix a, Matrix b);
void copyMatrix(Matrix out, Matrix in);
double random_number(void);
void my_range(int n, int *i1, int *i2);

int main(int argc, char *argv[])
{
    //Create a single matrix a
    Matrix a;
    Matrix buf;
    //Create a set of 100 matrix
    Matrix b[N];
    int i;
    int rank, i1, i2;
    double row1[2];
    double row2[2];
    double col1[2];
    double col2[2];


    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    my_range(N, &i1, &i2);

    for (i = 0; i < i1; i++) {
        randomMatrix(a);
    }
    for (i = i1; i <= i2; i++) {
    randomMatrix(b[i]);
    }
    for (i = i2 + 1; i < N; i++) {
        randomMatrix(a);
    }

    unitMatrix(a);

    for (i = i1; i <= i2; i++) {
    multMatrix(a, a, b[i]);
    MPI_Reduce(&a,&buf, 4, MPI_DOUBLE, MPI_PROD, 0,
             MPI_COMM_WORLD);
    }

    if (rank == 0) printMatrix(buf);

    MPI_Finalize();
    return 0;
}

//print a single matrix
void printMatrix(Matrix m) 
{
    printf("%26.18e %26.18e %26.18e %26.18e\n", 
       m[0][0], m[0][1], m[1][0], m[1][1]);
}

void unitMatrix(Matrix m)
{
    m[0][0] = 1.0;
    m[0][1] = 0.0;
    m[1][0] = 0.0;
    m[1][1] = 1.0;
}

void randomMatrix(Matrix m)
{
    m[0][0] = 1.0;
    m[0][1] = 1.0;
    m[1][0] = 1.0;
    m[1][1] = 1.0;
}

double random_number(void)
{
    const int mr = 714025;
    const int ia = 1366;
    const int ic = 150889;
    const double qdnorm = 1.0 / mr;
    static int irandom = 0;

    irandom = (ia * irandom + ic) % mr;
    return(irandom * qdnorm);
}

void multMatrix(Matrix r, Matrix a, Matrix b)
{
    // multMatrix(r, a, b) calculates  r = a * b
    // multMatrix(a, a, b) calculates  a = a * b
    // multMatrix(a, b, a) calculates  a = b * a

    Matrix tmp;

    tmp[0][0] = a[0][0] * b[0][0] + a[1][0] * b[0][1];
    tmp[0][1] = a[0][1] * b[0][0] + a[1][1] * b[0][1];
    tmp[1][0] = a[0][0] * b[1][0] + a[1][0] * b[1][1];
    tmp[1][1] = a[0][1] * b[1][0] + a[1][1] * b[1][1];

    copyMatrix(r, tmp);
}

void copyMatrix(Matrix out, Matrix in)
{
    out[0][0] = in[0][0];
    out[0][1] = in[0][1];
    out[1][0] = in[1][0];
    out[1][1] = in[1][1];
}

void my_range(int n, int *i1, int*i2)
{
    int size, rank, chunk, rest;

    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    chunk = n / size;
    rest = n % size;

    if (rank < rest) {
    chunk = chunk + 1;
    *i1 = chunk * rank;
    } else {
    *i1 = chunk * rank + rest;
    }

    *i2 = *i1 + chunk - 1;
}

Answer 1

您的代码使用逐元素矩阵乘法减少部分结果，即r[i][j] = a[i][j] * b[i][j]，因此给您错误的结果。正如 haraldkl 已经指出的那样，您可以将 MPI 的机制用于用户定义的 MPI 缩减运算符MPI_Op_create。您还应该创建用户定义的 MPI 数据类型，以便能够将每个数组作为单个矩阵实体处理。例如：

void myMatrixProd(Matrix *in, Matrix *inout, int *len, MPI_Datatype *dptr)
{
   int i;

   for (i = 0; i < *len; i++)
   {
      multMatrix(inout[i], in[i], inout[i]);
   }
}

...

MPI_Op multOp;
MPI_Datatype matrixType;

MPI_Type_contiguous(2*2, MPI_DOUBLE, &matrixType);
MPI_Type_commit(&matrixType);

MPI_Op_create(myMatrixProd, 0, &multOp);

Matrix a, buf;

// Compute partial product into a
multMatrix(...);

// Reduce the partial products to get the total into rank 0
MPI_Reduce(&a, &buf, 1, matrixType, multOp, 0, MPI_COMM_WORLD);

需要注意的重要一点是MPI_Op_create 的第二个参数是0。这是一个标志，指示归约运算符是否可交换。矩阵乘法不是可交换的（但仍然是所有 MPI 归约运算符所要求的关联），因此应该在此处指定 0。

Answer 2

据我了解您的问题，您正在寻找user defined MPI operation。