1*8be712e4SBarry Smith #include <petscmat.h>
2*8be712e4SBarry Smith #include <petsc/private/matorderimpl.h>
3*8be712e4SBarry Smith #include <metis.h>
4*8be712e4SBarry Smith
5*8be712e4SBarry Smith /*
6*8be712e4SBarry Smith MatGetOrdering_METISND - Find the nested dissection ordering of a given matrix.
7*8be712e4SBarry Smith */
MatGetOrdering_METISND(Mat mat,MatOrderingType type,IS * row,IS * col)8*8be712e4SBarry Smith PETSC_INTERN PetscErrorCode MatGetOrdering_METISND(Mat mat, MatOrderingType type, IS *row, IS *col)
9*8be712e4SBarry Smith {
10*8be712e4SBarry Smith PetscInt i, j, iptr, ival, nrow, *xadj, *adjncy, *perm, *iperm;
11*8be712e4SBarry Smith const PetscInt *ia, *ja;
12*8be712e4SBarry Smith int status;
13*8be712e4SBarry Smith Mat B = NULL;
14*8be712e4SBarry Smith idx_t options[METIS_NOPTIONS];
15*8be712e4SBarry Smith PetscBool done;
16*8be712e4SBarry Smith
17*8be712e4SBarry Smith PetscFunctionBegin;
18*8be712e4SBarry Smith PetscCall(MatGetRowIJ(mat, 0, PETSC_TRUE, PETSC_TRUE, &nrow, &ia, &ja, &done));
19*8be712e4SBarry Smith if (!done) {
20*8be712e4SBarry Smith PetscCall(MatConvert(mat, MATSEQAIJ, MAT_INITIAL_MATRIX, &B));
21*8be712e4SBarry Smith PetscCall(MatGetRowIJ(B, 0, PETSC_TRUE, PETSC_TRUE, &nrow, &ia, &ja, &done));
22*8be712e4SBarry Smith }
23*8be712e4SBarry Smith METIS_SetDefaultOptions(options);
24*8be712e4SBarry Smith options[METIS_OPTION_NUMBERING] = 0;
25*8be712e4SBarry Smith PetscOptionsBegin(PetscObjectComm((PetscObject)mat), ((PetscObject)mat)->prefix, "METISND Options", "Mat");
26*8be712e4SBarry Smith
27*8be712e4SBarry Smith ival = (PetscInt)options[METIS_OPTION_NSEPS];
28*8be712e4SBarry Smith PetscCall(PetscOptionsInt("-mat_ordering_metisnd_nseps", "number of different separators per level", "None", ival, &ival, NULL));
29*8be712e4SBarry Smith options[METIS_OPTION_NSEPS] = (idx_t)ival;
30*8be712e4SBarry Smith
31*8be712e4SBarry Smith ival = (PetscInt)options[METIS_OPTION_NITER];
32*8be712e4SBarry Smith PetscCall(PetscOptionsInt("-mat_ordering_metisnd_niter", "number of refinement iterations", "None", ival, &ival, NULL));
33*8be712e4SBarry Smith options[METIS_OPTION_NITER] = (idx_t)ival;
34*8be712e4SBarry Smith
35*8be712e4SBarry Smith ival = (PetscInt)options[METIS_OPTION_UFACTOR];
36*8be712e4SBarry Smith PetscCall(PetscOptionsInt("-mat_ordering_metisnd_ufactor", "maximum allowed imbalance", "None", ival, &ival, NULL));
37*8be712e4SBarry Smith options[METIS_OPTION_UFACTOR] = (idx_t)ival;
38*8be712e4SBarry Smith
39*8be712e4SBarry Smith ival = (PetscInt)options[METIS_OPTION_PFACTOR];
40*8be712e4SBarry Smith PetscCall(PetscOptionsInt("-mat_ordering_metisnd_pfactor", "minimum degree of vertices that will be ordered last", "None", ival, &ival, NULL));
41*8be712e4SBarry Smith options[METIS_OPTION_PFACTOR] = (idx_t)ival;
42*8be712e4SBarry Smith
43*8be712e4SBarry Smith PetscOptionsEnd();
44*8be712e4SBarry Smith
45*8be712e4SBarry Smith PetscCall(PetscMalloc4(nrow + 1, &xadj, ia[nrow], &adjncy, nrow, &perm, nrow, &iperm));
46*8be712e4SBarry Smith /* The adjacency list of a vertex should not contain the vertex itself.
47*8be712e4SBarry Smith */
48*8be712e4SBarry Smith iptr = 0;
49*8be712e4SBarry Smith xadj[iptr] = 0;
50*8be712e4SBarry Smith for (j = 0; j < nrow; j++) {
51*8be712e4SBarry Smith for (i = ia[j]; i < ia[j + 1]; i++) {
52*8be712e4SBarry Smith if (ja[i] != j) adjncy[iptr++] = ja[i];
53*8be712e4SBarry Smith }
54*8be712e4SBarry Smith xadj[j + 1] = iptr;
55*8be712e4SBarry Smith }
56*8be712e4SBarry Smith
57*8be712e4SBarry Smith status = METIS_NodeND(&nrow, (idx_t *)xadj, (idx_t *)adjncy, NULL, options, (idx_t *)perm, (idx_t *)iperm);
58*8be712e4SBarry Smith switch (status) {
59*8be712e4SBarry Smith case METIS_OK:
60*8be712e4SBarry Smith break;
61*8be712e4SBarry Smith case METIS_ERROR:
62*8be712e4SBarry Smith SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_LIB, "METIS returned with an unspecified error");
63*8be712e4SBarry Smith case METIS_ERROR_INPUT:
64*8be712e4SBarry Smith SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_LIB, "METIS received an invalid input");
65*8be712e4SBarry Smith case METIS_ERROR_MEMORY:
66*8be712e4SBarry Smith SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_MEM, "METIS could not compute ordering");
67*8be712e4SBarry Smith default:
68*8be712e4SBarry Smith SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_LIB, "Unexpected return value");
69*8be712e4SBarry Smith }
70*8be712e4SBarry Smith
71*8be712e4SBarry Smith if (B) {
72*8be712e4SBarry Smith PetscCall(MatRestoreRowIJ(B, 0, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done));
73*8be712e4SBarry Smith PetscCall(MatDestroy(&B));
74*8be712e4SBarry Smith } else {
75*8be712e4SBarry Smith PetscCall(MatRestoreRowIJ(mat, 0, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done));
76*8be712e4SBarry Smith }
77*8be712e4SBarry Smith
78*8be712e4SBarry Smith PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nrow, perm, PETSC_COPY_VALUES, row));
79*8be712e4SBarry Smith PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nrow, perm, PETSC_COPY_VALUES, col));
80*8be712e4SBarry Smith PetscCall(PetscFree4(xadj, adjncy, perm, iperm));
81*8be712e4SBarry Smith PetscFunctionReturn(PETSC_SUCCESS);
82*8be712e4SBarry Smith }
83