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