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: }