Actual source code: sorder.c

petsc-master 2020-09-24
Report Typos and Errors

  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: {
 17:   SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot do default flow ordering for matrix type");
 18:   return(0);
 19: }

 21: PETSC_INTERN PetscErrorCode MatGetOrdering_Natural(Mat mat,MatOrderingType type,IS *irow,IS *icol)
 22: {
 24:   PetscInt       n,i,*ii;
 25:   PetscBool      done;
 26:   MPI_Comm       comm;

 29:   PetscObjectGetComm((PetscObject)mat,&comm);
 30:   MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,NULL,NULL,&done);
 31:   MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,NULL,NULL,&done);
 32:   if (done) { /* matrix may be "compressed" in symbolic factorization, due to i-nodes or block storage */
 33:     /*
 34:       We actually create general index sets because this avoids mallocs to
 35:       to obtain the indices in the MatSolve() routines.
 36:       ISCreateStride(PETSC_COMM_SELF,n,0,1,irow);
 37:       ISCreateStride(PETSC_COMM_SELF,n,0,1,icol);
 38:     */
 39:     PetscMalloc1(n,&ii);
 40:     for (i=0; i<n; i++) ii[i] = i;
 41:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_COPY_VALUES,irow);
 42:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,icol);
 43:   } else {
 44:     PetscInt start,end;

 46:     MatGetOwnershipRange(mat,&start,&end);
 47:     ISCreateStride(comm,end-start,start,1,irow);
 48:     ISCreateStride(comm,end-start,start,1,icol);
 49:   }
 50:   ISSetIdentity(*irow);
 51:   ISSetIdentity(*icol);
 52:   return(0);
 53: }

 55: /*
 56:      Orders the rows (and columns) by the lengths of the rows.
 57:    This produces a symmetric Ordering but does not require a
 58:    matrix with symmetric non-zero structure.
 59: */
 60: PETSC_INTERN PetscErrorCode MatGetOrdering_RowLength(Mat mat,MatOrderingType type,IS *irow,IS *icol)
 61: {
 63:   PetscInt       n,*permr,*lens,i;
 64:   const PetscInt *ia,*ja;
 65:   PetscBool      done;

 68:   MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,&ia,&ja,&done);
 69:   if (!done) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get rows for matrix");

 71:   PetscMalloc2(n,&lens,n,&permr);
 72:   for (i=0; i<n; i++) {
 73:     lens[i]  = ia[i+1] - ia[i];
 74:     permr[i] = i;
 75:   }
 76:   MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,&ia,&ja,&done);

 78:   PetscSortIntWithPermutation(n,lens,permr);

 80:   ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,irow);
 81:   ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,icol);
 82:   PetscFree2(lens,permr);
 83:   return(0);
 84: }

 86: /*@C
 87:    MatOrderingRegister - Adds a new sparse matrix ordering to the matrix package.

 89:    Not Collective

 91:    Input Parameters:
 92: +  sname - name of ordering (for example MATORDERINGND)
 93: -  function - function pointer that creates the ordering

 95:    Level: developer

 97:    Sample usage:
 98: .vb
 99:    MatOrderingRegister("my_order", MyOrder);
100: .ve

102:    Then, your partitioner can be chosen with the procedural interface via
103: $     MatOrderingSetType(part,"my_order)
104:    or at runtime via the option
105: $     -pc_factor_mat_ordering_type my_order

107: .seealso: MatOrderingRegisterDestroy(), MatOrderingRegisterAll()
108: @*/
109: PetscErrorCode  MatOrderingRegister(const char sname[],PetscErrorCode (*function)(Mat,MatOrderingType,IS*,IS*))
110: {

114:   MatInitializePackage();
115:   PetscFunctionListAdd(&MatOrderingList,sname,function);
116:   return(0);
117: }

