1: #include <petscmat.h>
  2: #include <petsc/private/matorderimpl.h>
  3: #include <petsc/private/matmetisimpl.h>
  4: #include <metis.h>

  6: /*
  7:     MatGetOrdering_METISND - Find the nested dissection ordering of a given matrix.
  8: */
  9: PETSC_INTERN PetscErrorCode MatGetOrdering_METISND(Mat mat, MatOrderingType type, IS *row, IS *col)
 10: {
 11:   PetscInt        i, j, iptr, ival, nrow, *xadj, *adjncy, *perm, *iperm;
 12:   const PetscInt *ia, *ja;
 13:   Mat             B = NULL;
 14:   idx_t           options[METIS_NOPTIONS];
 15:   PetscBool       done;

 17:   PetscFunctionBegin;
 18:   PetscCall(MatGetRowIJ(mat, 0, PETSC_TRUE, PETSC_TRUE, &nrow, &ia, &ja, &done));
 19:   if (!done) {
 20:     PetscCall(MatConvert(mat, MATSEQAIJ, MAT_INITIAL_MATRIX, &B));
 21:     PetscCall(MatGetRowIJ(B, 0, PETSC_TRUE, PETSC_TRUE, &nrow, &ia, &ja, &done));
 22:   }
 23:   PetscCallMETIS(METIS_SetDefaultOptions, options);
 24:   options[METIS_OPTION_NUMBERING] = 0;
 25:   PetscOptionsBegin(PetscObjectComm((PetscObject)mat), ((PetscObject)mat)->prefix, "METISND Options", "Mat");

 27:   ival = (PetscInt)options[METIS_OPTION_NSEPS];
 28:   PetscCall(PetscOptionsInt("-mat_ordering_metisnd_nseps", "number of different separators per level", "None", ival, &ival, NULL));
 29:   options[METIS_OPTION_NSEPS] = (idx_t)ival;

 31:   ival = (PetscInt)options[METIS_OPTION_NITER];
 32:   PetscCall(PetscOptionsInt("-mat_ordering_metisnd_niter", "number of refinement iterations", "None", ival, &ival, NULL));
 33:   options[METIS_OPTION_NITER] = (idx_t)ival;

 35:   ival = (PetscInt)options[METIS_OPTION_UFACTOR];
 36:   PetscCall(PetscOptionsInt("-mat_ordering_metisnd_ufactor", "maximum allowed imbalance", "None", ival, &ival, NULL));
 37:   options[METIS_OPTION_UFACTOR] = (idx_t)ival;

 39:   ival = (PetscInt)options[METIS_OPTION_PFACTOR];
 40:   PetscCall(PetscOptionsInt("-mat_ordering_metisnd_pfactor", "minimum degree of vertices that will be ordered last", "None", ival, &ival, NULL));
 41:   options[METIS_OPTION_PFACTOR] = (idx_t)ival;

 43:   PetscOptionsEnd();

 45:   PetscCall(PetscMalloc4(nrow + 1, &xadj, ia[nrow], &adjncy, nrow, &perm, nrow, &iperm));
 46:   /* The adjacency list of a vertex should not contain the vertex itself.
 47:   */
 48:   iptr       = 0;
 49:   xadj[iptr] = 0;
 50:   for (j = 0; j < nrow; j++) {
 51:     for (i = ia[j]; i < ia[j + 1]; i++) {
 52:       if (ja[i] != j) adjncy[iptr++] = ja[i];
 53:     }
 54:     xadj[j + 1] = iptr;
 55:   }

 57:   PetscCallMETIS(METIS_NodeND, &nrow, (idx_t *)xadj, (idx_t *)adjncy, NULL, options, (idx_t *)perm, (idx_t *)iperm);

 59:   if (B) {
 60:     PetscCall(MatRestoreRowIJ(B, 0, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done));
 61:     PetscCall(MatDestroy(&B));
 62:   } else {
 63:     PetscCall(MatRestoreRowIJ(mat, 0, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done));
 64:   }

 66:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nrow, perm, PETSC_COPY_VALUES, row));
 67:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nrow, perm, PETSC_COPY_VALUES, col));
 68:   PetscCall(PetscFree4(xadj, adjncy, perm, iperm));
 69:   PetscFunctionReturn(PETSC_SUCCESS);
 70: }