MPI多进程分配求解矩阵行列式的实现问题咨询
现有一段采用Cramer's rule求解任意n×n矩阵行列式的C语言函数代码,已完成初步开发,经多组测试用例验证计算结果准确无报错。现计划对该代码进行MPI并行化改造,调度多个进程分别承担n阶矩阵的分段计算任务。
举例来说,若给定5 x 5矩阵、进程总数size=5,希望分配0~4号进程分别计算Cramer分解得到的前5个4 x 4子矩阵的行列式。目前正在自学MPI,该思路实现存在较多难点,因此求助。
行列式计算核心代码如下:
int determinate_solver(int r, int* ptr, int rank, int size) { int ans = 0, a = 0, b = 0, c = 0, d = 0; int inner_sol = 0, inner_det = 0; // 尝试根据进程数划分计算范围 int upperLimit = r; int start_val = rank * ceil(upperLimit / size) + 1, end_val; /* 不确定start_val和end_val的正确实现方式 */ if (rank == (size - 1)) { end_val = upperLimit; } else { end_val = start_val + ceil(upperLimit / size) - 1; } // 1阶、2阶矩阵直接计算 if (r == 1 || r == 2) { if (r == 1) { ans = ptr[0]; } else { a = ptr[0]; b = ptr[1]; c = ptr[2]; d = ptr[3]; ans = (a * d) - (b * c); } } else { int i, j, k, l, n = 0, sign = +1, basic, element; // 存储余子式的数组 int* q = (int*)calloc(((r - 1) * (r - 1)), sizeof(int)); for (int i = 0; i < r; i++) { l = 0; n = 0; basic = *(ptr + i); for (int j = 0; j < r; j++) { for (k = 0; k < r; k++) { element = *(ptr + l); if ((j == 0) || (i == k)); else { *(q + n) = element; n = n + 1; } l = l + 1; } } inner_det = determinate_solver (r - 1, q, rank, size); inner_sol = sign * basic * inner_det; ans = ans + inner_sol; sign = sign * -1; printf("The final solution is: %d from rank %d \n \n", ans, rank); } free(q); // 原代码漏了内存释放,会泄漏 } return ans; }
整个程序总代码量不足150行,为精简篇幅,本问题中省略了无关代码。
已尝试的方案
在上述代码中加入了end_val和start_val参数,计划通过这两个变量划定每个进程的有效计算范围。最初将第一层for循环修改为如下形式:
for(int i = start_val; i < end_val; i++){ /* 计算逻辑 */ }
但该方案运行失败,ans返回值错误。start_val和end_val在一维数组范围划定场景下表现正常,但在二维计算场景下存在适配问题。此外还尝试将代码中所有i实例替换为start_val,仍未生效。目前判断start_val和end_val可有效实现任务划分需求,但不确定正确的实现方式。
期望的答复方向
- 可提供可行的实现思路、技术说明
- 不使用Cramer's rule的替代实现方案也可接受
代码并行失败是三个核心逻辑错误导致的,按以下步骤修改即可正常运行:
拆分并行逻辑和串行递归逻辑
目前把MPI进程划分的逻辑写进了递归调用的行列式函数里,会导致每一层递归(比如计算4阶、3阶、2阶子矩阵时)都会重新给所有进程拆分计算任务,最后每个进程拿到的是多层递归里零散的计算片段,根本无法拼接成正确结果。
正确做法是单独写一个纯串行的行列式计算函数,不传入rank/size参数,只负责输入矩阵和阶数返回行列式值,所有递归计算都走这个串行函数,仅在最顶层的n阶矩阵计算时做MPI任务拆分。修正任务范围的索引和计算逻辑
原范围计算有两个低级错误:- C语言循环是0基索引,写的
start_val = rank * ceil(...) +1多了+1,会直接漏掉第一个展开项、甚至导致索引越界 ceil(upperLimit/size)完全不会生效:upperLimit和size都是int类型时,/是整数除法,会先向下取整再传给ceil,等于没算。
正确的块大小和范围计算用整数运算即可,不需要引入浮点函数:
int block_size = (n + size - 1) / size; // 等价于ceil(n/size)的整数实现 int start = rank * block_size; int end = (rank+1)*block_size > n ? n : (rank+1)*block_size;循环直接用
for(int i = start; i < end; i++)遍历当前进程负责的展开项即可。- C语言循环是0基索引,写的
补充结果汇总逻辑
每个进程只计算了自己负责的若干个代数余子式的和,必须通过MPI的规约操作把所有进程的局部结果累加,才能得到最终的行列式值。所有进程计算完局部和后,调用MPI_Reduce(&local_ans, &global_ans, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD),即可把所有进程的局部结果求和汇总到0号进程的global_ans变量中。
额外优化提示
- 原代码中
calloc申请的余子式数组q在循环结束后没有free,递归深度大时会出现严重内存泄漏,记得用完释放。 - 原代码余子式生成逻辑里
if ((j == 0) || (i == k));末尾的分号很容易引发误操作,建议改成显式的判断分支,比如if(j!=0 && i!=k) { 赋值逻辑 },可读性更强。 - 如果矩阵阶数较大,克莱姆法则的时间复杂度是O(n!),性能会极差,可替换为LU分解法计算行列式,时间复杂度可降到O(n³),并行改造的逻辑是通用的,只需要把串行行列式函数换成LU分解的实现即可。
内容的提问来源于stack exchange,提问作者Ammons_k

