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: }