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 */ 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