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
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
91: $ MatOrderingSetType(part, "my_order)
92: or at runtime via the option
93: $ -pc_factor_mat_ordering_type my_order
95: .seealso: `MatOrderingRegisterAll()`, `MatGetOrdering()`
96: @*/
97: PetscErrorCode MatOrderingRegister(const char sname[], PetscErrorCode (*function)(Mat, MatOrderingType, IS *, IS *))
98: {
99: PetscFunctionBegin;
100: PetscCall(MatInitializePackage());
101: PetscCall(PetscFunctionListAdd(&MatOrderingList, sname, function));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: #include <../src/mat/impls/aij/mpi/mpiaij.h>
106: /*@C
107: MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
108: improve numerical stability of LU factorization.
110: Collective
112: Input Parameters:
113: + mat - the matrix
114: - type - type of reordering, one of the following
115: .vb
116: MATORDERINGNATURAL_OR_ND - Nested dissection unless matrix is SBAIJ then it is natural
117: MATORDERINGNATURAL - Natural
118: MATORDERINGND - Nested Dissection
119: MATORDERING1WD - One-way Dissection
120: MATORDERINGRCM - Reverse Cuthill-McKee
121: MATORDERINGQMD - Quotient Minimum Degree
122: MATORDERINGEXTERNAL - Use an ordering internal to the factorzation package and do not compute or use PETSc's
123: .ve
125: Output Parameters:
126: + rperm - row permutation indices
127: - cperm - column permutation indices
129: Options Database Key:
130: + -mat_view_ordering draw - plots matrix nonzero structure in new ordering
131: - -pc_factor_mat_ordering_type <nd,natural,..> - ordering to use with `PC`s based on factorization, `MATLU`, `MATILU`, MATCHOLESKY`, `MATICC`
133: Level: intermediate
135: Notes:
136: This DOES NOT actually reorder the matrix; it merely returns two index sets
137: that define a reordering. This is usually not used directly, rather use the
138: options `PCFactorSetMatOrderingType()`
140: The user can define additional orderings; see `MatOrderingRegister()`.
142: These are generally only implemented for sequential sparse matrices.
144: Some external packages that PETSc can use for direct factorization such as SuperLU_DIST do not accept orderings provided by
145: this call.
147: If `MATORDERINGEXTERNAL` is used then PETSc does not compute an ordering and utilizes one built into the factorization package
149: .seealso: `MatOrderingRegister()`, `PCFactorSetMatOrderingType()`, `MatColoring`, `MatColoringCreate()`, `MatOrderingType`, `Mat`
150: @*/
151: PetscErrorCode MatGetOrdering(Mat mat, MatOrderingType type, IS *rperm, IS *cperm)
152: {
153: PetscInt mmat, nmat, mis;
154: PetscErrorCode (*r)(Mat, MatOrderingType, IS *, IS *);
155: PetscBool flg, ismpiaij;
157: PetscFunctionBegin;
159: PetscAssertPointer(rperm, 3);
160: PetscAssertPointer(cperm, 4);
161: PetscCheck(mat->assembled, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
162: PetscCheck(!mat->factortype, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
163: PetscCheck(type, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Ordering type cannot be null");
165: PetscCall(PetscStrcmp(type, MATORDERINGEXTERNAL, &flg));
166: if (flg) {
167: *rperm = NULL;
168: *cperm = NULL;
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: /* This code is terrible. MatGetOrdering() multiple dispatch should use matrix and this code should move to impls/aij/mpi. */
173: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIAIJ, &ismpiaij));
174: if (ismpiaij) { /* Reorder using diagonal block */
175: Mat Ad, Ao;
176: const PetscInt *colmap;
177: IS lrowperm, lcolperm;
178: PetscInt i, rstart, rend, *idx;
179: const PetscInt *lidx;
181: PetscCall(MatMPIAIJGetSeqAIJ(mat, &Ad, &Ao, &colmap));
182: PetscCall(MatGetOrdering(Ad, type, &lrowperm, &lcolperm));
183: PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
184: /* Remap row index set to global space */
185: PetscCall(ISGetIndices(lrowperm, &lidx));
186: PetscCall(PetscMalloc1(rend - rstart, &idx));
187: for (i = 0; i + rstart < rend; i++) idx[i] = rstart + lidx[i];
188: PetscCall(ISRestoreIndices(lrowperm, &lidx));
189: PetscCall(ISDestroy(&lrowperm));
190: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), rend - rstart, idx, PETSC_OWN_POINTER, rperm));
191: PetscCall(ISSetPermutation(*rperm));
192: /* Remap column index set to global space */
193: PetscCall(ISGetIndices(lcolperm, &lidx));
194: PetscCall(PetscMalloc1(rend - rstart, &idx));
195: for (i = 0; i + rstart < rend; i++) idx[i] = rstart + lidx[i];
196: PetscCall(ISRestoreIndices(lcolperm, &lidx));
197: PetscCall(ISDestroy(&lcolperm));
198: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), rend - rstart, idx, PETSC_OWN_POINTER, cperm));
199: PetscCall(ISSetPermutation(*cperm));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: if (!mat->rmap->N) { /* matrix has zero rows */
204: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, cperm));
205: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, rperm));
206: PetscCall(ISSetIdentity(*cperm));
207: PetscCall(ISSetIdentity(*rperm));
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: PetscCall(MatGetLocalSize(mat, &mmat, &nmat));
212: PetscCheck(mmat == nmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must be square matrix, rows %" PetscInt_FMT " columns %" PetscInt_FMT, mmat, nmat);
214: PetscCall(MatOrderingRegisterAll());
215: PetscCall(PetscFunctionListFind(MatOrderingList, type, &r));
216: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown or unregistered type: %s", type);
218: PetscCall(PetscLogEventBegin(MAT_GetOrdering, mat, 0, 0, 0));
219: PetscCall((*r)(mat, type, rperm, cperm));
220: PetscCall(ISSetPermutation(*rperm));
221: PetscCall(ISSetPermutation(*cperm));
222: /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
223: PetscCall(ISGetLocalSize(*rperm, &mis));
224: if (mmat > mis) PetscCall(MatInodeAdjustForInodes(mat, rperm, cperm));
225: PetscCall(PetscLogEventEnd(MAT_GetOrdering, mat, 0, 0, 0));
227: PetscCall(PetscOptionsHasName(((PetscObject)mat)->options, ((PetscObject)mat)->prefix, "-mat_view_ordering", &flg));
228: if (flg) {
229: Mat tmat;
230: PetscCall(MatPermute(mat, *rperm, *cperm, &tmat));
231: PetscCall(MatViewFromOptions(tmat, (PetscObject)mat, "-mat_view_ordering"));
232: PetscCall(MatDestroy(&tmat));
233: }
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
238: {
239: PetscFunctionBegin;
240: *list = MatOrderingList;
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }