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: extern PetscErrorCode MatGetOrdering_Flow_SeqAIJ(Mat, MatOrderingType, IS *, IS *);
14: PetscErrorCode MatGetOrdering_Flow(Mat mat, MatOrderingType type, IS *irow, IS *icol)
15: {
16: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot do default flow ordering for matrix type");
17: }
19: PETSC_INTERN PetscErrorCode MatGetOrdering_Natural(Mat mat, MatOrderingType type, IS *irow, IS *icol)
20: {
21: PetscInt n, i, *ii;
22: PetscBool done;
23: MPI_Comm comm;
25: PetscObjectGetComm((PetscObject)mat, &comm);
26: MatGetRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, &n, NULL, NULL, &done);
27: 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: ISCreateStride(PETSC_COMM_SELF,n,0,1,irow);
33: ISCreateStride(PETSC_COMM_SELF,n,0,1,icol);
34: */
35: PetscMalloc1(n, &ii);
36: for (i = 0; i < n; i++) ii[i] = i;
37: ISCreateGeneral(PETSC_COMM_SELF, n, ii, PETSC_COPY_VALUES, irow);
38: ISCreateGeneral(PETSC_COMM_SELF, n, ii, PETSC_OWN_POINTER, icol);
39: } else {
40: PetscInt start, end;
42: MatGetOwnershipRange(mat, &start, &end);
43: ISCreateStride(comm, end - start, start, 1, irow);
44: ISCreateStride(comm, end - start, start, 1, icol);
45: }
46: ISSetIdentity(*irow);
47: ISSetIdentity(*icol);
48: return 0;
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: MatGetRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, &n, &ia, &ja, &done);
65: PetscMalloc2(n, &lens, n, &permr);
66: for (i = 0; i < n; i++) {
67: lens[i] = ia[i + 1] - ia[i];
68: permr[i] = i;
69: }
70: MatRestoreRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, NULL, &ia, &ja, &done);
72: PetscSortIntWithPermutation(n, lens, permr);
74: ISCreateGeneral(PETSC_COMM_SELF, n, permr, PETSC_COPY_VALUES, irow);
75: ISCreateGeneral(PETSC_COMM_SELF, n, permr, PETSC_COPY_VALUES, icol);
76: PetscFree2(lens, permr);
77: return 0;
78: }
80: /*@C
81: MatOrderingRegister - Adds a new sparse matrix ordering to the matrix package.
83: Not Collective
85: Input Parameters:
86: + sname - name of ordering (for example `MATORDERINGND`)
87: - function - function pointer that creates the ordering
89: Level: developer
91: Sample usage:
92: .vb
93: MatOrderingRegister("my_order", MyOrder);
94: .ve
96: Then, your partitioner can be chosen with the procedural interface via
97: $ MatOrderingSetType(part,"my_order)
98: or at runtime via the option
99: $ -pc_factor_mat_ordering_type my_order
101: .seealso: `MatOrderingRegisterAll()`, `MatGetOrdering()`
102: @*/
103: PetscErrorCode MatOrderingRegister(const char sname[], PetscErrorCode (*function)(Mat, MatOrderingType, IS *, IS *))
104: {
105: MatInitializePackage();
106: PetscFunctionListAdd(&MatOrderingList, sname, function);
107: return 0;
108: }
110: #include <../src/mat/impls/aij/mpi/mpiaij.h>
111: /*@C
112: MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
113: improve numerical stability of LU factorization.
115: Collective
117: Input Parameters:
118: + mat - the matrix
119: - type - type of reordering, one of the following
121: .vb
122: MATORDERINGNATURAL_OR_ND - Nested dissection unless matrix is SBAIJ then it is natural
123: MATORDERINGNATURAL - Natural
124: MATORDERINGND - Nested Dissection
125: MATORDERING1WD - One-way Dissection
126: MATORDERINGRCM - Reverse Cuthill-McKee
127: MATORDERINGQMD - Quotient Minimum Degree
128: MATORDERINGEXTERNAL - Use an ordering internal to the factorzation package and do not compute or use PETSc's
129: .ve
131: Output Parameters:
132: + rperm - row permutation indices
133: - cperm - column permutation indices
135: Options Database Key:
136: + -mat_view_ordering draw - plots matrix nonzero structure in new ordering
137: - -pc_factor_mat_ordering_type <nd,natural,..> - ordering to use with `PC`s based on factorization, `MATLU`, `MATILU`, MATCHOLESKY`, `MATICC`
139: Level: intermediate
141: Notes:
142: This DOES NOT actually reorder the matrix; it merely returns two index sets
143: that define a reordering. This is usually not used directly, rather use the
144: options `PCFactorSetMatOrderingType()`
146: The user can define additional orderings; see MatOrderingRegister().
148: These are generally only implemented for sequential sparse matrices.
150: Some external packages that PETSc can use for direct factorization such as SuperLU do not accept orderings provided by
151: this call.
153: If `MATORDERINGEXTERNAL` is used then PETSc does not compute an ordering and utilizes one built into the factorization package
155: .seealso: `MatOrderingRegister()`, `PCFactorSetMatOrderingType()`, `MatColoring`, `MatColoringCreate()`, `MatOrderingType`, `Mat`
156: @*/
157: PetscErrorCode MatGetOrdering(Mat mat, MatOrderingType type, IS *rperm, IS *cperm)
158: {
159: PetscInt mmat, nmat, mis;
160: PetscErrorCode (*r)(Mat, MatOrderingType, IS *, IS *);
161: PetscBool flg, ismpiaij;
170: PetscStrcmp(type, MATORDERINGEXTERNAL, &flg);
171: if (flg) {
172: *rperm = NULL;
173: *cperm = NULL;
174: return 0;
175: }
177: /* This code is terrible. MatGetOrdering() multiple dispatch should use matrix and this code should move to impls/aij/mpi. */
178: PetscObjectTypeCompare((PetscObject)mat, MATMPIAIJ, &ismpiaij);
179: if (ismpiaij) { /* Reorder using diagonal block */
180: Mat Ad, Ao;
181: const PetscInt *colmap;
182: IS lrowperm, lcolperm;
183: PetscInt i, rstart, rend, *idx;
184: const PetscInt *lidx;
186: MatMPIAIJGetSeqAIJ(mat, &Ad, &Ao, &colmap);
187: MatGetOrdering(Ad, type, &lrowperm, &lcolperm);
188: MatGetOwnershipRange(mat, &rstart, &rend);
189: /* Remap row index set to global space */
190: ISGetIndices(lrowperm, &lidx);
191: PetscMalloc1(rend - rstart, &idx);
192: for (i = 0; i + rstart < rend; i++) idx[i] = rstart + lidx[i];
193: ISRestoreIndices(lrowperm, &lidx);
194: ISDestroy(&lrowperm);
195: ISCreateGeneral(PetscObjectComm((PetscObject)mat), rend - rstart, idx, PETSC_OWN_POINTER, rperm);
196: ISSetPermutation(*rperm);
197: /* Remap column index set to global space */
198: ISGetIndices(lcolperm, &lidx);
199: PetscMalloc1(rend - rstart, &idx);
200: for (i = 0; i + rstart < rend; i++) idx[i] = rstart + lidx[i];
201: ISRestoreIndices(lcolperm, &lidx);
202: ISDestroy(&lcolperm);
203: ISCreateGeneral(PetscObjectComm((PetscObject)mat), rend - rstart, idx, PETSC_OWN_POINTER, cperm);
204: ISSetPermutation(*cperm);
205: return 0;
206: }
208: if (!mat->rmap->N) { /* matrix has zero rows */
209: ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, cperm);
210: ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, rperm);
211: ISSetIdentity(*cperm);
212: ISSetIdentity(*rperm);
213: return 0;
214: }
216: MatGetLocalSize(mat, &mmat, &nmat);
219: MatOrderingRegisterAll();
220: PetscFunctionListFind(MatOrderingList, type, &r);
223: PetscLogEventBegin(MAT_GetOrdering, mat, 0, 0, 0);
224: (*r)(mat, type, rperm, cperm);
225: ISSetPermutation(*rperm);
226: ISSetPermutation(*cperm);
227: /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
228: ISGetLocalSize(*rperm, &mis);
229: if (mmat > mis) MatInodeAdjustForInodes(mat, rperm, cperm);
230: PetscLogEventEnd(MAT_GetOrdering, mat, 0, 0, 0);
232: PetscOptionsHasName(((PetscObject)mat)->options, ((PetscObject)mat)->prefix, "-mat_view_ordering", &flg);
233: if (flg) {
234: Mat tmat;
235: MatPermute(mat, *rperm, *cperm, &tmat);
236: MatViewFromOptions(tmat, (PetscObject)mat, "-mat_view_ordering");
237: MatDestroy(&tmat);
238: }
239: return 0;
240: }
242: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
243: {
244: *list = MatOrderingList;
245: return 0;
246: }