# PETSc - Manipulate BAIJ matrix locally

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

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)

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