Actual source code: sorder.c

petsc-master 2017-10-18
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              = 0;
 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: #if !defined(PETSC_USE_DEBUG)
 19:   return(0);
 20: #endif
 21: }

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

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

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

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

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

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

 80:   PetscSortIntWithPermutation(n,lens,permr);

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

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

 91:    Not Collective

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

 97:    Level: developer

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

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

109: .keywords: matrix, ordering, register

111: .seealso: MatOrderingRegisterDestroy(), MatOrderingRegisterAll()
112: @*/
113: PetscErrorCode  MatOrderingRegister(const char sname[],PetscErrorCode (*function)(Mat,MatOrderingType,IS*,IS*))
114: {

118:   PetscFunctionListAdd(&MatOrderingList,sname,function);
119:   return(0);
120: }

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

127:    Collective on Mat

129:    Input Parameters:
130: +  mat - the matrix
131: -  type - type of reordering, one of the following:
132: $      MATORDERINGNATURAL - Natural
133: $      MATORDERINGND - Nested Dissection
134: $      MATORDERING1WD - One-way Dissection
135: $      MATORDERINGRCM - Reverse Cuthill-McKee
136: $      MATORDERINGQMD - Quotient Minimum Degree

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


143:    Options Database Key:
144: . -mat_view_ordering draw - plots matrix nonzero structure in new ordering

146:    Level: intermediate

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

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

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

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


161: .keywords: matrix, set, ordering, factorization, direct, ILU, LU,
162:            fill, reordering, natural, Nested Dissection,
163:            One-way Dissection, Cholesky, Reverse Cuthill-McKee,
164:            Quotient Minimum Degree

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

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

182:   /* This code is terrible. MatGetOrdering() multiple dispatch should use matrix and this code should move to impls/aij/mpi. */
183:   PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJ,&ismpiaij);
184:   if (ismpiaij) {               /* Reorder using diagonal block */
185:     Mat            Ad,Ao;
186:     const PetscInt *colmap;
187:     IS             lrowperm,lcolperm;
188:     PetscInt       i,rstart,rend,*idx;
189:     const PetscInt *lidx;

191:     MatMPIAIJGetSeqAIJ(mat,&Ad,&Ao,&colmap);
192:     MatGetOrdering(Ad,type,&lrowperm,&lcolperm);
193:     MatGetOwnershipRange(mat,&rstart,&rend);
194:     /* Remap row index set to global space */
195:     ISGetIndices(lrowperm,&lidx);
196:     PetscMalloc1(rend-rstart,&idx);
197:     for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
198:     ISRestoreIndices(lrowperm,&lidx);
199:     ISDestroy(&lrowperm);
200:     ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,rperm);
201:     ISSetPermutation(*rperm);
202:     /* Remap column index set to global space */
203:     ISGetIndices(lcolperm,&lidx);
204:     PetscMalloc1(rend-rstart,&idx);
205:     for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
206:     ISRestoreIndices(lcolperm,&lidx);
207:     ISDestroy(&lcolperm);
208:     ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,cperm);
209:     ISSetPermutation(*cperm);
210:     return(0);
211:   }

213:   /* this chunk of code is REALLY bad, should maybe get the ordering from the factor matrix,
214:      then those that don't support orderings will handle their cases themselves. */
215:   PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);
216:   PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense);
217:   PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJCUSP,&ismpiaijcusp);
218:   PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJCUSPARSE,&ismpiaijcusparse);
219:   PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&ismpibaij);
220:   PetscObjectTypeCompare((PetscObject)mat,MATMPISBAIJ,&ismpisbaij);
221:   PetscObjectTypeCompare((PetscObject)mat,MATELEMENTAL,&iselemental);
222:   if (isseqdense || ismpidense || ismpibaij || ismpisbaij || ismpiaijcusp || ismpiaijcusparse || iselemental) {
223:     MatGetLocalSize(mat,&m,NULL);
224:     /*
225:        These matrices only give natural ordering
226:     */
227:     ISCreateStride(PETSC_COMM_SELF,m,0,1,cperm);
228:     ISCreateStride(PETSC_COMM_SELF,m,0,1,rperm);
229:     ISSetIdentity(*cperm);
230:     ISSetIdentity(*rperm);
231:     return(0);
232:   }

234:   if (!mat->rmap->N) { /* matrix has zero rows */
235:     ISCreateStride(PETSC_COMM_SELF,0,0,1,cperm);
236:     ISCreateStride(PETSC_COMM_SELF,0,0,1,rperm);
237:     ISSetIdentity(*cperm);
238:     ISSetIdentity(*rperm);
239:     return(0);
240:   }

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

245:   MatOrderingRegisterAll();
246:   PetscFunctionListFind(MatOrderingList,type,&r);
247:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);

249:   PetscLogEventBegin(MAT_GetOrdering,mat,0,0,0);
250:   (*r)(mat,type,rperm,cperm);
251:   ISSetPermutation(*rperm);
252:   ISSetPermutation(*cperm);
253:   /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
254:   ISGetLocalSize(*rperm,&mis);
255:   if (mmat > mis) {MatInodeAdjustForInodes(mat,rperm,cperm);}
256:   PetscLogEventEnd(MAT_GetOrdering,mat,0,0,0);


259:   PetscOptionsHasName(((PetscObject)mat)->options,((PetscObject)mat)->prefix,"-mat_view_ordering",&flg);
260:   if (flg) {
261:     Mat tmat;
262:     MatPermute(mat,*rperm,*cperm,&tmat);
263:     MatViewFromOptions(tmat,(PetscObject)mat,"-mat_view_ordering");
264:     MatDestroy(&tmat);
265:   }
266:   return(0);
267: }

269: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
270: {
272:   *list = MatOrderingList;
273:   return(0);
274: }