Actual source code: sorder.c

  1: /*$Id: sorder.c,v 1.90 2001/06/21 21:17:30 bsmith Exp $*/
  2: /*
  3:      Provides the code that allows PETSc users to register their own
  4:   sequential matrix Ordering routines.
  5: */
 6:  #include src/mat/matimpl.h
 7:  #include petscsys.h

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

 12: EXTERN int MatOrdering_Flow_SeqAIJ(Mat,MatOrderingType,IS *,IS *);

 14: #undef __FUNCT__  
 16: int MatOrdering_Flow(Mat mat,MatOrderingType type,IS *irow,IS *icol)
 17: {
 19:   SETERRQ(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
 26: #undef __FUNCT__  
 28: int MatOrdering_Natural(Mat mat,MatOrderingType type,IS *irow,IS *icol)
 29: {
 30:   int        n,ierr,i,*ii;
 31:   PetscTruth done;
 32:   MPI_Comm   comm;

 35:   PetscObjectGetComm((PetscObject)mat,&comm);
 36:   MatGetRowIJ(mat,0,PETSC_FALSE,&n,PETSC_NULL,PETSC_NULL,&done);
 37:   MatRestoreRowIJ(mat,0,PETSC_FALSE,&n,PETSC_NULL,PETSC_NULL,&done);
 38:   if (done) { /* matrix may be "compressed" in symbolic factorization, due to i-nodes or block storage */
 39:     /*
 40:       We actually create general index sets because this avoids mallocs to
 41:       to obtain the indices in the MatSolve() routines.
 42:       ISCreateStride(PETSC_COMM_SELF,n,0,1,irow);
 43:       ISCreateStride(PETSC_COMM_SELF,n,0,1,icol);
 44:     */
 45:     PetscMalloc(n*sizeof(int),&ii);
 46:     for (i=0; i<n; i++) ii[i] = i;
 47:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,irow);
 48:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,icol);
 49:     PetscFree(ii);
 50:   } else {
 51:     int 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: */
 69: #undef __FUNCT__  
 71: int MatOrdering_RowLength(Mat mat,MatOrderingType type,IS *irow,IS *icol)
 72: {
 73:   int        ierr,n,*ia,*ja,*permr,*lens,i;
 74:   PetscTruth done;

 77:   MatGetRowIJ(mat,0,PETSC_FALSE,&n,&ia,&ja,&done);
 78:   if (!done) SETERRQ(PETSC_ERR_SUP,"Cannot get rows for matrix");

 80:   ierr  = PetscMalloc(2*n*sizeof(int),&lens);
 81:   permr = lens + n;
 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,&n,&ia,&ja,&done);

 88:   PetscSortIntWithPermutation(n,lens,permr);

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

 97: /*MC
 98:    MatOrderingRegisterDynamic - Adds a new sparse matrix ordering to the 
 99:                                matrix package. 

101:    Synopsis:
102:    int MatOrderingRegisterDynamic(char *name_ordering,char *path,char *name_create,int (*routine_create)(MatOrdering))

104:    Not Collective

106:    Input Parameters:
107: +  sname - name of ordering (for example MATORDERING_ND)
108: .  path - location of library where creation routine is 
109: .  name - name of function that creates the ordering type,a string
110: -  function - function pointer that creates the ordering

112:    Level: developer

114:    If dynamic libraries are used, then the fourth input argument (function)
115:    is ignored.

117:    Sample usage:
118: .vb
119:    MatOrderingRegisterDynamic("my_order",/home/username/my_lib/lib/libO/solaris/mylib.a,
120:                "MyOrder",MyOrder);
121: .ve

123:    Then, your partitioner can be chosen with the procedural interface via
124: $     MatOrderingSetType(part,"my_order)
125:    or at runtime via the option
126: $     -pc_ilu_mat_ordering_type my_order
127: $     -pc_lu_mat_ordering_type my_order

129:    ${PETSC_ARCH} and ${BOPT} occuring in pathname will be replaced with appropriate values.

131: .keywords: matrix, ordering, register

133: .seealso: MatOrderingRegisterDestroy(), MatOrderingRegisterAll()
134: M*/

136: #undef __FUNCT__  
138: int MatOrderingRegister(char *sname,char *path,char *name,int (*function)(Mat,MatOrderingType,IS*,IS*))
139: {
140:   int  ierr;
141:   char fullname[256];

144:   PetscFListConcat(path,name,fullname);
145:   PetscFListAdd(&MatOrderingList,sname,fullname,(void (*)(void))function);
146:   return(0);
147: }

149: #undef __FUNCT__  
151: /*@C
152:    MatOrderingRegisterDestroy - Frees the list of ordering routines.

154:    Not collective

156:    Level: developer
157:    
158: .keywords: matrix, register, destroy

160: .seealso: MatOrderingRegisterDynamic(), MatOrderingRegisterAll()
161: @*/
162: int MatOrderingRegisterDestroy(void)
163: {

167:   if (MatOrderingList) {
168:     PetscFListDestroy(&MatOrderingList);
169:     MatOrderingList = 0;
170:   }
171:   return(0);
172: }

174: EXTERN int MatAdjustForInodes(Mat,IS *,IS *);

176:  #include src/mat/impls/aij/mpi/mpiaij.h
177: #undef __FUNCT__  
179: /*@C
180:    MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
181:    improve numerical stability of LU factorization.

183:    Collective on Mat

185:    Input Parameters:
186: +  mat - the matrix
187: -  type - type of reordering, one of the following:
188: $      MATORDERING_NATURAL - Natural
189: $      MATORDERING_ND - Nested Dissection
190: $      MATORDERING_1WD - One-way Dissection
191: $      MATORDERING_RCM - Reverse Cuthill-McKee
192: $      MATORDERING_QMD - Quotient Minimum Degree

194:    Output Parameters:
195: +  rperm - row permutation indices
196: -  cperm - column permutation indices


199:    Options Database Key:
200: . -mat_view_ordering_draw - plots matrix nonzero structure in new ordering

202:    Level: intermediate
203:    
204:    Notes:
205:       This DOES NOT actually reorder the matrix; it merely returns two index sets
206:    that define a reordering. This is usually not used directly, rather use the 
207:    options PCLUSetMatOrdering() or PCILUSetMatOrdering().

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

211: .keywords: matrix, set, ordering, factorization, direct, ILU, LU,
212:            fill, reordering, natural, Nested Dissection,
213:            One-way Dissection, Cholesky, Reverse Cuthill-McKee, 
214:            Quotient Minimum Degree

216: .seealso:   MatOrderingRegisterDynamic(), PCLUSetMatOrdering(), PCILUSetMatOrdering()
217: @*/
218: int MatGetOrdering(Mat mat,MatOrderingType type,IS *rperm,IS *cperm)
219: {
220:   int         ierr,mmat,nmat,mis,m;
221:   int         (*r)(Mat,MatOrderingType,IS*,IS*);
222:   PetscTruth  flg,isseqdense,ismpidense;

226:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
227:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

229:   PetscTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);
230:   PetscTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense);
231:   if (isseqdense || ismpidense) {
232:     MatGetLocalSize(mat,&m,PETSC_NULL);
233:     /*
234:        Dense matrices only give natural ordering
235:     */
236:     ISCreateStride(PETSC_COMM_SELF,0,m,1,cperm);
237:     ISCreateStride(PETSC_COMM_SELF,0,m,1,rperm);
238:     ISSetIdentity(*cperm);
239:     ISSetIdentity(*rperm);
240:     ISSetPermutation(*rperm);
241:     ISSetPermutation(*cperm);
242:     return(0);
243:   }

245:   if (!mat->M) { /* matrix has zero rows */
246:     ISCreateStride(PETSC_COMM_SELF,0,0,1,cperm);
247:     ISCreateStride(PETSC_COMM_SELF,0,0,1,rperm);
248:     ISSetIdentity(*cperm);
249:     ISSetIdentity(*rperm);
250:     ISSetPermutation(*rperm);
251:     ISSetPermutation(*cperm);
252:     return(0);
253:   }

255:   if (!MatOrderingRegisterAllCalled) {
256:     MatOrderingRegisterAll(PETSC_NULL);
257:   }

259:   PetscLogEventBegin(MAT_GetOrdering,mat,0,0,0);
260:    PetscFListFind(mat->comm,MatOrderingList,type,(void (**)(void)) &r);
261:   if (!r) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);}

263:   (*r)(mat,type,rperm,cperm);
264:   ISSetPermutation(*rperm);
265:   ISSetPermutation(*cperm);

267:   /*
268:       Adjust for inode (reduced matrix ordering) only if row permutation
269:     is smaller then matrix size
270:   */
271:   MatGetLocalSize(mat,&mmat,&nmat);
272:   ISGetLocalSize(*rperm,&mis);
273:   if (mmat > mis) {
274:     MatAdjustForInodes(mat,rperm,cperm);
275:   }

277:   PetscLogEventEnd(MAT_GetOrdering,mat,0,0,0);

279:   PetscOptionsHasName(PETSC_NULL,"-mat_view_ordering_draw",&flg);
280:   if (flg) {
281:     Mat tmat;
282:     PetscOptionsHasName(PETSC_NULL,"-mat_view_contour",&flg);
283:     if (flg) {
284:       PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);
285:     }
286:     MatPermute(mat,*rperm,*cperm,&tmat);
287:     MatView(tmat,PETSC_VIEWER_DRAW_(mat->comm));
288:     if (flg) {
289:       PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));
290:     }
291:     MatDestroy(tmat);
292:   }

294:   return(0);
295: }

297: #undef __FUNCT__  
299: int MatGetOrderingList(PetscFList *list)
300: {
302:   *list = MatOrderingList;
303:   return(0);
304: }