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