119: #include <../src/mat/impls/aij/mpi/mpiaij.h>
120: /*@C
121:    MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
122:    improve numerical stability of LU factorization.

124:    Collective on Mat

126:    Input Parameters:
127: +  mat - the matrix
128: -  type - type of reordering, one of the following:
129: $      MATORDERINGNATURAL_OR_ND - Nested dissection unless matrix is SBAIJ then it is natural
130: $      MATORDERINGNATURAL - Natural
131: $      MATORDERINGND - Nested Dissection
132: $      MATORDERING1WD - One-way Dissection
133: $      MATORDERINGRCM - Reverse Cuthill-McKee
134: $      MATORDERINGQMD - Quotient Minimum Degree
135: $      MATORDERINGEXTERNAL - Use an ordering internal to the factorzation package and do not compute or use PETSc's

137:    Output Parameters:
138: +  rperm - row permutation indices
139: -  cperm - column permutation indices

141:    Options Database Key:
142: + -mat_view_ordering draw - plots matrix nonzero structure in new ordering
143: - -pc_factor_mat_ordering_type <nd,natural,..> - ordering to use with PCs based on factorization, LU, ILU, Cholesky, ICC

145:    Level: intermediate

147:    Notes:
148:       This DOES NOT actually reorder the matrix; it merely returns two index sets
149:    that define a reordering. This is usually not used directly, rather use the
150:    options PCFactorSetMatOrderingType()

152:    The user can define additional orderings; see MatOrderingRegister().

154:    These are generally only implemented for sequential sparse matrices.

156:    Some external packages that PETSc can use for direct factorization such as SuperLU do not accept orderings provided by
157:    this call.

159:    If MATORDERINGEXTERNAL is used then PETSc does not compute an ordering and utilizes one built into the factorization package

161:            fill, reordering, natural, Nested Dissection,
162:            One-way Dissection, Cholesky, Reverse Cuthill-McKee,
163:            Quotient Minimum Degree

165: .seealso:   MatOrderingRegister(), PCFactorSetMatOrderingType(), MatColoring, MatColoringCreate()
166: @*/
167: PetscErrorCode  MatGetOrdering(Mat mat,MatOrderingType type,IS *rperm,IS *cperm)
168: {
170:   PetscInt       mmat,nmat,mis,m;
171:   PetscErrorCode (*r)(Mat,MatOrderingType,IS*,IS*);
172:   PetscBool      flg = PETSC_FALSE,isseqdense,ismpidense,ismpiaij,ismpibaij,ismpisbaij,ismpiaijcusparse,iselemental,isscalapack,flg1;

178:   if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
179:   if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

181:   PetscStrcmp(type,MATORDERINGEXTERNAL,&flg1);
182:   if (flg1) {
183:     *rperm = NULL;
184:     *cperm = NULL;
185:     return(0);
186:   }

188:   PetscStrcmp(type,MATORDERINGNATURAL_OR_ND,&flg1);
189:   if (flg1) {
190:     PetscBool isseqsbaij;
191:     PetscObjectTypeCompareAny((PetscObject)mat,&isseqsbaij,MATSEQSBAIJ,MATSEQBAIJ,NULL);
192:     if (isseqsbaij) {
193:       type = MATORDERINGNATURAL;
194:     } else {
195:       type = MATORDERINGND;
196:     }
197:   }

199:   /* This code is terrible. MatGetOrdering() multiple dispatch should use matrix and this code should move to impls/aij/mpi. */
200:   PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJ,&ismpiaij);
201:   if (ismpiaij) {               /* Reorder using diagonal block */
202:     Mat            Ad,Ao;
203:     const PetscInt *colmap;
204:     IS             lrowperm,lcolperm;
205:     PetscInt       i,rstart,rend,*idx;
206:     const PetscInt *lidx;

208:     MatMPIAIJGetSeqAIJ(mat,&Ad,&Ao,&colmap);
209:     MatGetOrdering(Ad,type,&lrowperm,&lcolperm);
210:     MatGetOwnershipRange(mat,&rstart,&rend);
211:     /* Remap row index set to global space */
212:     ISGetIndices(lrowperm,&lidx);
213:     PetscMalloc1(rend-rstart,&idx);
214:     for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
215:     ISRestoreIndices(lrowperm,&lidx);
216:     ISDestroy(&lrowperm);
217:     ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,rperm);
218:     ISSetPermutation(*rperm);
219:     /* Remap column index set to global space */
220:     ISGetIndices(lcolperm,&lidx);
221:     PetscMalloc1(rend-rstart,&idx);
222:     for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
223:     ISRestoreIndices(lcolperm,&lidx);
224:     ISDestroy(&lcolperm);
225:     ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,cperm);
226:     ISSetPermutation(*cperm);
227:     return(0);
228:   }

230:   /* this chunk of code is REALLY bad, should maybe get the ordering from the factor matrix,
231:      then those that don't support orderings will handle their cases themselves. */
232:   PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);
233:   PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense);
234:   PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJCUSPARSE,&ismpiaijcusparse);
235:   PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&ismpibaij);
236:   PetscObjectTypeCompare((PetscObject)mat,MATMPISBAIJ,&ismpisbaij);
237:   PetscObjectTypeCompare((PetscObject)mat,MATELEMENTAL,&iselemental);
238:   PetscObjectTypeCompare((PetscObject)mat,MATSCALAPACK,&isscalapack);
239:   if (isseqdense || ismpidense || ismpibaij || ismpisbaij || ismpiaijcusparse || iselemental || isscalapack) {
240:     MatGetLocalSize(mat,&m,NULL);
241:     /*
242:        These matrices only give natural ordering
243:     */
244:     ISCreateStride(PETSC_COMM_SELF,m,0,1,cperm);
245:     ISCreateStride(PETSC_COMM_SELF,m,0,1,rperm);
246:     ISSetIdentity(*cperm);
247:     ISSetIdentity(*rperm);
248:     return(0);
249:   }

251:   if (!mat->rmap->N) { /* matrix has zero rows */
252:     ISCreateStride(PETSC_COMM_SELF,0,0,1,cperm);
253:     ISCreateStride(PETSC_COMM_SELF,0,0,1,rperm);
254:     ISSetIdentity(*cperm);
255:     ISSetIdentity(*rperm);
256:     return(0);
257:   }

259:   MatGetLocalSize(mat,&mmat,&nmat);
260:   if (mmat != nmat) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",mmat,nmat);

262:   MatOrderingRegisterAll();
263:   PetscFunctionListFind(MatOrderingList,type,&r);
264:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);

266:   PetscLogEventBegin(MAT_GetOrdering,mat,0,0,0);
267:   (*r)(mat,type,rperm,cperm);
268:   ISSetPermutation(*rperm);
269:   ISSetPermutation(*cperm);
270:   /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
271:   ISGetLocalSize(*rperm,&mis);
272:   if (mmat > mis) {MatInodeAdjustForInodes(mat,rperm,cperm);}
273:   PetscLogEventEnd(MAT_GetOrdering,mat,0,0,0);


276:   PetscOptionsHasName(((PetscObject)mat)->options,((PetscObject)mat)->prefix,"-mat_view_ordering",&flg);
277:   if (flg) {
278:     Mat tmat;
279:     MatPermute(mat,*rperm,*cperm,&tmat);
280:     MatViewFromOptions(tmat,(PetscObject)mat,"-mat_view_ordering");
281:     MatDestroy(&tmat);
282:   }
283:   return(0);
284: }

286: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
287: {
289:   *list = MatOrderingList;
290:   return(0);
291: }