基于Windows MPI的Cannon矩阵乘法实现遭遇死锁问题排查

环境与问题现象 #

• 环境：C++14、Windows MPI、Visual Studio 2022
• 运行命令：mpiexec -np 16 HW1.exe
• 实现逻辑：基于4x4笛卡尔拓扑的Cannon矩阵乘法，初始对齐阶段尝试按行坐标偏移量移位矩阵A
• 问题：仅首行进程能执行MPI_Cart_shift(comm_2d, 0, -mycoords[0], &shiftsource, &shiftdest);及后续MPI_Sendrecv_replace，其余进程陷入死锁；将偏移量改为固定-1时，所有进程正常执行。

死锁原因分析 #

死锁核心是通信配对不匹配：
当使用-mycoords[0]作为移位偏移量时，不同行进程的移位目标/来源无法形成双向配对的通信关系：

• 行0进程（mycoords[0]=0）：偏移量为0，MPI_Cart_shift返回的发送/接收目标都是自身，可正常完成通信
• 行1进程（mycoords[0]=1）：偏移量为-1，需向左移1步，但行0进程的通信对象是自己，不会接收行1的消息
• 行2/3进程（mycoords[0]=2/3）：偏移量为-2/-3，发送目标均为行0进程，但行0进程不会处理这些外来消息，导致行2/3的MPI_Sendrecv_replace因无人接收而阻塞，最终引发死锁。

而固定偏移-1时，所有进程均执行左移1步操作，每个进程的发送目标是左侧进程，接收来源是右侧进程，形成闭环通信，无配对缺失，因此不会死锁。

修正方案 #

Cannon算法的初始对齐要求第i行进程左移i步、第j列进程上移j步，正确的实现方式是分多次执行单步移位，而非一次性移位i步，确保每次通信都有配对对象：

1. 修正矩阵A的初始对齐代码 #

将原来的单次移位替换为循环执行mycoords[0]次左移1步：

// 原代码：
// MPI_Cart_shift(comm_2d, 0, -mycoords[0], &shiftsource, &shiftdest);
// MPI_Sendrecv_replace(a, nlocal * nlocal, MPI_DOUBLE, shiftdest, 1, shiftsource, 1, comm_2d, &status);

// 修正后：
for (int s = 0; s < mycoords[0]; ++s) {
    MPI_Cart_shift(comm_2d, 0, -1, &shiftsource, &shiftdest);
    MPI_Sendrecv_replace(a, nlocal * nlocal, MPI_DOUBLE, shiftdest, 1, shiftsource, 1, comm_2d, &status);
}

2. 修正矩阵B的初始对齐代码 #

同理，将矩阵B的单次移位替换为循环执行mycoords[1]次上移1步：

// 原代码：
// MPI_Cart_shift(comm_2d, 1, -mycoords[1], &shiftsource, &shiftdest);
// MPI_Sendrecv_replace(b, nlocal * nlocal, MPI_DOUBLE, shiftdest, 1, shiftsource, 1, comm_2d, &status);

// 修正后：
for (int s = 0; s < mycoords[1]; ++s) {
    MPI_Cart_shift(comm_2d, 1, -1, &shiftsource, &shiftdest);
    MPI_Sendrecv_replace(b, nlocal * nlocal, MPI_DOUBLE, shiftdest, 1, shiftsource, 1, comm_2d, &status);
}

3. 额外代码问题修正 #

• 转置逻辑错误：主函数中的双重循环会将矩阵转置两次，等于未转置，应改为仅交换上三角元素：

// 原转置代码：
// for (int i = 0; i < 4; ++i) {
//     for (int j = 0; j < 4; ++j) {
//         ...交换操作...
//     }
// }

// 修正后：
for (int i = 0; i < 4; ++i) {
    for (int j = i + 1; j < 4; ++j) {
        double t = a[i * 4 + j];
        a[i*4+j] = a[j*4+i];
        a[j * 4 + i] = t;

        t = b[i * 4 + j];
        b[i*4+j] = b[j*4+i];
        b[j * 4 + i] = t;
    }
}

• 整数转换安全：std::sqrt(npes)返回double类型，赋值给int时需显式转换：

dims[0] = dims[1] = static_cast<int>(std::sqrt(npes));

完整修正后代码 #

#include <mpi.h>
#include <stdio.h>
#include <cmath>

