Actual source code: sorder.c
petsc-3.6.1 2015-08-06
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: PetscFunctionList MatOrderingList = 0;
10: PetscBool MatOrderingRegisterAllCalled = PETSC_FALSE;
12: extern PetscErrorCode MatGetOrdering_Flow_SeqAIJ(Mat,MatOrderingType,IS*,IS*);
16: PetscErrorCode MatGetOrdering_Flow(Mat mat,MatOrderingType type,IS *irow,IS *icol)
17: {
19: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot do default flow ordering for matrix type");
20: #if !defined(PETSC_USE_DEBUG)
21: return(0);
22: #endif
23: }
27: PETSC_EXTERN PetscErrorCode MatGetOrdering_Natural(Mat mat,MatOrderingType type,IS *irow,IS *icol)
28: {
30: PetscInt n,i,*ii;
31: PetscBool done;
32: MPI_Comm comm;
35: PetscObjectGetComm((PetscObject)mat,&comm);
36: MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,NULL,NULL,&done);
37: MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,NULL,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: PetscMalloc1(n,&ii);
46: for (i=0; i<n; i++) ii[i] = i;
47: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_COPY_VALUES,irow);
48: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,icol);
49: } else {
50: PetscInt start,end;
52: MatGetOwnershipRange(mat,&start,&end);
53: ISCreateStride(comm,end-start,start,1,irow);
54: ISCreateStride(comm,end-start,start,1,icol);
55: }
56: ISSetIdentity(*irow);
57: ISSetIdentity(*icol);
58: return(0);
59: }
61: /*
62: Orders the rows (and columns) by the lengths of the rows.
63: This produces a symmetric Ordering but does not require a
64: matrix with symmetric non-zero structure.
65: */
68: PETSC_EXTERN PetscErrorCode MatGetOrdering_RowLength(Mat mat,MatOrderingType type,IS *irow,IS *icol)
69: {
71: PetscInt n,*permr,*lens,i;
72: const PetscInt *ia,*ja;
73: PetscBool done;
76: MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,&ia,&ja,&done);
77: if (!done) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get rows for matrix");
79: PetscMalloc2(n,&lens,n,&permr);
80: for (i=0; i<n; i++) {
81: lens[i] = ia[i+1] - ia[i];
82: permr[i] = i;
83: }
84: MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,&ia,&ja,&done);
86: PetscSortIntWithPermutation(n,lens,permr);
88: ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,irow);
89: ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,icol);
90: PetscFree2(lens,permr);
91: return(0);
92: }
96: /*@C
97: MatOrderingRegister - Adds a new sparse matrix ordering to the matrix package.
99: Not Collective
101: Input Parameters:
102: + sname - name of ordering (for example MATORDERINGND)
103: - function - function pointer that creates the ordering
105: Level: developer
107: Sample usage:
108: .vb
109: MatOrderingRegister("my_order", MyOrder);
110: .ve
112: Then, your partitioner can be chosen with the procedural interface via
113: $ MatOrderingSetType(part,"my_order)
114: or at runtime via the option
115: $ -pc_factor_mat_ordering_type my_order
117: .keywords: matrix, ordering, register
119: .seealso: MatOrderingRegisterDestroy(), MatOrderingRegisterAll()
120: @*/
121: PetscErrorCode MatOrderingRegister(const char sname[],PetscErrorCode (*function)(Mat,MatOrderingType,IS*,IS*))
122: {
126: PetscFunctionListAdd(&MatOrderingList,sname,function);
127: return(0);
128: }
130: #include <../src/mat/impls/aij/mpi/mpiaij.h>
133: /*@C
134: MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
135: improve numerical stability of LU factorization.
137: Collective on Mat
139: Input Parameters:
140: + mat - the matrix
141: - type - type of reordering, one of the following:
142: $ MATORDERINGNATURAL - Natural
143: $ MATORDERINGND - Nested Dissection
144: $ MATORDERING1WD - One-way Dissection
145: $ MATORDERINGRCM - Reverse Cuthill-McKee
146: $ MATORDERINGQMD - Quotient Minimum Degree
148: Output Parameters:
149: + rperm - row permutation indices
150: - cperm - column permutation indices
153: Options Database Key:
154: . -mat_view_ordering draw - plots matrix nonzero structure in new ordering
156: Level: intermediate
158: Notes:
159: This DOES NOT actually reorder the matrix; it merely returns two index sets
160: that define a reordering. This is usually not used directly, rather use the
161: options PCFactorSetMatOrderingType()
163: The user can define additional orderings; see MatOrderingRegister().
165: These are generally only implemented for sequential sparse matrices.
167: The external packages that PETSc can use for direct factorization such as SuperLU do not accept orderings provided by
168: this call.
171: .keywords: matrix, set, ordering, factorization, direct, ILU, LU,
172: fill, reordering, natural, Nested Dissection,
173: One-way Dissection, Cholesky, Reverse Cuthill-McKee,
174: Quotient Minimum Degree
176: .seealso: MatOrderingRegister(), PCFactorSetMatOrderingType()
177: @*/
178: PetscErrorCode MatGetOrdering(Mat mat,MatOrderingType type,IS *rperm,IS *cperm)
179: {
181: PetscInt mmat,nmat,mis,m;
182: PetscErrorCode (*r)(Mat,MatOrderingType,IS*,IS*);
183: PetscBool flg = PETSC_FALSE,isseqdense,ismpidense,ismpiaij,ismpibaij,ismpisbaij,ismpiaijcusp,ismpiaijcusparse,iselemental;
189: if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
190: if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
192: /* This code is terrible. MatGetOrdering() multiple dispatch should use matrix and this code should move to impls/aij/mpi. */
193: PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJ,&ismpiaij);
194: if (ismpiaij) { /* Reorder using diagonal block */
195: Mat Ad,Ao;
196: const PetscInt *colmap;
197: IS lrowperm,lcolperm;
198: PetscInt i,rstart,rend,*idx;
199: const PetscInt *lidx;
201: MatMPIAIJGetSeqAIJ(mat,&Ad,&Ao,&colmap);
202: MatGetOrdering(Ad,type,&lrowperm,&lcolperm);
203: MatGetOwnershipRange(mat,&rstart,&rend);
204: /* Remap row index set to global space */
205: ISGetIndices(lrowperm,&lidx);
206: PetscMalloc1(rend-rstart,&idx);
207: for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
208: ISRestoreIndices(lrowperm,&lidx);
209: ISDestroy(&lrowperm);
210: ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,rperm);
211: ISSetPermutation(*rperm);
212: /* Remap column index set to global space */
213: ISGetIndices(lcolperm,&lidx);
214: PetscMalloc1(rend-rstart,&idx);
215: for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
216: ISRestoreIndices(lcolperm,&lidx);
217: ISDestroy(&lcolperm);
218: ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,cperm);
219: ISSetPermutation(*cperm);
220: return(0);
221: }
223: /* this chunk of code is REALLY bad, should maybe get the ordering from the factor matrix,
224: then those that don't support orderings will handle their cases themselves. */
225: PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);
226: PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense);
227: PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJCUSP,&ismpiaijcusp);
228: PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJCUSPARSE,&ismpiaijcusparse);
229: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&ismpibaij);
230: PetscObjectTypeCompare((PetscObject)mat,MATMPISBAIJ,&ismpisbaij);
231: PetscObjectTypeCompare((PetscObject)mat,MATELEMENTAL,&iselemental);
232: if (isseqdense || ismpidense || ismpibaij || ismpisbaij || ismpiaijcusp || ismpiaijcusparse || iselemental) {
233: MatGetLocalSize(mat,&m,NULL);
234: /*
235: These matrices only give natural ordering
236: */
237: ISCreateStride(PETSC_COMM_SELF,m,0,1,cperm);
238: ISCreateStride(PETSC_COMM_SELF,m,0,1,rperm);
239: ISSetIdentity(*cperm);
240: ISSetIdentity(*rperm);
241: return(0);
242: }
244: if (!mat->rmap->N) { /* matrix has zero rows */
245: ISCreateStride(PETSC_COMM_SELF,0,0,1,cperm);
246: ISCreateStride(PETSC_COMM_SELF,0,0,1,rperm);
247: ISSetIdentity(*cperm);
248: ISSetIdentity(*rperm);
249: return(0);
250: }
252: MatGetLocalSize(mat,&mmat,&nmat);
253: if (mmat != nmat) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",mmat,nmat);
255: MatOrderingRegisterAll();
256: PetscFunctionListFind(MatOrderingList,type,&r);
257: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);
259: PetscLogEventBegin(MAT_GetOrdering,mat,0,0,0);
260: (*r)(mat,type,rperm,cperm);
261: ISSetPermutation(*rperm);
262: ISSetPermutation(*cperm);
263: /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
264: ISGetLocalSize(*rperm,&mis);
265: if (mmat > mis) {MatInodeAdjustForInodes(mat,rperm,cperm);}
266: PetscLogEventEnd(MAT_GetOrdering,mat,0,0,0);
269: PetscOptionsHasName(((PetscObject)mat)->prefix,"-mat_view_ordering",&flg);
270: if (flg) {
271: Mat tmat;
272: MatPermute(mat,*rperm,*cperm,&tmat);
273: MatViewFromOptions(tmat,(PetscObject)mat,"-mat_view_ordering");
274: MatDestroy(&tmat);
275: }
276: return(0);
277: }
281: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
282: {
284: *list = MatOrderingList;
285: return(0);
286: }