使用OpenMP做1000x1000矩阵乘法无输出，10x10可正常运行是什么原因

问题根因及解决方案 #

1. 核心问题：栈内存溢出 #

你将三个1000x1000的double类型数组声明为局部变量，而局部变量默认分配在栈空间：

• 单个double占8字节，单个1000*1000的double数组占8MB，三个数组总占用24MB
• 主流操作系统默认栈空间远小于该阈值：Linux默认栈大小为8MB，Windows默认仅1MB，大数组直接撑爆栈内存，会导致程序直接崩溃无输出

10x10规模下数组总占用仅2.4KB，远小于栈容量，因此可以正常运行

解决方法：

• 把三个数组声明移到函数外作为全局变量（全局变量分配在堆空间，不受栈大小限制）
• 或者使用malloc动态申请数组内存

2. 次要问题：运行效率极低，易被误认为无输出 #

你当前使用的i-j-k三重循环顺序，对矩阵b的访问是按列跳跃访问，cache命中率极低，1000阶矩阵乘法在优化不足的情况下可能运行数分钟都出不来结果，容易被误认为程序无响应。
解决方法：

• 调整循环顺序为i-k-j，让两个矩阵的访问都符合行优先的存储顺序，cache命中率可提升数十倍：

#pragma omp for
for (i=0; i<n; i++) {
    for (k=0; k<n; k++){
        double tmp = a[i][k];
        for (j=0; j<n; j++) {
            mul[i][j] += tmp * b[k][j]; 
        } 
    }
}

• 取消矩阵乘法循环内的printf打印，逐行打印会极大拖慢运行速度

3. 隐藏问题：线程安全风险 #

rand()不是线程安全函数，多个线程同时调用会出现数据竞争，不仅生成的随机数不符合预期，还可能引发未知异常。
解决方法：

• 改用线程安全的随机数生成函数，或者在并行区外单线程生成两个矩阵的随机值

4. 其他小问题：变量未初始化 #

norm变量没有初始化为0，最终计算的Frobenius范数结果会出错。

内容的提问来源于stack exchange，提问作者P.Chian