Netlib ScaLAPACK example1.f执行结果不稳定问题求助

1. 验证进程网格与MPI映射的一致性 #

• 确保代码中进程数校验逻辑生效：原example1包含IF( NPROCS.NE.NPROW*NPCOL ) THEN判断，修改NPROW/NPCOL后需确认该分支能正确捕获进程数不匹配的情况（尽管你指定了-n 4，但可避免意外场景）。
• 打印进程网格坐标：在BLACS_GRIDINIT后添加代码，输出每个进程的MYROW和MYCOL，检查每次运行的坐标分配是否一致，排除网格映射异常：

  CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
  WRITE(*,*) 'Process ', MYPID, ': MYROW=', MYROW, ', MYCOL=', MYCOL
  

• 尝试调整网格映射模式：将BLACS_GRIDINIT的第二个参数从'Row'改为'Column'，测试结果是否稳定。

2. 调整编译链接顺序与方式 #

• 修正库链接顺序：链接器按从左到右解析符号，正确顺序需保证依赖库在被依赖库之后，建议改为：

  mpif77 example1.f -o example1 libscalapack.a liblapack.a librefblas.a libmpi.a
  

• 改用编译器自带的库链接选项：避免手动指定静态库，让mpif77自动处理MPI依赖：

  mpif77 example1.f -o example1 -lscalapack -llapack -lrefblas
  

• 校验库编译一致性：确保ScaLAPACK、LAPACK、BLAS均使用与mpif77相同的编译器、浮点选项（如64位模式、精度设置）编译，避免ABI不兼容。

3. 统一浮点环境与编译优化选项 #

• 显式指定浮点编译选项：添加选项保证所有编译单元浮点行为一致，例如GCC下：

  mpif77 -O2 -fno-trapping-math -ffast-math example1.f -o example1 -lscalapack -llapack -lrefblas
  

• 统一优化级别：确保示例程序与所有依赖库使用相同的优化级别（如均用-O2），避免因优化差异导致浮点计算行为不一致。
• 抑制无关浮点警告：IEEE非规格化数和下溢通常不影响计算正确性，可通过编译选项抑制相关输出，避免干扰结果判断。

4. 检查数据初始化与残差计算逻辑 #

• 验证局部数据大小：打印每个进程的局部矩阵维度LOCNR、LOCNC和LDA，确认修改进程网格后，局部数据分配和初始化正确，无未初始化内存：

  WRITE(*,*) 'Process ', MYPID, ': LOCNR=', LOCNR, ', LOCNC=', LOCNC, ', LDA=', LDA
  

• 检查PDGESV返回状态：调用PDGESV后立即打印INFO参数，错误运行时若INFO≠0，可根据ScaLAPACK文档定位问题（如矩阵奇异、数据分布错误）：

  CALL PDGESV( N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB, INFO )
  WRITE(*,*) 'Process ', MYPID, ': PDGESV INFO=', INFO
  

• 排查残差计算的MPI通信：原示例残差通过MPI_REDUCE全局求和得到，检查通信的根进程、操作类型（MPI_SUM）是否正确，避免通信错误引入异常值。

5. 对比PDGESV调用参数与测试用例xdsvd #

• 调整矩阵分块大小：原示例MB=16、NB=16，对于N=100的矩阵会出现非整数块，尝试改为MB=10、NB=10（整除100），测试结果是否稳定。
• 对齐xdsvd的调用参数：对比xdsvd中PDGESV的输入参数（如DESCA设置、IPIV数组大小、矩阵存储格式），确保example1的参数设置与测试用例一致。

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