Actual source code: sorder.c
1: /*
2: Provides the code that allows PETSc users to register their own
3: sequential matrix Ordering routines.
4: */
5: #include <petsc/private/matimpl.h>
6: #include <petscmat.h>
8: PetscFunctionList MatOrderingList = NULL;
9: PetscBool MatOrderingRegisterAllCalled = PETSC_FALSE;
11: PETSC_INTERN PetscErrorCode MatGetOrdering_Natural(Mat mat, MatOrderingType type, IS *irow, IS *icol)
12: {
13: PetscInt n, i, *ii;
14: PetscBool done;
15: MPI_Comm comm;
17: PetscFunctionBegin;
18: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
19: PetscCall(MatGetRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, &n, NULL, NULL, &done));
20: PetscCall(MatRestoreRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, NULL, NULL, NULL, &done));
21: if (done) { /* matrix may be "compressed" in symbolic factorization, due to i-nodes or block storage */
22: /*
23: We actually create general index sets because this avoids mallocs to
24: to obtain the indices in the MatSolve() routines.
25: PetscCall(ISCreateStride(PETSC_COMM_SELF,n,0,1,irow));
26: PetscCall(ISCreateStride(PETSC_COMM_SELF,n,0,1,icol));
27: */
28: PetscCall(PetscMalloc1(n, &ii));
29: for (i = 0; i < n; i++) ii[i] = i;
30: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, ii, PETSC_COPY_VALUES, irow));
31: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, ii, PETSC_OWN_POINTER, icol));
32: } else {
33: PetscInt start, end;
35: PetscCall(MatGetOwnershipRange(mat, &start, &end));
36: PetscCall(ISCreateStride(comm, end - start, start, 1, irow));
37: PetscCall(ISCreateStride(comm, end - start, start, 1, icol));
38: }
39: PetscCall(ISSetIdentity(*irow));
40: PetscCall(ISSetIdentity(*icol));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: /*
45: Orders the rows (and columns) by the lengths of the rows.
46: This produces a symmetric Ordering but does not require a
47: matrix with symmetric non-zero structure.
48: */
49: PETSC_INTERN PetscErrorCode MatGetOrdering_RowLength(Mat mat, MatOrderingType type, IS *irow, IS *icol)
50: {
51: PetscInt n, *permr, *lens, i;
52: const PetscInt *ia, *ja;
53: PetscBool done;
55: PetscFunctionBegin;
56: PetscCall(MatGetRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, &n, &ia, &ja, &done));
57: PetscCheck(done, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot get rows for matrix");
59: PetscCall(PetscMalloc2(n, &lens, n, &permr));
60: for (i = 0; i < n; i++) {
61: lens[i] = ia[i + 1] - ia[i];
62: permr[i] = i;
63: }
64: PetscCall(MatRestoreRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, NULL, &ia, &ja, &done));
66: PetscCall(PetscSortIntWithPermutation(n, lens, permr));
68: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, permr, PETSC_COPY_VALUES, irow));
69: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, permr, PETSC_COPY_VALUES, icol));
70: PetscCall(PetscFree2(lens, permr));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: /*@C
75: MatOrderingRegister - Adds a new sparse matrix ordering to the matrix package.
77: Not Collective, No Fortran Support
79: Input Parameters:
80: + sname - name of ordering (for example `MATORDERINGND`)
81: - function - function pointer that creates the ordering
83: Level: developer
85: Example Usage:
86: .vb
87: MatOrderingRegister("my_order", MyOrder);
88: .ve
90: Then, your partitioner can be chosen with the procedural interface via `MatOrderingSetType(part, "my_order)` or at runtime via the option
91: `-pc_factor_mat_ordering_type my_order`
93: .seealso: `Mat`, `MatOrderingType`, `MatOrderingRegisterAll()`, `MatGetOrdering()`
94: @*/
95: PetscErrorCode MatOrderingRegister(const char sname[], PetscErrorCode (*function)(Mat, MatOrderingType, IS *, IS *))
96: {
97: PetscFunctionBegin;
98: PetscCall(MatInitializePackage());
99: PetscCall(PetscFunctionListAdd(&MatOrderingList, sname, function));
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: #include <../src/mat/impls/aij/mpi/mpiaij.h>
104: /*@
105: MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
106: improve numerical stability of LU factorization.
108: Collective
110: Input Parameters:
111: + mat - the matrix
112: - type - type of reordering, one of the following
113: .vb
114: MATORDERINGNATURAL_OR_ND - Nested dissection unless matrix is SBAIJ then it is natural
115: MATORDERINGNATURAL - Natural
116: MATORDERINGND - Nested Dissection
117: MATORDERING1WD - One-way Dissection
118: MATORDERINGRCM - Reverse Cuthill-McKee
119: MATORDERINGQMD - Quotient Minimum Degree
120: MATORDERINGEXTERNAL - Use an ordering internal to the factorzation package and do not compute or use PETSc's
121: .ve
123: Output Parameters:
124: + rperm - row permutation indices
125: - cperm - column permutation indices
127: Options Database Key:
128: + -mat_view_ordering draw - plots matrix nonzero structure in new ordering
129: - -pc_factor_mat_ordering_type type - ordering to use with `PC`s based on factorization see `MatOrderingType`
131: Level: intermediate
133: Notes:
134: This DOES NOT actually reorder the matrix; it merely returns two index sets
135: that define a reordering. This is usually not used directly, rather use the
136: options `PCFactorSetMatOrderingType()`
138: The user can define additional orderings; see `MatOrderingRegister()`.
140: These are generally only implemented for sequential sparse matrices.
142: Some external packages that PETSc can use for direct factorization such as SuperLU_DIST do not accept orderings provided by
143: this call.
145: If `MATORDERINGEXTERNAL` is used then PETSc does not compute an ordering and utilizes one built into the factorization package
147: .seealso: `MatOrderingRegister()`, `PCFactorSetMatOrderingType()`, `MatColoring`, `MatColoringCreate()`, `MatOrderingType`, `Mat`
148: @*/
149: PetscErrorCode MatGetOrdering(Mat mat, MatOrderingType type, IS *rperm, IS *cperm)
150: {
151: PetscBool flg;
153: PetscFunctionBegin;
155: PetscAssertPointer(rperm, 3);
156: PetscAssertPointer(cperm, 4);
157: PetscCheck(mat->assembled, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
158: PetscCheck(!mat->factortype, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
159: PetscCheck(type, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Ordering type cannot be null");
161: PetscCall(PetscStrcmp(type, MATORDERINGEXTERNAL, &flg));
162: if (flg) {
163: *rperm = NULL;
164: *cperm = NULL;
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: if (!mat->rmap->N) { /* matrix has zero rows */
169: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, cperm));
170: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, rperm));
171: PetscCall(ISSetIdentity(*cperm));
172: PetscCall(ISSetIdentity(*rperm));
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: if (!mat->ops->getordering) {
177: PetscInt mmat, nmat, mis;
178: PetscErrorCode (*r)(Mat, MatOrderingType, IS *, IS *);
180: PetscCall(MatGetLocalSize(mat, &mmat, &nmat));
181: PetscCheck(mmat == nmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must be square matrix, rows %" PetscInt_FMT " columns %" PetscInt_FMT, mmat, nmat);
183: PetscCall(MatOrderingRegisterAll());
184: PetscCall(PetscFunctionListFind(MatOrderingList, type, &r));
185: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown or unregistered type: %s", type);
187: PetscCall(PetscLogEventBegin(MAT_GetOrdering, mat, 0, 0, 0));
188: PetscCall((*r)(mat, type, rperm, cperm));
189: PetscCall(ISSetPermutation(*rperm));
190: PetscCall(ISSetPermutation(*cperm));
191: /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
192: PetscCall(ISGetLocalSize(*rperm, &mis));
193: if (mmat > mis) PetscCall(MatInodeAdjustForInodes(mat, rperm, cperm));
194: PetscCall(PetscLogEventEnd(MAT_GetOrdering, mat, 0, 0, 0));
196: PetscCall(PetscOptionsHasName(((PetscObject)mat)->options, ((PetscObject)mat)->prefix, "-mat_view_ordering", &flg));
197: if (flg) {
198: Mat tmat;
199: PetscCall(MatPermute(mat, *rperm, *cperm, &tmat));
200: PetscCall(MatViewFromOptions(tmat, (PetscObject)mat, "-mat_view_ordering"));
201: PetscCall(MatDestroy(&tmat));
202: }
203: } else PetscUseTypeMethod(mat, getordering, type, rperm, cperm);
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
208: {
209: PetscFunctionBegin;
210: *list = MatOrderingList;
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }