Actual source code: sorder.c
2: /*
3: Provides the code that allows PETSc users to register their own
4: sequential matrix Ordering routines.
5: */
6: #include <petsc/private/matimpl.h>
7: #include <petscmat.h>
9: PetscFunctionList MatOrderingList = NULL;
10: PetscBool MatOrderingRegisterAllCalled = PETSC_FALSE;
12: PetscErrorCode MatGetOrdering_Flow(Mat mat, MatOrderingType type, IS *irow, IS *icol)
13: {
14: PetscFunctionBegin;
15: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot do default flow ordering for matrix type");
16: }
18: PETSC_INTERN PetscErrorCode MatGetOrdering_Natural(Mat mat, MatOrderingType type, IS *irow, IS *icol)
19: {
20: PetscInt n, i, *ii;
21: PetscBool done;
22: MPI_Comm comm;
24: PetscFunctionBegin;
25: PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
26: PetscCall(MatGetRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, &n, NULL, NULL, &done));
27: PetscCall(MatRestoreRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, NULL, NULL, NULL, &done));
28: if (done) { /* matrix may be "compressed" in symbolic factorization, due to i-nodes or block storage */
29: /*
30: We actually create general index sets because this avoids mallocs to
31: to obtain the indices in the MatSolve() routines.
32: PetscCall(ISCreateStride(PETSC_COMM_SELF,n,0,1,irow));
33: PetscCall(ISCreateStride(PETSC_COMM_SELF,n,0,1,icol));
34: */
35: PetscCall(PetscMalloc1(n, &ii));
36: for (i = 0; i < n; i++) ii[i] = i;
37: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, ii, PETSC_COPY_VALUES, irow));
38: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, ii, PETSC_OWN_POINTER, icol));
39: } else {
40: PetscInt start, end;
42: PetscCall(MatGetOwnershipRange(mat, &start, &end));
43: PetscCall(ISCreateStride(comm, end - start, start, 1, irow));
44: PetscCall(ISCreateStride(comm, end - start, start, 1, icol));
45: }
46: PetscCall(ISSetIdentity(*irow));
47: PetscCall(ISSetIdentity(*icol));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: /*
52: Orders the rows (and columns) by the lengths of the rows.
53: This produces a symmetric Ordering but does not require a
54: matrix with symmetric non-zero structure.
55: */
56: PETSC_INTERN PetscErrorCode MatGetOrdering_RowLength(Mat mat, MatOrderingType type, IS *irow, IS *icol)
57: {
58: PetscInt n, *permr, *lens, i;
59: const PetscInt *ia, *ja;
60: PetscBool done;
62: PetscFunctionBegin;
63: PetscCall(MatGetRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, &n, &ia, &ja, &done));
64: PetscCheck(done, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot get rows for matrix");
66: PetscCall(PetscMalloc2(n, &lens, n, &permr));
67: for (i = 0; i < n; i++) {
68: lens[i] = ia[i + 1] - ia[i];
69: permr[i] = i;
70: }
71: PetscCall(MatRestoreRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, NULL, &ia, &ja, &done));
73: PetscCall(PetscSortIntWithPermutation(n, lens, permr));
75: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, permr, PETSC_COPY_VALUES, irow));
76: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, permr, PETSC_COPY_VALUES, icol));
77: PetscCall(PetscFree2(lens, permr));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: /*@C
82: MatOrderingRegister - Adds a new sparse matrix ordering to the matrix package.
84: Not Collective
86: Input Parameters:
87: + sname - name of ordering (for example `MATORDERINGND`)
88: - function - function pointer that creates the ordering
90: Level: developer
92: Sample usage:
93: .vb
94: MatOrderingRegister("my_order", MyOrder);
95: .ve
97: Then, your partitioner can be chosen with the procedural interface via
98: $ MatOrderingSetType(part, "my_order)
99: or at runtime via the option
100: $ -pc_factor_mat_ordering_type my_order
102: .seealso: `MatOrderingRegisterAll()`, `MatGetOrdering()`
103: @*/
104: PetscErrorCode MatOrderingRegister(const char sname[], PetscErrorCode (*function)(Mat, MatOrderingType, IS *, IS *))
105: {
106: PetscFunctionBegin;
107: PetscCall(MatInitializePackage());
108: PetscCall(PetscFunctionListAdd(&MatOrderingList, sname, function));
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: #include <../src/mat/impls/aij/mpi/mpiaij.h>
113: /*@C
114: MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
115: improve numerical stability of LU factorization.
117: Collective
119: Input Parameters:
120: + mat - the matrix
121: - type - type of reordering, one of the following
122: .vb
123: MATORDERINGNATURAL_OR_ND - Nested dissection unless matrix is SBAIJ then it is natural
124: MATORDERINGNATURAL - Natural
125: MATORDERINGND - Nested Dissection
126: MATORDERING1WD - One-way Dissection
127: MATORDERINGRCM - Reverse Cuthill-McKee
128: MATORDERINGQMD - Quotient Minimum Degree
129: MATORDERINGEXTERNAL - Use an ordering internal to the factorzation package and do not compute or use PETSc's
130: .ve
132: Output Parameters:
133: + rperm - row permutation indices
134: - cperm - column permutation indices
136: Options Database Key:
137: + -mat_view_ordering draw - plots matrix nonzero structure in new ordering
138: - -pc_factor_mat_ordering_type <nd,natural,..> - ordering to use with `PC`s based on factorization, `MATLU`, `MATILU`, MATCHOLESKY`, `MATICC`
140: Level: intermediate
142: Notes:
143: This DOES NOT actually reorder the matrix; it merely returns two index sets
144: that define a reordering. This is usually not used directly, rather use the
145: options `PCFactorSetMatOrderingType()`
147: The user can define additional orderings; see `MatOrderingRegister()`.
149: These are generally only implemented for sequential sparse matrices.
151: Some external packages that PETSc can use for direct factorization such as SuperLU_DIST do not accept orderings provided by
152: this call.
154: If `MATORDERINGEXTERNAL` is used then PETSc does not compute an ordering and utilizes one built into the factorization package
156: .seealso: `MatOrderingRegister()`, `PCFactorSetMatOrderingType()`, `MatColoring`, `MatColoringCreate()`, `MatOrderingType`, `Mat`
157: @*/
158: PetscErrorCode MatGetOrdering(Mat mat, MatOrderingType type, IS *rperm, IS *cperm)
159: {
160: PetscInt mmat, nmat, mis;
161: PetscErrorCode (*r)(Mat, MatOrderingType, IS *, IS *);
162: PetscBool flg, ismpiaij;
164: PetscFunctionBegin;
168: PetscCheck(mat->assembled, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
169: PetscCheck(!mat->factortype, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
170: PetscCheck(type, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Ordering type cannot be null");
172: PetscCall(PetscStrcmp(type, MATORDERINGEXTERNAL, &flg));
173: if (flg) {
174: *rperm = NULL;
175: *cperm = NULL;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: /* This code is terrible. MatGetOrdering() multiple dispatch should use matrix and this code should move to impls/aij/mpi. */
180: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIAIJ, &ismpiaij));
181: if (ismpiaij) { /* Reorder using diagonal block */
182: Mat Ad, Ao;
183: const PetscInt *colmap;
184: IS lrowperm, lcolperm;
185: PetscInt i, rstart, rend, *idx;
186: const PetscInt *lidx;
188: PetscCall(MatMPIAIJGetSeqAIJ(mat, &Ad, &Ao, &colmap));
189: PetscCall(MatGetOrdering(Ad, type, &lrowperm, &lcolperm));
190: PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
191: /* Remap row index set to global space */
192: PetscCall(ISGetIndices(lrowperm, &lidx));
193: PetscCall(PetscMalloc1(rend - rstart, &idx));
194: for (i = 0; i + rstart < rend; i++) idx[i] = rstart + lidx[i];
195: PetscCall(ISRestoreIndices(lrowperm, &lidx));
196: PetscCall(ISDestroy(&lrowperm));
197: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), rend - rstart, idx, PETSC_OWN_POINTER, rperm));
198: PetscCall(ISSetPermutation(*rperm));
199: /* Remap column index set to global space */
200: PetscCall(ISGetIndices(lcolperm, &lidx));
201: PetscCall(PetscMalloc1(rend - rstart, &idx));
202: for (i = 0; i + rstart < rend; i++) idx[i] = rstart + lidx[i];
203: PetscCall(ISRestoreIndices(lcolperm, &lidx));
204: PetscCall(ISDestroy(&lcolperm));
205: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), rend - rstart, idx, PETSC_OWN_POINTER, cperm));
206: PetscCall(ISSetPermutation(*cperm));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: if (!mat->rmap->N) { /* matrix has zero rows */
211: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, cperm));
212: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, rperm));
213: PetscCall(ISSetIdentity(*cperm));
214: PetscCall(ISSetIdentity(*rperm));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: PetscCall(MatGetLocalSize(mat, &mmat, &nmat));
219: PetscCheck(mmat == nmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must be square matrix, rows %" PetscInt_FMT " columns %" PetscInt_FMT, mmat, nmat);
221: PetscCall(MatOrderingRegisterAll());
222: PetscCall(PetscFunctionListFind(MatOrderingList, type, &r));
223: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown or unregistered type: %s", type);
225: PetscCall(PetscLogEventBegin(MAT_GetOrdering, mat, 0, 0, 0));
226: PetscCall((*r)(mat, type, rperm, cperm));
227: PetscCall(ISSetPermutation(*rperm));
228: PetscCall(ISSetPermutation(*cperm));
229: /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
230: PetscCall(ISGetLocalSize(*rperm, &mis));
231: if (mmat > mis) PetscCall(MatInodeAdjustForInodes(mat, rperm, cperm));
232: PetscCall(PetscLogEventEnd(MAT_GetOrdering, mat, 0, 0, 0));
234: PetscCall(PetscOptionsHasName(((PetscObject)mat)->options, ((PetscObject)mat)->prefix, "-mat_view_ordering", &flg));
235: if (flg) {
236: Mat tmat;
237: PetscCall(MatPermute(mat, *rperm, *cperm, &tmat));
238: PetscCall(MatViewFromOptions(tmat, (PetscObject)mat, "-mat_view_ordering"));
239: PetscCall(MatDestroy(&tmat));
240: }
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
245: {
246: PetscFunctionBegin;
247: *list = MatOrderingList;
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }