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