多进程Fortran保存PETSc二进制矩阵，Python读取报错如何解决？

问题背景 #

Fortran多进程保存稀疏PETSc矩阵后，Python读取时出现数据不一致错误；但单线程保存矩阵时，Python可正常读取。

Fortran保存代码 #

CALL PetscViewerBinaryOpen(PETSC_COMM_SELF, TRIM(ADJUSTL(filename)), FILE_MODE_WRITE, viewer2, ier)
CALL MatView(aa_symmetric, viewer2, ier)
CALL PetscViewerDestroy(viewer2, ier)

单线程正常读取的Python代码 #

import sys
import petsc4py
from petsc4py import PETSc

# Initialize PETSc
petsc4py.init(sys.argv)

# Create a viewer for reading the binary file
viewer = PETSc.Viewer().createBinary(filename, mode='r', comm=PETSc.COMM_WORLD)

# Create a matrix and load data from the binary file
A_petsc = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
A_petsc.setType(PETSc.Mat.Type.AIJ)
A_petsc.setFromOptions()
A_petsc.load(viewer)

# Finalize PETSc
PETSc.Finalize()

多进程保存后的错误信息 #

Traceback (most recent call last): 
File "/home/Python/test_matrixImport_binary.py", line 80, in <module>
        A_petsc.load(viewer)   File "petsc4py/PETSc/Mat.pyx", line 2025, in 
petsc4py.PETSc.Mat.load petsc4py.PETSc.Error: error code 79 [0] MatLoad() at 
/home/lib/petsc-3.21.0/src/mat/interface/matrix.c:1344 [0] MatLoad_SeqAIJ() at 
/home/lib/petsc-3.21.0/src/mat/impls/aij/seq/aij.c:5091 [0] MatLoad_SeqAIJ_Binary() at 
/home/lib/petsc-3.21.0/src/mat/impls/aij/seq/aij.c:5142 [0] Unexpected data in file [0] 
Inconsistent matrix data in file: nonzeros = 460, sum-row-lengths = 761

解决方案 #

问题根源是Fortran保存时使用PETSC_COMM_SELF，多进程下每个进程会独立将自身的矩阵分片写入同一文件，导致文件内容被覆盖、混杂，损坏数据结构。解决步骤如下：

1. 修改Fortran保存代码，使用全局通信子 #

将PETSC_COMM_SELF替换为全局通信器PETSC_COMM_WORLD，让PETSc统一协调多进程的矩阵写入，生成符合并行格式的二进制文件：

CALL PetscViewerBinaryOpen(PETSC_COMM_WORLD, TRIM(ADJUSTL(filename)), FILE_MODE_WRITE, viewer2, ier)
CALL MatView(aa_symmetric, viewer2, ier)
CALL PetscViewerDestroy(viewer2, ier)

2. 适配Python读取逻辑 #

• 单进程读取并行保存的矩阵：创建viewer时指定comm=PETSc.COMM_SELF，确保进程读取完整的并行矩阵数据：

  import sys
  import petsc4py
  from petsc4py import PETSc
  
  petsc4py.init(sys.argv)
  
  # 单进程读取并行格式的二进制文件
  viewer = PETSc.Viewer().createBinary(filename, mode='r', comm=PETSc.COMM_SELF)
  A_petsc = PETSc.Mat().create(comm=PETSc.COMM_SELF)
  A_petsc.setType(PETSc.Mat.Type.AIJ)
  A_petsc.setFromOptions()
  A_petsc.load(viewer)
  
  PETSc.Finalize()
  

• 多进程读取并行保存的矩阵：保持comm=PETSc.COMM_WORLD，用mpirun -n 2 python your_script.py启动脚本，PETSc会自动分配各进程读取对应的矩阵分片。

关键原理 #

PETSc的二进制格式在并行场景下，需要通过全局通信器协调各进程的写入操作，确保文件包含完整的并行矩阵元数据与分片信息。使用PETSC_COMM_SELF会让每个进程独立写入，破坏文件结构，导致读取时数据校验失败。

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