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