C语言下ScaLAPACK矩阵列优先处理的阶段疑问

核心前提 #

ScaLAPACK底层基于Fortran实现，因此无论你用C还是Fortran编程，进程本地持有的矩阵块集合必须以列优先方式连续存储，这是调用ScaLAPACK核心函数的硬性要求。结合你的输入-计算-输出流程，各阶段的处理规则如下：

输入阶段 #

• 文件存储格式不强制列优先，但读取后必须确保进程本地的矩阵块是列优先布局：
  • 若输入文件是行优先（如C语言常规存储：1 2 3 4 5 6 7 8 9），需通过MPI-IO的MPI_Type_create_subarray创建子数组类型，跳过行优先的存储间隔，直接将列元素映射到本地列优先的数组中；或读取后在本地完成行优先到列优先的转换。
  • 若输入文件本身是列优先（如1 4 7 2 5 8 3 6 9），可直接按连续块读取到本地，无需额外转换，效率更高。
• 进程间的分布：ScaLAPACK采用二维块循环分布（依赖BLACS进程网格），每个进程持有的块在本地内存中必须列优先连续。比如进程0持有原矩阵第一列的1,4,7，在C数组中需存储为[1,4,7]（列优先下同一列的元素连续）。

计算阶段 #

• 只要调用ScaLAPACK的核心运算函数（如矩阵乘法pdgemm），必须保证本地矩阵块是列优先存储，这是函数的输入约定，和你用C还是Fortran无关。
• 若你自行实现矩阵运算（不依赖ScaLAPACK），可选择行优先存储，但如果后续要和ScaLAPACK交互，仍需转换为列优先布局。另外，i-k-j循环优化是针对行优先存储的缓存友好性优化，和ScaLAPACK的列优先存储要求不冲突——你可以在保持本地列优先布局的前提下，调整循环顺序来优化计算性能。

输出阶段 #

• 与输入阶段对称，需根据目标文件的使用场景决定：
  • 若输出给ScaLAPACK/Fortran程序使用，可直接将本地列优先的块按分布规则写入文件，无需转换。
  • 若输出给C程序使用，可先在本地将列优先块转换为行优先，再通过MPI-IO写入；或直接用MPI_Type_create_subarray创建行优先的文件视图，将本地列优先块直接映射写入到文件的对应位置，避免本地转换开销。

总结 #

整个流程的核心约束点在计算阶段：只要涉及ScaLAPACK的调用，进程本地的矩阵块必须是列优先存储。输入和输出阶段则可根据文件的存储格式、目标使用场景灵活处理，核心是确保进入计算阶段前数据符合ScaLAPACK的布局要求，输出时按需转换格式即可。

内容的提问来源于stack exchange，提问作者G. Ianni