Scalapack pdgels最小二乘多MPI任务运行异常求助

问题现象 #

使用Scalapack的pdgels最小二乘例程时，单MPI任务运行正常，B_loc能正确更新为最小二乘解；但4个MPI任务（2x2 BLACS网格）运行时，仅rank 0的B_loc得到更新，其余进程的B_loc初始为无效垃圾值且未被更新。预期是将8x8全局矩阵A和向量B分发到4个任务，每个任务对应4x4本地矩阵A_loc及对应B_loc。

运行输出 #

单任务运行输出 #

BEFORE: Local B size 8 on rank: 0 values: 1 2 3 4 5 6 7 8 
AFTER: Local B size 8 on rank: 0 values: 0.960918 0.0252256 0.0222086 0.0319953 -0.0471391 0.0205009 -0.0163878 0.00267855 

四任务运行输出 #

BEFORE: Local B size 4 on rank: 0 values: 1 2 3 4 
BEFORE: Local B size 4 on rank: 1 values: 4.68603e-310 0 2.122e-314 2.122e-314 
BEFORE: Local B size 4 on rank: 2 values: 5 6 7 8 
BEFORE: Local B size 4 on rank: 3 values: 4.67566e-310 0 2.122e-314 2.122e-314 
AFTER: Local B size 4 on rank: 1 values: 4.68603e-310 0 2.122e-314 2.122e-314 
AFTER: Local B size 4 on rank: 2 values: 5 6 7 8 
AFTER: Local B size 4 on rank: 3 values: 4.67566e-310 0 2.122e-314 2.122e-314 
AFTER: Local B size 4 on rank: 0 values: 0.321957 1.12458 -0.21503 -0.231507 

问题分析与修正 #

1. pdgels参数传递错误 #

pdgels要求前两个参数M和N是全局矩阵的行数和列数，而非本地矩阵的大小。原代码错误传递m_loc和n_loc，导致多进程下无法正确识别全局问题规模。

修正：将pdgels调用中的&m_loc和&n_loc替换为全局的&m和&n：

pdgels_(&trans, &m, &n, &nrhs, A_loc, &ia, &ja, descA_loc,
    B_loc, &ib, &jb, descB_loc, &wkopt, &lwork, &info);

// 第二次调用同样修正
pdgels_(&trans, &m, &n, &nrhs, A_loc, &ia, &ja, descA_loc,
    B_loc, &ib, &jb, descB_loc, work, &lwork, &info);

2. 全局矩阵初始化与分发逻辑错误 #

原代码中所有进程都初始化了全局矩阵A和B，但根据BLACS块描述（mb=m、nb=n），只有rank 0（BLACS网格的(0,0)进程）拥有全局矩阵的有效块，其他进程的全局矩阵是未初始化的垃圾值，导致pdgemr2d_分发时读取错误数据。

修正：仅在rank 0进程初始化全局矩阵，其他进程无需分配或初始化全局矩阵：

// 全局矩阵仅在rank0分配初始化
double *A = nullptr, *B = nullptr;
if (rank == 0) {
    A = new double[m*n];
    B = new double[m];
    // 生成列优先矩阵
    for (int i = 0; i < m*n; i++) {
        A[i] = i+1;
    }
    // 生成目标向量
    for (int i = 0; i < m; i++) {
        B[i] = (i+1);
    }
}

3. BLACS与MPI变量冲突 #

原代码调用blacs_pinfo_(&myid, &size)时，覆盖了之前通过MPI_Comm_size获取的size变量，可能导致BLACS网格初始化错误。

修正：使用独立变量存储BLACS进程总数：

int blacs_myid, blacs_size;
blacs_pinfo_( &blacs_myid, &blacs_size );

4. 内存释放语句错误 #

原代码使用delete [] work, A , B, A_loc, B_loc;，逗号运算符只会释放work，其他内存无法正确释放，造成内存泄漏。

修正：分别释放每个动态分配的数组：

delete [] work;
if (rank == 0) {
    delete [] A;
    delete [] B;
}
delete [] A_loc;
delete [] B_loc;

5. 描述子初始化验证 #

确保descinit_调用后检查info值，确认描述子初始化成功：

descinit_( descA,  &m, &n, &mb, &nb, &izero, &izero, &context, &m, &info);
if (info != 0) {
    printf("Rank %d: descinit for descA failed with info %d\n", rank, info);
    MPI_Abort(MPI_COMM_WORLD, 1);
}
// 其他descinit调用同理添加检查

修正后关键逻辑说明 #

1. 仅rank 0初始化全局矩阵，通过pdgemr2d_将数据分发到各个进程的本地矩阵A_loc和B_loc。
2. pdgels使用全局问题规模参数，确保多进程下协同计算最小二乘解。
3. 修复内存管理问题，避免内存泄漏和未初始化数据访问。

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