高斯约旦法求矩阵逆的循环展开优化错误排查

你的循环展开版本逻辑完全偏离了高斯约旦法的核心步骤，导致计算结果错误，核心问题如下：

• 重复归一化：在i循环的每次迭代中都重新对k+1行执行归一化，违背了每行仅需归一化一次的原则，导致k+1行被多次缩放，数值完全失真。
• 主元选择错误：归一化k+1行时，你误用了original[k][k]作为主元，正确的主元应该是original[k+1][k+1]（k+1行的对角线元素）。
• 消元顺序混乱：高斯约旦法要求先完成第k行的所有消元操作（确保所有行的k列仅k行是1，其余为0），再处理k+1行。你的代码将两行的处理混在i循环中，消元操作互相干扰，破坏了矩阵的中间状态。
• 边界未处理：当矩阵大小为奇数时，最后一行会被完全忽略，因为k循环步长为2，无法覆盖到最后一个索引。

修正后的循环展开实现（保留缓存优化） #

以下代码严格遵循高斯约旦法步骤，同时保持行优先的内存访问顺序（缓存友好），实现两步循环展开：

// 处理偶数部分，每次批量处理两行
for(k=0; k < sizeOfMatrix - 1; k += 2) {
    // -------------------------- 处理第k行 --------------------------
    // 1. 归一化k行（主元为k行对角线元素）
    double pivot_k = original[k][k];
    for(j=0; j < sizeOfMatrix; j++) {
        original[k][j] /= pivot_k;
        inverse[k][j] /= pivot_k;
    }
    // 2. 用k行消去所有其他行的k列
    for(i=0; i < sizeOfMatrix; i++) {
        if(i != k) {
            double pivot = original[i][k];
            for(j=0; j < sizeOfMatrix; j++) {
                original[i][j] -= original[k][j] * pivot;
                inverse[i][j] -= inverse[k][j] * pivot;
            }
        }
    }

    // -------------------------- 处理第k+1行 --------------------------
    // 1. 归一化k+1行（此时k+1行的k列已被消为0，主元为自身对角线元素）
    double pivot_k1 = original[k+1][k+1];
    for(j=0; j < sizeOfMatrix; j++) {
        original[k+1][j] /= pivot_k1;
        inverse[k+1][j] /= pivot_k1;
    }
    // 2. 用k+1行消去所有其他行的k+1列
    for(i=0; i < sizeOfMatrix; i++) {
        if(i != k+1) {
            double pivot = original[i][k+1];
            for(j=0; j < sizeOfMatrix; j++) {
                original[i][j] -= original[k+1][j] * pivot;
                inverse[i][j] -= inverse[k+1][j] * pivot;
            }
        }
    }
}

// 处理奇数大小矩阵的最后一行（避免遗漏）
if(sizeOfMatrix % 2 != 0) {
    k = sizeOfMatrix - 1;
    double pivot = original[k][k];
    for(j=0; j < sizeOfMatrix; j++) {
        original[k][j] /= pivot;
        inverse[k][j] /= pivot;
    }
    for(i=0; i < sizeOfMatrix; i++) {
        if(i != k) {
            double pivot_i = original[i][k];
            for(j=0; j < sizeOfMatrix; j++) {
                original[i][j] -= original[k][j] * pivot_i;
                inverse[i][j] -= inverse[k][j] * pivot_i;
            }
        }
    }
}

进一步缓存优化建议 #

如果想进一步提升缓存命中率，可以对j循环也进行展开（连续内存访问的循环展开收益更高），示例如下：

// j循环4倍展开（需处理剩余元素）
for(j=0; j < sizeOfMatrix - 3; j +=4) {
    original[k][j] /= pivot_k;
    original[k][j+1] /= pivot_k;
    original[k][j+2] /= pivot_k;
    original[k][j+3] /= pivot_k;
    
    inverse[k][j] /= pivot_k;
    inverse[k][j+1] /= pivot_k;
    inverse[k][j+2] /= pivot_k;
    inverse[k][j+3] /= pivot_k;
}
// 处理剩余的1-3个元素
for(; j < sizeOfMatrix; j++) {
    original[k][j] /= pivot_k;
    inverse[k][j] /= pivot_k;
}

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