1: #include <../src/mat/impls/adj/mpi/mpiadj.h>

  3: #include <metis.h>

  5: #define PetscCallMETIS(n, func) \
  6:   do { \
  7:     PetscCheck(n != METIS_ERROR_INPUT, PETSC_COMM_SELF, PETSC_ERR_LIB, "METIS error due to wrong inputs and/or options for %s", func); \
  8:     PetscCheck(n != METIS_ERROR_MEMORY, PETSC_COMM_SELF, PETSC_ERR_LIB, "METIS error due to insufficient memory in %s", func); \
  9:     PetscCheck(n != METIS_ERROR, PETSC_COMM_SELF, PETSC_ERR_LIB, "METIS general error in %s", func); \
 10:   } while (0)

 12: #define PetscCallMetis_(name, func, args) \
 13:   do { \
 14:     PetscStackPushExternal(name); \
 15:     int status = func args; \
 16:     PetscStackPop; \
 17:     PetscCallMETIS(status, name); \
 18:   } while (0)

 20: #define PetscCallMetis(func, args) PetscCallMetis_(PetscStringize(func), func, args)

 22: PETSC_EXTERN PetscErrorCode MatMeshToCellGraph_Metis(Mat mesh, PetscInt ncommonnodes, Mat *dual)
 23: {
 24:   PetscInt   *newxadj, *newadjncy;
 25:   PetscInt    numflag = 0;
 26:   Mat_MPIAdj *adj     = (Mat_MPIAdj *)mesh->data, *newadj;
 27:   PetscBool   flg;
 28:   MPI_Comm    comm;
 29:   PetscMPIInt size;

 31:   PetscFunctionBegin;
 32:   PetscCall(PetscObjectTypeCompare((PetscObject)mesh, MATMPIADJ, &flg));
 33:   PetscCall(PetscObjectGetComm((PetscObject)mesh, &comm));
 34:   PetscCheck(flg, comm, PETSC_ERR_SUP, "Must use MPIAdj matrix type");

 36:   PetscCallMPI(MPI_Comm_size(comm, &size));
 37:   PetscCheck(size == 1, comm, PETSC_ERR_SUP, "MatMeshToCellGraph_Metis() requires a sequential matrix (communicator size must be 1)");

 39:   {
 40:     idx_t ne = mesh->rmap->N;
 41:     idx_t nn = mesh->cmap->N;

 43:     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));
 44:   }

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

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

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