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:   MatGetRowIJ(mat, 0, PETSC_TRUE, PETSC_TRUE, &nrow, &ia, &ja, &done);
 19:   if (!done) {
 20:     MatConvert(mat, MATSEQAIJ, MAT_INITIAL_MATRIX, &B);
 21:     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:   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:   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:   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:   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:   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:     MatRestoreRowIJ(B, 0, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done);
 73:     MatDestroy(&B);
 74:   } else {
 75:     MatRestoreRowIJ(mat, 0, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done);
 76:   }

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