Actual source code: sorder.c
petsc-3.8.4 2018-03-24
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 = 0;
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: {
17: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot do default flow ordering for matrix type");
18: #if !defined(PETSC_USE_DEBUG)
19: return(0);
20: #endif
21: }
23: PETSC_INTERN PetscErrorCode MatGetOrdering_Natural(Mat mat,MatOrderingType type,IS *irow,IS *icol)
24: {
26: PetscInt n,i,*ii;
27: PetscBool done;
28: MPI_Comm comm;
31: PetscObjectGetComm((PetscObject)mat,&comm);
32: MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,NULL,NULL,&done);
33: MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,NULL,NULL,&done);
34: if (done) { /* matrix may be "compressed" in symbolic factorization, due to i-nodes or block storage */
35: /*
36: We actually create general index sets because this avoids mallocs to
37: to obtain the indices in the MatSolve() routines.
38: ISCreateStride(PETSC_COMM_SELF,n,0,1,irow);
39: ISCreateStride(PETSC_COMM_SELF,n,0,1,icol);
40: */
41: PetscMalloc1(n,&ii);
42: for (i=0; i<n; i++) ii[i] = i;
43: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_COPY_VALUES,irow);
44: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,icol);
45: } else {
46: PetscInt start,end;
48: MatGetOwnershipRange(mat,&start,&end);
49: ISCreateStride(comm,end-start,start,1,irow);
50: ISCreateStride(comm,end-start,start,1,icol);
51: }
52: ISSetIdentity(*irow);
53: ISSetIdentity(*icol);
54: return(0);
55: }
57: /*
58: Orders the rows (and columns) by the lengths of the rows.
59: This produces a symmetric Ordering but does not require a
60: matrix with symmetric non-zero structure.
61: */
62: PETSC_INTERN PetscErrorCode MatGetOrdering_RowLength(Mat mat,MatOrderingType type,IS *irow,IS *icol)
63: {
65: PetscInt n,*permr,*lens,i;
66: const PetscInt *ia,*ja;
67: PetscBool done;
70: MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,&ia,&ja,&done);
71: if (!done) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get rows for matrix");
73: PetscMalloc2(n,&lens,n,&permr);
74: for (i=0; i<n; i++) {
75: lens[i] = ia[i+1] - ia[i];
76: permr[i] = i;
77: }
78: MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,&ia,&ja,&done);
80: PetscSortIntWithPermutation(n,lens,permr);
82: ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,irow);
83: ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,icol);
84: PetscFree2(lens,permr);
85: return(0);
86: }
88: /*@C
89: MatOrderingRegister - Adds a new sparse matrix ordering to the matrix package.
91: Not Collective
93: Input Parameters:
94: + sname - name of ordering (for example MATORDERINGND)
95: - function - function pointer that creates the ordering
97: Level: developer
99: Sample usage:
100: .vb
101: MatOrderingRegister("my_order", MyOrder);
102: .ve
104: Then, your partitioner can be chosen with the procedural interface via
105: $ MatOrderingSetType(part,"my_order)
106: or at runtime via the option
107: $ -pc_factor_mat_ordering_type my_order
109: .keywords: matrix, ordering, register
111: .seealso: MatOrderingRegisterDestroy(), MatOrderingRegisterAll()
112: @*/
113: PetscErrorCode MatOrderingRegister(const char sname[],PetscErrorCode (*function)(Mat,MatOrderingType,IS*,IS*))
114: {
118: PetscFunctionListAdd(&MatOrderingList,sname,function);
119: return(0);
120: }
122: #include <../src/mat/impls/aij/mpi/mpiaij.h>
123: /*@C
124: MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
125: improve numerical stability of LU factorization.
127: Collective on Mat
129: Input Parameters:
130: + mat - the matrix
131: - type - type of reordering, one of the following:
132: $ MATORDERINGNATURAL - Natural
133: $ MATORDERINGND - Nested Dissection
134: $ MATORDERING1WD - One-way Dissection
135: $ MATORDERINGRCM - Reverse Cuthill-McKee
136: $ MATORDERINGQMD - Quotient Minimum Degree
138: Output Parameters:
139: + rperm - row permutation indices
140: - cperm - column permutation indices
143: Options Database Key:
144: . -mat_view_ordering draw - plots matrix nonzero structure in new ordering
146: Level: intermediate
148: Notes:
149: This DOES NOT actually reorder the matrix; it merely returns two index sets
150: that define a reordering. This is usually not used directly, rather use the
151: options PCFactorSetMatOrderingType()
153: The user can define additional orderings; see MatOrderingRegister().
155: These are generally only implemented for sequential sparse matrices.
157: The external packages that PETSc can use for direct factorization such as SuperLU do not accept orderings provided by
158: this call.
161: .keywords: matrix, set, ordering, factorization, direct, ILU, LU,
162: fill, reordering, natural, Nested Dissection,
163: One-way Dissection, Cholesky, Reverse Cuthill-McKee,
164: Quotient Minimum Degree
166: .seealso: MatOrderingRegister(), PCFactorSetMatOrderingType()
167: @*/
168: PetscErrorCode MatGetOrdering(Mat mat,MatOrderingType type,IS *rperm,IS *cperm)
169: {
171: PetscInt mmat,nmat,mis,m;
172: PetscErrorCode (*r)(Mat,MatOrderingType,IS*,IS*);
173: PetscBool flg = PETSC_FALSE,isseqdense,ismpidense,ismpiaij,ismpibaij,ismpisbaij,ismpiaijcusp,ismpiaijcusparse,iselemental;
179: if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
180: if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
182: /* This code is terrible. MatGetOrdering() multiple dispatch should use matrix and this code should move to impls/aij/mpi. */
183: PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJ,&ismpiaij);
184: if (ismpiaij) { /* Reorder using diagonal block */
185: Mat Ad,Ao;
186: const PetscInt *colmap;
187: IS lrowperm,lcolperm;
188: PetscInt i,rstart,rend,*idx;
189: const PetscInt *lidx;
191: MatMPIAIJGetSeqAIJ(mat,&Ad,&Ao,&colmap);
192: MatGetOrdering(Ad,type,&lrowperm,&lcolperm);
193: MatGetOwnershipRange(mat,&rstart,&rend);
194: /* Remap row index set to global space */
195: ISGetIndices(lrowperm,&lidx);
196: PetscMalloc1(rend-rstart,&idx);
197: for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
198: ISRestoreIndices(lrowperm,&lidx);
199: ISDestroy(&lrowperm);
200: ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,rperm);
201: ISSetPermutation(*rperm);
202: /* Remap column index set to global space */
203: ISGetIndices(lcolperm,&lidx);
204: PetscMalloc1(rend-rstart,&idx);
205: for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
206: ISRestoreIndices(lcolperm,&lidx);
207: ISDestroy(&lcolperm);
208: ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,cperm);
209: ISSetPermutation(*cperm);
210: return(0);
211: }
213: /* this chunk of code is REALLY bad, should maybe get the ordering from the factor matrix,
214: then those that don't support orderings will handle their cases themselves. */
215: PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);
216: PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense);
217: PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJCUSP,&ismpiaijcusp);
218: PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJCUSPARSE,&ismpiaijcusparse);
219: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&ismpibaij);
220: PetscObjectTypeCompare((PetscObject)mat,MATMPISBAIJ,&ismpisbaij);
221: PetscObjectTypeCompare((PetscObject)mat,MATELEMENTAL,&iselemental);
222: if (isseqdense || ismpidense || ismpibaij || ismpisbaij || ismpiaijcusp || ismpiaijcusparse || iselemental) {
223: MatGetLocalSize(mat,&m,NULL);
224: /*
225: These matrices only give natural ordering
226: */
227: ISCreateStride(PETSC_COMM_SELF,m,0,1,cperm);
228: ISCreateStride(PETSC_COMM_SELF,m,0,1,rperm);
229: ISSetIdentity(*cperm);
230: ISSetIdentity(*rperm);
231: return(0);
232: }
234: if (!mat->rmap->N) { /* matrix has zero rows */
235: ISCreateStride(PETSC_COMM_SELF,0,0,1,cperm);
236: ISCreateStride(PETSC_COMM_SELF,0,0,1,rperm);
237: ISSetIdentity(*cperm);
238: ISSetIdentity(*rperm);
239: return(0);
240: }
242: MatGetLocalSize(mat,&mmat,&nmat);
243: if (mmat != nmat) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %D columns %D",mmat,nmat);
245: MatOrderingRegisterAll();
246: PetscFunctionListFind(MatOrderingList,type,&r);
247: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);
249: PetscLogEventBegin(MAT_GetOrdering,mat,0,0,0);
250: (*r)(mat,type,rperm,cperm);
251: ISSetPermutation(*rperm);
252: ISSetPermutation(*cperm);
253: /* Adjust for inode (reduced matrix ordering) only if row permutation is smaller the matrix size */
254: ISGetLocalSize(*rperm,&mis);
255: if (mmat > mis) {MatInodeAdjustForInodes(mat,rperm,cperm);}
256: PetscLogEventEnd(MAT_GetOrdering,mat,0,0,0);
259: PetscOptionsHasName(((PetscObject)mat)->options,((PetscObject)mat)->prefix,"-mat_view_ordering",&flg);
260: if (flg) {
261: Mat tmat;
262: MatPermute(mat,*rperm,*cperm,&tmat);
263: MatViewFromOptions(tmat,(PetscObject)mat,"-mat_view_ordering");
264: MatDestroy(&tmat);
265: }
266: return(0);
267: }
269: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
270: {
272: *list = MatOrderingList;
273: return(0);
274: }