Actual source code: sorder.c
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 = NULL;
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: {
16: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot do default flow ordering for matrix type");
17: }
19: PETSC_INTERN PetscErrorCode MatGetOrdering_Natural(Mat mat,MatOrderingType type,IS *irow,IS *icol)
20: {
21: PetscInt n,i,*ii;
22: PetscBool done;
23: MPI_Comm comm;
25: PetscObjectGetComm((PetscObject)mat,&comm);
26: MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,NULL,NULL,&done);
27: MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,NULL,NULL,&done);
28: if (done) { /* matrix may be "compressed" in symbolic factorization, due to i-nodes or block storage */
29: /*
30: We actually create general index sets because this avoids mallocs to
31: to obtain the indices in the MatSolve() routines.
32: ISCreateStride(PETSC_COMM_SELF,n,0,1,irow);
33: ISCreateStride(PETSC_COMM_SELF,n,0,1,icol);
34: */
35: PetscMalloc1(n,&ii);
36: for (i=0; i<n; i++) ii[i] = i;
37: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_COPY_VALUES,irow);
38: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,icol);
39: } else {
40: PetscInt start,end;
42: MatGetOwnershipRange(mat,&start,&end);
43: ISCreateStride(comm,end-start,start,1,irow);
44: ISCreateStride(comm,end-start,start,1,icol);
45: }
46: ISSetIdentity(*irow);
47: ISSetIdentity(*icol);
48: return 0;
49: }
51: /*
52: Orders the rows (and columns) by the lengths of the rows.
53: This produces a symmetric Ordering but does not require a
54: matrix with symmetric non-zero structure.
55: */
56: PETSC_INTERN PetscErrorCode MatGetOrdering_RowLength(Mat mat,MatOrderingType type,IS *irow,IS *icol)
57: {
58: PetscInt n,*permr,*lens,i;
59: const PetscInt *ia,*ja;
60: PetscBool done;
62: MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,&ia,&ja,&done);
65: PetscMalloc2(n,&lens,n,&permr);
66: for (i=0; i<n; i++) {
67: lens[i] = ia[i+1] - ia[i];
68: permr[i] = i;
69: }
70: MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,&ia,&ja,&done);
72: PetscSortIntWithPermutation(n,lens,permr);
74: ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,irow);
75: ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,icol);
76: PetscFree2(lens,permr);
77: return 0;
78: }
80: /*@C
81: MatOrderingRegister - Adds a new sparse matrix ordering to the matrix package.
83: Not Collective
85: Input Parameters:
86: + sname - name of ordering (for example MATORDERINGND)
87: - function - function pointer that creates the ordering
89: Level: developer
91: Sample usage:
92: .vb
93: MatOrderingRegister("my_order", MyOrder);
94: .ve
96: Then, your partitioner can be chosen with the procedural interface via
97: $ MatOrderingSetType(part,"my_order)
98: or at runtime via the option
99: $ -pc_factor_mat_ordering_type my_order
101: .seealso: MatOrderingRegisterAll()
102: @*/
103: PetscErrorCode MatOrderingRegister(const char sname[],PetscErrorCode (*function)(Mat,MatOrderingType,IS*,IS*))
104: {
105: MatInitializePackage();
106: PetscFunctionListAdd(&MatOrderingList,sname,function);
107: return 0;
108: }
110: #include <../src/mat/impls/aij/mpi/mpiaij.h>
111: /*@C
112: MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
113: improve numerical stability of LU factorization.
115: Collective on Mat
117: Input Parameters:
118: + mat - the matrix
119: - type - type of reordering, one of the following:
120: $ MATORDERINGNATURAL_OR_ND - Nested dissection unless matrix is SBAIJ then it is natural
121: $ MATORDERINGNATURAL - Natural
122: $ MATORDERINGND - Nested Dissection
123: $ MATORDERING1WD - One-way Dissection
124: $ MATORDERINGRCM - Reverse Cuthill-McKee
125: $ MATORDERINGQMD - Quotient Minimum Degree
126: $ MATORDERINGEXTERNAL - Use an ordering internal to the factorzation package and do not compute or use PETSc's
128: Output Parameters:
129: + rperm - row permutation indices
130: - cperm - column permutation indices
132: Options Database Key:
133: + -mat_view_ordering draw - plots matrix nonzero structure in new ordering
134: - -pc_factor_mat_ordering_type <nd,natural,..> - ordering to use with PCs based on factorization, LU, ILU, Cholesky, ICC
136: Level: intermediate
138: Notes:
139: This DOES NOT actually reorder the matrix; it merely returns two index sets
140: that define a reordering. This is usually not used directly, rather use the
141: options PCFactorSetMatOrderingType()
143: The user can define additional orderings; see MatOrderingRegister().
145: These are generally only implemented for sequential sparse matrices.
147: Some external packages that PETSc can use for direct factorization such as SuperLU do not accept orderings provided by
148: this call.
150: If MATORDERINGEXTERNAL is used then PETSc does not compute an ordering and utilizes one built into the factorization package
152: fill, reordering, natural, Nested Dissection,
153: One-way Dissection, Cholesky, Reverse Cuthill-McKee,
154: Quotient Minimum Degree
156: .seealso: MatOrderingRegister(), PCFactorSetMatOrderingType(), MatColoring, MatColoringCreate()
157: @*/
158: PetscErrorCode MatGetOrdering(Mat mat,MatOrderingType type,IS *rperm,IS *cperm)
159: {
160: PetscInt mmat,nmat,mis;
161: PetscErrorCode (*r)(Mat,MatOrderingType,IS*,IS*);
162: PetscBool flg,ismpiaij;
171: PetscStrcmp(type,MATORDERINGEXTERNAL,&flg);
172: if (flg) {
173: *rperm = NULL;
174: *cperm = NULL;
175: return 0;
176: }
178: /* This code is terrible. MatGetOrdering() multiple dispatch should use matrix and this code should move to impls/aij/mpi. */
179: PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJ,&ismpiaij);
180: if (ismpiaij) { /* Reorder using diagonal block */
181: Mat Ad,Ao;
182: const PetscInt *colmap;
183: IS lrowperm,lcolperm;
184: PetscInt i,rstart,rend,*idx;
185: const PetscInt *lidx;
187: MatMPIAIJGetSeqAIJ(mat,&Ad,&Ao,&colmap);
188: MatGetOrdering(Ad,type,&lrowperm,&lcolperm);
189: MatGetOwnershipRange(mat,&rstart,&rend);
190: /* Remap row index set to global space */
191: ISGetIndices(lrowperm,&lidx);
192: PetscMalloc1(rend-rstart,&idx);
193: for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
194: ISRestoreIndices(lrowperm,&lidx);
195: ISDestroy(&lrowperm);
196: ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,rperm);
197: ISSetPermutation(*rperm);
198: /* Remap column index set to global space */
199: ISGetIndices(lcolperm,&lidx);
200: PetscMalloc1(rend-rstart,&idx);
201: for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
202: ISRestoreIndices(lcolperm,&lidx);
203: ISDestroy(&lcolperm);
204: ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,cperm);
205: ISSetPermutation(*cperm);
206: return 0;
207: }
209: if (!mat->rmap->N) { /* matrix has zero rows */
210: ISCreateStride(PETSC_COMM_SELF,0,0,1,cperm);
211: ISCreateStride(PETSC_COMM_SELF,0,0,1,rperm);
212: ISSetIdentity(*cperm);
213: ISSetIdentity(*rperm);
214: return 0;
215: }
217: MatGetLocalSize(mat,&mmat,&nmat);
220: MatOrderingRegisterAll();
221: PetscFunctionListFind(MatOrderingList,type,&r);
224: PetscLogEventBegin(MAT_GetOrdering,mat,0,0,0);
225: (*r)(mat,type,rperm,cperm);
226: ISSetPermutation(*rperm);
227: ISSetPermutation(*cperm);
228: /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
229: ISGetLocalSize(*rperm,&mis);
230: if (mmat > mis) MatInodeAdjustForInodes(mat,rperm,cperm);
231: PetscLogEventEnd(MAT_GetOrdering,mat,0,0,0);
233: PetscOptionsHasName(((PetscObject)mat)->options,((PetscObject)mat)->prefix,"-mat_view_ordering",&flg);
234: if (flg) {
235: Mat tmat;
236: MatPermute(mat,*rperm,*cperm,&tmat);
237: MatViewFromOptions(tmat,(PetscObject)mat,"-mat_view_ordering");
238: MatDestroy(&tmat);
239: }
240: return 0;
241: }
243: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
244: {
245: *list = MatOrderingList;
246: return 0;
247: }