CODEWITHSUNDEEP
pythonfortranpetsc
Artur
Artur

Reputation: 447

Reading PETSc binary matrix in Python

I am trying to read a sparse PETSc matrix saved from a Fortran code into a binary file like so:

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

When I run the Fortran code single-threaded, all works fine and I can read the matrix without any issues using the following Python code:

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()

However, when I run the Fortran code on more processors ("mpirun -n 2 my_exe"), I get the following error on the Python side:

    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

How can I fix this?

Upvotes: 0

Views: 133

Answers (1)

Artur
Artur

Reputation: 447

It seems that the issue was on the export side, not on the PETSc side. I originally wrote the export routine for symmetric matrices but then applied it to one that's not symmetric, leading to confused addressing in the binary file. I hadn't caught it earlier because I was using my own ASCII import routine on the Python side that didn't care about symmetry. The snippet in the question should work fine. I am now using the following two functions that I think are a bit cleaner:

import numpy as np
from scipy.sparse import coo_matrix, csr_matrix 
from petsc4py import PETSc

def readMat(filename, toScipy=False):
    viewer = PETSc.Viewer().createBinary(filename, mode='r', comm=PETSc.COMM_WORLD)
    A_petsc = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
    A_petsc.setType(PETSc.Mat.Type.AIJ)
    A_petsc.setFromOptions()
    A_petsc.load(viewer)
    if toScipy:
        I, J, V = A_petsc.getValuesCSR()
        return csr_matrix((V, J, I)).tocoo()
    else:
        return A_petsc

def readVec(filename, toNumpy=False):
    viewer = PETSc.Viewer().createBinary(filename, mode='r', comm=PETSc.COMM_WORLD)
    v_petsc = PETSc.Vec().create(comm=PETSc.COMM_WORLD)
    v_petsc.setFromOptions()
    v_petsc.load(viewer)
    if toNumpy:
        return np.array(v_petsc)
    else:
        return v_petsc

Upvotes: 0

Related Questions