1: #include <../src/mat/impls/adj/mpi/mpiadj.h>
  2: #include <petsc/private/matmetisimpl.h>
  3: #include <metis.h>

  5: PETSC_EXTERN PetscErrorCode MatMeshToCellGraph_METIS(Mat mesh, PetscInt ncommonnodes, Mat *dual)
  6: {
  7:   PetscInt   *newxadj, *newadjncy;
  8:   PetscInt    numflag = 0;
  9:   Mat_MPIAdj *adj     = (Mat_MPIAdj *)mesh->data, *newadj;
 10:   PetscBool   flg;
 11:   MPI_Comm    comm;
 12:   PetscMPIInt size;

 14:   PetscFunctionBegin;
 15:   PetscCall(PetscObjectTypeCompare((PetscObject)mesh, MATMPIADJ, &flg));
 16:   PetscCall(PetscObjectGetComm((PetscObject)mesh, &comm));
 17:   PetscCheck(flg, comm, PETSC_ERR_SUP, "Must use MPIAdj matrix type");

 19:   PetscCallMPI(MPI_Comm_size(comm, &size));
 20:   PetscCheck(size == 1, comm, PETSC_ERR_WRONG_MPI_SIZE, "MatMeshToCellGraph_METIS() requires a sequential matrix (communicator size must be 1)");

 22:   {
 23:     idx_t ne = mesh->rmap->N;
 24:     idx_t nn = mesh->cmap->N;

 26:     PetscCallMETIS(METIS_MeshToDual, &ne, &nn, (idx_t *)adj->i, (idx_t *)adj->j, (idx_t *)&ncommonnodes, (idx_t *)&numflag, (idx_t **)&newxadj, (idx_t **)&newadjncy);
 27:   }

 29:   for (PetscInt i = 0; i < mesh->rmap->N; i++) PetscCall(PetscSortInt(newxadj[i + 1] - newxadj[i], newadjncy + newxadj[i]));

 31:   PetscCall(MatCreateMPIAdj(PetscObjectComm((PetscObject)mesh), mesh->rmap->n, mesh->rmap->N, newxadj, newadjncy, NULL, dual));
 32:   newadj = (Mat_MPIAdj *)(*dual)->data;

 34:   newadj->freeaijwithfree = PETSC_TRUE; /* signal the matrix should be freed with system free since space was allocated by METIS */
 35:   PetscFunctionReturn(PETSC_SUCCESS);
 36: }