3
$\begingroup$

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)

$\endgroup$
2
$\begingroup$

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);
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.