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