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 <nd,natural,..> - ordering to use with `PC`s based on factorization, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
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: PetscInt mmat, nmat, mis;
152: PetscErrorCode (*r)(Mat, MatOrderingType, IS *, IS *);
153: PetscBool flg, ismpiaij;
155: PetscFunctionBegin;
157: PetscAssertPointer(rperm, 3);
158: PetscAssertPointer(cperm, 4);
159: PetscCheck(mat->assembled, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
160: PetscCheck(!mat->factortype, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
161: PetscCheck(type, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Ordering type cannot be null");
163: PetscCall(PetscStrcmp(type, MATORDERINGEXTERNAL, &flg));
164: if (flg) {
165: *rperm = NULL;
166: *cperm = NULL;
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: /* This code is terrible. MatGetOrdering() multiple dispatch should use matrix and this code should move to impls/aij/mpi. */
171: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIAIJ, &ismpiaij));
172: if (ismpiaij) { /* Reorder using diagonal block */
173: Mat Ad, Ao;
174: const PetscInt *colmap;
175: IS lrowperm, lcolperm;
176: PetscInt i, rstart, rend, *idx;
177: const PetscInt *lidx;
179: PetscCall(MatMPIAIJGetSeqAIJ(mat, &Ad, &Ao, &colmap));
180: PetscCall(MatGetOrdering(Ad, type, &lrowperm, &lcolperm));
181: PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
182: /* Remap row index set to global space */
183: PetscCall(ISGetIndices(lrowperm, &lidx));
184: PetscCall(PetscMalloc1(rend - rstart, &idx));
185: for (i = 0; i + rstart < rend; i++) idx[i] = rstart + lidx[i];
186: PetscCall(ISRestoreIndices(lrowperm, &lidx));
187: PetscCall(ISDestroy(&lrowperm));
188: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), rend - rstart, idx, PETSC_OWN_POINTER, rperm));
189: PetscCall(ISSetPermutation(*rperm));
190: /* Remap column index set to global space */
191: PetscCall(ISGetIndices(lcolperm, &lidx));
192: PetscCall(PetscMalloc1(rend - rstart, &idx));
193: for (i = 0; i + rstart < rend; i++) idx[i] = rstart + lidx[i];
194: PetscCall(ISRestoreIndices(lcolperm, &lidx));
195: PetscCall(ISDestroy(&lcolperm));
196: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), rend - rstart, idx, PETSC_OWN_POINTER, cperm));
197: PetscCall(ISSetPermutation(*cperm));
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: if (!mat->rmap->N) { /* matrix has zero rows */
202: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, cperm));
203: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, rperm));
204: PetscCall(ISSetIdentity(*cperm));
205: PetscCall(ISSetIdentity(*rperm));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: PetscCall(MatGetLocalSize(mat, &mmat, &nmat));
210: PetscCheck(mmat == nmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must be square matrix, rows %" PetscInt_FMT " columns %" PetscInt_FMT, mmat, nmat);
212: PetscCall(MatOrderingRegisterAll());
213: PetscCall(PetscFunctionListFind(MatOrderingList, type, &r));
214: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown or unregistered type: %s", type);
216: PetscCall(PetscLogEventBegin(MAT_GetOrdering, mat, 0, 0, 0));
217: PetscCall((*r)(mat, type, rperm, cperm));
218: PetscCall(ISSetPermutation(*rperm));
219: PetscCall(ISSetPermutation(*cperm));
220: /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
221: PetscCall(ISGetLocalSize(*rperm, &mis));
222: if (mmat > mis) PetscCall(MatInodeAdjustForInodes(mat, rperm, cperm));
223: PetscCall(PetscLogEventEnd(MAT_GetOrdering, mat, 0, 0, 0));
225: PetscCall(PetscOptionsHasName(((PetscObject)mat)->options, ((PetscObject)mat)->prefix, "-mat_view_ordering", &flg));
226: if (flg) {
227: Mat tmat;
228: PetscCall(MatPermute(mat, *rperm, *cperm, &tmat));
229: PetscCall(MatViewFromOptions(tmat, (PetscObject)mat, "-mat_view_ordering"));
230: PetscCall(MatDestroy(&tmat));
231: }
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
236: {
237: PetscFunctionBegin;
238: *list = MatOrderingList;
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }