Actual source code: spnd.c
1: #include <petscmat.h>
2: #include <petsc/private/matorderimpl.h>
4: /*
5: MatGetOrdering_ND - Find the nested dissection ordering of a given matrix.
6: */
7: PETSC_INTERN PetscErrorCode MatGetOrdering_ND(Mat mat, MatOrderingType type, IS *row, IS *col)
8: {
9: PetscInt i, *mask, *xls, *ls, nrow, *perm;
10: const PetscInt *ia, *ja;
11: PetscBool done;
12: Mat B = NULL;
14: PetscFunctionBegin;
15: PetscCall(MatGetRowIJ(mat, 1, PETSC_TRUE, PETSC_TRUE, &nrow, &ia, &ja, &done));
16: if (!done) {
17: PetscCall(MatConvert(mat, MATAIJ, MAT_INITIAL_MATRIX, &B));
18: PetscCall(MatGetRowIJ(B, 1, PETSC_TRUE, PETSC_TRUE, &nrow, &ia, &ja, &done));
19: }
21: PetscCall(PetscMalloc4(nrow, &mask, nrow, &perm, nrow + 1, &xls, nrow, &ls));
22: PetscCall(SPARSEPACKgennd(&nrow, ia, ja, mask, perm, xls, ls));
23: if (B) {
24: PetscCall(MatRestoreRowIJ(B, 1, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done));
25: PetscCall(MatDestroy(&B));
26: } else {
27: PetscCall(MatRestoreRowIJ(mat, 1, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done));
28: }
30: /* shift because Sparsepack indices start at one */
31: for (i = 0; i < nrow; i++) perm[i]--;
33: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nrow, perm, PETSC_COPY_VALUES, row));
34: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nrow, perm, PETSC_COPY_VALUES, col));
35: PetscCall(PetscFree4(mask, perm, xls, ls));
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }