Actual source code: sorder.c

petsc-3.3-p7 2013-05-11
  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>  /*I "petscmat.h" I*/

  9: PetscFList      MatOrderingList = 0;
 10: PetscBool  MatOrderingRegisterAllCalled = PETSC_FALSE;

 12: extern PetscErrorCode MatGetOrdering_Flow_SeqAIJ(Mat,const MatOrderingType,IS *,IS *);

 16: PetscErrorCode MatGetOrdering_Flow(Mat mat,const MatOrderingType type,IS *irow,IS *icol)
 17: {
 19:   SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Cannot do default flow ordering for matrix type");
 20: #if !defined(PETSC_USE_DEBUG)
 21:   return(0);
 22: #endif
 23: }

 25: EXTERN_C_BEGIN
 28: PetscErrorCode  MatGetOrdering_Natural(Mat mat,const MatOrderingType type,IS *irow,IS *icol)
 29: {
 31:   PetscInt       n,i,*ii;
 32:   PetscBool      done;
 33:   MPI_Comm       comm;

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

 53:     MatGetOwnershipRange(mat,&start,&end);
 54:     ISCreateStride(comm,end-start,start,1,irow);
 55:     ISCreateStride(comm,end-start,start,1,icol);
 56:   }
 57:   ISSetIdentity(*irow);
 58:   ISSetIdentity(*icol);
 59:   return(0);
 60: }
 61: EXTERN_C_END

 63: EXTERN_C_BEGIN
 64: /*
 65:      Orders the rows (and columns) by the lengths of the rows. 
 66:    This produces a symmetric Ordering but does not require a 
 67:    matrix with symmetric non-zero structure.
 68: */
 71: PetscErrorCode  MatGetOrdering_RowLength(Mat mat,const MatOrderingType type,IS *irow,IS *icol)
 72: {
 74:   PetscInt       n,*ia,*ja,*permr,*lens,i;
 75:   PetscBool      done;

 78:   MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,&ia,&ja,&done);
 79:   if (!done) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_SUP,"Cannot get rows for matrix");

 81:   PetscMalloc2(n,PetscInt,&lens,n,PetscInt,&permr);
 82:   for (i=0; i<n; i++) {
 83:     lens[i]  = ia[i+1] - ia[i];
 84:     permr[i] = i;
 85:   }
 86:   MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,&ia,&ja,&done);

 88:   PetscSortIntWithPermutation(n,lens,permr);

 90:   ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,irow);
 91:   ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,icol);
 92:   PetscFree2(lens,permr);
 93:   return(0);
 94: }
 95: EXTERN_C_END

 99: PetscErrorCode  MatOrderingRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(Mat,const MatOrderingType,IS*,IS*))
100: {
102:   char           fullname[PETSC_MAX_PATH_LEN];

105:   PetscFListConcat(path,name,fullname);
106:   PetscFListAdd(&MatOrderingList,sname,fullname,(void (*)(void))function);
107:   return(0);
108: }

112: /*@
113:    MatOrderingRegisterDestroy - Frees the list of ordering routines.

115:    Not collective

117:    Level: developer
118:    
119: .keywords: matrix, register, destroy

121: .seealso: MatOrderingRegisterDynamic(), MatOrderingRegisterAll()
122: @*/
123: PetscErrorCode  MatOrderingRegisterDestroy(void)
124: {

128:   PetscFListDestroy(&MatOrderingList);
129:   MatOrderingRegisterAllCalled = PETSC_FALSE;
130:   return(0);
131: }

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

140:    Collective on Mat

142:    Input Parameters:
143: +  mat - the matrix
144: -  type - type of reordering, one of the following:
145: $      MATORDERINGNATURAL - Natural
146: $      MATORDERINGND - Nested Dissection
147: $      MATORDERING1WD - One-way Dissection
148: $      MATORDERINGRCM - Reverse Cuthill-McKee
149: $      MATORDERINGQMD - Quotient Minimum Degree

151:    Output Parameters:
152: +  rperm - row permutation indices
153: -  cperm - column permutation indices


156:    Options Database Key:
157: . -mat_view_ordering_draw - plots matrix nonzero structure in new ordering

159:    Level: intermediate
160:    
161:    Notes:
162:       This DOES NOT actually reorder the matrix; it merely returns two index sets
163:    that define a reordering. This is usually not used directly, rather use the 
164:    options PCFactorSetMatOrderingType()

166:    The user can define additional orderings; see MatOrderingRegisterDynamic().

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

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


