Actual source code: metisnd.c

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

  5: /*
  6:     MatGetOrdering_METISND - Find the nested dissection ordering of a given matrix.
  7: */
  8: PETSC_INTERN PetscErrorCode MatGetOrdering_METISND(Mat mat, MatOrderingType type, IS *row, IS *col)
  9: {
 10:   PetscInt        i, j, iptr, ival, nrow, *xadj, *adjncy, *perm, *iperm;
 11:   const PetscInt *ia, *ja;
 12:   int             status;
 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:   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:   status = METIS_NodeND(&nrow, (idx_t *)xadj, (idx_t *)adjncy, NULL, options, (idx_t *)perm, (idx_t *)iperm);
 58:   switch (status) {
 59:   case METIS_OK:
 60:     break;
 61:   case METIS_ERROR:
 62:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_LIB, "METIS returned with an unspecified error");
 63:   case METIS_ERROR_INPUT:
 64:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_LIB, "METIS received an invalid input");
 65:   case METIS_ERROR_MEMORY:
 66:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_MEM, "METIS could not compute ordering");
 67:   default:
 68:     SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_LIB, "Unexpected return value");
 69:   }

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

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