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