Actual source code: metisnd.c


  2: #include <petscmat.h>
  3: #include <petsc/private/matorderimpl.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:   int             status;
 14:   Mat             B = NULL;
 15:   idx_t           options[METIS_NOPTIONS];
 16:   PetscBool       done;

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

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

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

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

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

 44:   PetscOptionsEnd();

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

 58:   status = METIS_NodeND(&nrow, (idx_t *)xadj, (idx_t *)adjncy, NULL, options, (idx_t *)perm, (idx_t *)iperm);
 59:   switch (status) {
 60:   case METIS_OK:
 61:     break;
 62:   case METIS_ERROR:
 63:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_LIB, "METIS returned with an unspecified error");
 64:   case METIS_ERROR_INPUT:
 65:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_LIB, "METIS received an invalid input");
 66:   case METIS_ERROR_MEMORY:
 67:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_MEM, "METIS could not compute ordering");
 68:   default:
 69:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_LIB, "Unexpected return value");
 70:   }

 72:   if (B) {
 73:     PetscCall(MatRestoreRowIJ(B, 0, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done));
 74:     PetscCall(MatDestroy(&B));
 75:   } else {
 76:     PetscCall(MatRestoreRowIJ(mat, 0, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done));
 77:   }

 79:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nrow, perm, PETSC_COPY_VALUES, row));
 80:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nrow, perm, PETSC_COPY_VALUES, col));
 81:   PetscCall(PetscFree4(xadj, adjncy, perm, iperm));
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }