PETSc - Manipulate BAIJ matrix locally

Question

My program loads a parallel PETSc matrix $A$ on several MPI processes, each holding a block submatrix $A_i$.

I would like to retrieve the local submatrix $A_i$, the one corresponding to the current rank.

I am not trying to alter the matrix $A_i$, I only need to read it and compute some information (determinant, eigenvalues, etc).

int main(int argc, char **args) {
    PetscErrorCode  ierr;
    PetscViewer     fd;
    Mat A;

    ierr = SlepcInitialize(&argc, &args, (char *) nullptr, nullptr);
    if (ierr) return ierr;

    // Open file
    PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.mat", FILE_MODE_READ, &fd);

    // Load into Petsc object
    MatCreate(PETSC_COMM_WORLD, &A);
    MatSetType(A, MATMPIBAIJ);
    MatLoad(A, fd);

    PetscViewerDestroy(&fd);

    // Here I would like to manipulate the local submatrix of A
    // For example to show information with MatView


    // Clear workspace
    MatDestroy(&A);
    SlepcFinalize();
    return 0;
}

What would be the right way to retrieve the local submatrix for the current rank?

Note: I am using SlepcInitialize in order to use SLEPc for eigenvalues afterwards.

Edit :

I have pointers towards VecGetArray, I think I need an equivalent for matrices, but MatGetArray is not supposed to be used according to the documentation (and another post)

Answer 1

The solution I found for PETSc 3.7 was to use MatGetSubMatrices.

Warning : this method is not present in the last API.

I discovered that BAIJ is actually not intended for block matrices but rather to keep rows by group when distributing the matrix over processors.

Thus I used AIJ format. (MPIAIJ to get a parallel matrix).

First have the rights rows and columns indices (this is not PETSc related).

My methods split_indices_rows and split_indices_rows split the indices so that each MPI rank know which part of the matrix to handle by calling indices_rows_rank[mpi_rank] and indices_cols_rank[mpi_rank].

// indices_rows contain all rows indices of original matrix : {O , .. , height}
vector<int> indices_rows;
for (int i = 0; i < matrix_height; i++)
    indices_rows.push_back(i);

// indices_columns contain all columns indices of original matrix : {0 , .. , width}
vector<int> indices_columns;
for (int i = 0; i < matrix_width; i++)
    indices_columns.push_back(i);

vector<vector<int> > indices_rows_rank = split_indices_rows(indices_rows);
vector<vector<int> > indices_cols_rank = split_indices_cols(indices_cols);

Then you can create index sets (IS).

int mpi_rank;
MPI_Comm_rank(PETSC_COMM_WORLD, &mpi_rank);

// Create index sets
vector<int> local_rows = indices_rows_rank[mpi_rank];
IS local_row_indexes;
ISCreateGeneral(PETSC_COMM_WORLD, local_rows.size(), &local_rows[0], PETSC_COPY_VALUES, &local_row_indexes);

vector<int> local_cols = indices_cols_rank[mpi_rank];
IS local_col_indexes;
ISCreateGeneral(PETSC_COMM_WORLD, local_cols.size(), &local_cols[0], PETSC_COPY_VALUES, &local_col_indexes);

And finally call MatGetSubMatrices

Mat *local_submatrices;
MatGetSubMatrices(A, 1, &local_row_indexes, &local_col_indexes, MAT_INITIAL_MATRIX, &local_submatrices);