1 #include <petscmat.h>
2 #include <petsc/private/matorderimpl.h>
3 #include <metis.h>
4
5 /*
6 MatGetOrdering_METISND - Find the nested dissection ordering of a given matrix.
7 */
MatGetOrdering_METISND(Mat mat,MatOrderingType type,IS * row,IS * col)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;
16
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");
26
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;
30
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;
34
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;
38
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;
42
43 PetscOptionsEnd();
44
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 }
56
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 }
70
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 }
77
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 }
83