174: .keywords: matrix, set, ordering, factorization, direct, ILU, LU,
175:            fill, reordering, natural, Nested Dissection,
176:            One-way Dissection, Cholesky, Reverse Cuthill-McKee, 
177:            Quotient Minimum Degree

179: .seealso:   MatOrderingRegisterDynamic(), PCFactorSetMatOrderingType()
180: @*/
181: PetscErrorCode  MatGetOrdering(Mat mat,const MatOrderingType type,IS *rperm,IS *cperm)
182: {
184:   PetscInt       mmat,nmat,mis,m;
185:   PetscErrorCode (*r)(Mat,const MatOrderingType,IS*,IS*);
186:   PetscBool      flg = PETSC_FALSE,isseqdense,ismpidense,ismpiaij,ismpibaij,ismpisbaij,ismpiaijcusp;

192:   if (!mat->assembled) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
193:   if (mat->factortype) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

195:   /* this chunk of code is REALLY bad, should maybe get the ordering from the factor matrix,
196:      then those that don't support orderings will handle their cases themselfs. */
197:   PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);
198:   PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense);
199:   PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJ,&ismpiaij);
200:   PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJCUSP,&ismpiaijcusp);
201:   PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&ismpibaij);
202:   PetscObjectTypeCompare((PetscObject)mat,MATMPISBAIJ,&ismpisbaij);
203:   if (isseqdense || ismpidense || ismpiaij || ismpibaij || ismpisbaij || ismpiaijcusp) {
204:     MatGetLocalSize(mat,&m,PETSC_NULL);
205:     /*
206:        Dense matrices only give natural ordering
207:     */
208:     ISCreateStride(PETSC_COMM_SELF,0,m,1,cperm);
209:     ISCreateStride(PETSC_COMM_SELF,0,m,1,rperm);
210:     ISSetIdentity(*cperm);
211:     ISSetIdentity(*rperm);
212:     ISSetPermutation(*rperm);
213:     ISSetPermutation(*cperm);
214:     return(0);
215:   }

217:   if (!mat->rmap->N) { /* matrix has zero rows */
218:     ISCreateStride(PETSC_COMM_SELF,0,0,1,cperm);
219:     ISCreateStride(PETSC_COMM_SELF,0,0,1,rperm);
220:     ISSetIdentity(*cperm);
221:     ISSetIdentity(*rperm);
222:     ISSetPermutation(*rperm);
223:     ISSetPermutation(*cperm);
224:     return(0);
225:   }
226: 
227:   MatGetLocalSize(mat,&mmat,&nmat);
228:   if (mmat != nmat) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",mmat,nmat);

230:   if (!MatOrderingRegisterAllCalled) {MatOrderingRegisterAll(PETSC_NULL);}
231:   PetscFListFind(MatOrderingList,((PetscObject)mat)->comm,type,PETSC_TRUE,(void (**)(void)) &r);
232:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);

234:   PetscLogEventBegin(MAT_GetOrdering,mat,0,0,0);
235:   (*r)(mat,type,rperm,cperm);
236:   ISSetPermutation(*rperm);
237:   ISSetPermutation(*cperm);
238:   /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
239:   ISGetLocalSize(*rperm,&mis);
240:   if (mmat > mis) {MatInodeAdjustForInodes(mat,rperm,cperm);}
241:   PetscLogEventEnd(MAT_GetOrdering,mat,0,0,0);

243:   PetscOptionsGetBool(((PetscObject)mat)->prefix,"-mat_view_ordering_draw",&flg,PETSC_NULL);
244:   if (flg) {
245:     Mat tmat;
246:     flg  = PETSC_FALSE;
247:     PetscOptionsGetBool(((PetscObject)mat)->prefix,"-mat_view_contour",&flg,PETSC_NULL);
248:     if (flg) {
249:       PetscViewerPushFormat(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm),PETSC_VIEWER_DRAW_CONTOUR);
250:     }
251:     MatPermute(mat,*rperm,*cperm,&tmat);
252:     MatView(tmat,PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
253:     if (flg) {
254:       PetscViewerPopFormat(PETSC_VIEWER_DRAW_(((PetscObject)mat)->comm));
255:     }
256:     MatDestroy(&tmat);
257:   }

259:   return(0);
260: }

264: PetscErrorCode MatGetOrderingList(PetscFList *list)
265: {
267:   *list = MatOrderingList;
268:   return(0);
269: }