void MatrixMatrixMultiply(int n, double* a, double* b, double* c, MPI_Comm comm)
{
    int i;
    int nlocal;
    int npes, dims[2], periods[2];
    int myrank, my2drank, mycoords[2];
    int uprank, downrank, leftrank, rightrank;
    int shiftsource, shiftdest;
    MPI_Status status;
    MPI_Comm comm_2d;

    /* Get the communicator related information */
    MPI_Comm_size(comm, &npes);
    MPI_Comm_rank(comm, &myrank);

    /* Set up the Cartesian topology */
    dims[0] = dims[1] = static_cast<int>(std::sqrt(npes));

    /* Set the periods for wraparound connections */
    periods[0] = periods[1] = 1;

    /* Create the Cartesian topology, with rank reordering */
    MPI_Cart_create(comm, 2, dims, periods, 1, &comm_2d);

    /* Get the rank and coordinates with respect to the new topology */
    MPI_Comm_rank(comm_2d, &my2drank);
    MPI_Cart_coords(comm_2d, my2drank, 2, mycoords);

    /* Compute ranks of the up and left shifts */
    MPI_Cart_shift(comm_2d, 0, -1, &rightrank, &leftrank);
    MPI_Cart_shift(comm_2d, 1, -1, &downrank, &uprank);

    /* Determine the dimension of the local matrix block */
    nlocal = n / dims[0];

    /* Perform the initial matrix alignment. First for A */
    for (int s = 0; s < mycoords[0]; ++s) {
        MPI_Cart_shift(comm_2d, 0, -1, &shiftsource, &shiftdest);
        MPI_Sendrecv_replace(a, nlocal * nlocal, MPI_DOUBLE, shiftdest, 1, shiftsource, 1, comm_2d, &status);
    }

    fprintf(stdout, "rank: %d\n", myrank);
    fflush(stdout);

    /* Then for B */
    for (int s = 0; s < mycoords[1]; ++s) {
        MPI_Cart_shift(comm_2d, 1, -1, &shiftsource, &shiftdest);
        MPI_Sendrecv_replace(b, nlocal * nlocal, MPI_DOUBLE, shiftdest, 1, shiftsource, 1, comm_2d, &status);
    }

    /* Get into the main computation loop */
    for (i = 0; i < dims[0]; i++)
    {
        // 模拟本地矩阵乘法（实际需替换为正确的块乘逻辑）
        c[0] += a[0] * b[0];

        /* Shift matrix a left by one */
        MPI_Sendrecv_replace(a, nlocal * nlocal, MPI_DOUBLE, leftrank, 1, rightrank, 1, comm_2d, &status);

        /* Shift matrix b up by one */
        MPI_Sendrecv_replace(b, nlocal * nlocal, MPI_DOUBLE, uprank, 1, downrank, 1, comm_2d, &status);
    }

    /* Restore the original distribution of a and b */
    for (int s = 0; s < mycoords[0]; ++s) {
        MPI_Cart_shift(comm_2d, 0, 1, &shiftsource, &shiftdest);
        MPI_Sendrecv_replace(a, nlocal * nlocal, MPI_DOUBLE, shiftdest, 1, shiftsource, 1, comm_2d, &status);
    }

    for (int s = 0; s < mycoords[1]; ++s) {
        MPI_Cart_shift(comm_2d, 1, 1, &shiftsource, &shiftdest);
        MPI_Sendrecv_replace(b, nlocal * nlocal, MPI_DOUBLE, shiftdest, 1, shiftsource, 1, comm_2d, &status);
    }

    MPI_Comm_free(&comm_2d); /* Free up communicator */
}

int main(int argc, char** argv)
{
    MPI_Init(NULL, NULL);

    double a[16] = 
    {
        2.0, 3.0, 4.0, 5.0,
        9.0, 8.0, 7.0, 6.0,
        5.0, 4.0, 2.0, 3.0,
        8.0, 7.0, 3.0, 4.0
    };
    double b[16] = {
        3.0, 5.0, 7.0, 6.0,
        2.0, 7.0, 6.0, 3.0,
        7.0, 5.0, 3.0, 2.0,
        4.0, 3.0, 2.0, 5.0
    };
    double c[16] = {0.0};
    
    // 正确转置矩阵（仅交换上三角元素）
    for (int i = 0; i < 4; ++i) {
        for (int j = i + 1; j < 4; ++j) {
            double t = a[i * 4 + j];
            a[i*4+j] = a[j*4+i];
            a[j * 4 + i] = t;

            t = b[i * 4 + j];
            b[i*4+j] = b[j*4+i];
            b[j * 4 + i] = t;
        }
    }

    MatrixMatrixMultiply(4, a, b, c, MPI_COMM_WORLD);

    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    if (rank == 0)
    {
        for (int i = 0; i < 16; ++i)
        {
            printf("%f ", c[i]);
            if ((i + 1) % 4 == 0)
                printf("\n");
        }
    }

    MPI_Finalize();

    return 0;
}

内容的提问来源于stack exchange，提问作者Stuart