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: {
17: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot do default flow ordering for matrix type");
18: }
20: PETSC_INTERN PetscErrorCode MatGetOrdering_Natural(Mat mat,MatOrderingType type,IS *irow,IS *icol)
21: {
23: PetscInt n,i,*ii;
24: PetscBool done;
25: MPI_Comm comm;
28: PetscObjectGetComm((PetscObject)mat,&comm);
29: MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,NULL,NULL,&done);
30: MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,NULL,NULL,&done);
31: if (done) { /* matrix may be "compressed" in symbolic factorization, due to i-nodes or block storage */
32: /*
33: We actually create general index sets because this avoids mallocs to
34: to obtain the indices in the MatSolve() routines.
35: ISCreateStride(PETSC_COMM_SELF,n,0,1,irow);
36: ISCreateStride(PETSC_COMM_SELF,n,0,1,icol);
37: */
38: PetscMalloc1(n,&ii);
39: for (i=0; i<n; i++) ii[i] = i;
40: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_COPY_VALUES,irow);
41: ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,icol);
42: } else {
43: PetscInt start,end;
45: MatGetOwnershipRange(mat,&start,&end);
46: ISCreateStride(comm,end-start,start,1,irow);
47: ISCreateStride(comm,end-start,start,1,icol);
48: }
49: ISSetIdentity(*irow);
50: ISSetIdentity(*icol);
51: return(0);
52: }
54: /*
55: Orders the rows (and columns) by the lengths of the rows.
56: This produces a symmetric Ordering but does not require a
57: matrix with symmetric non-zero structure.
58: */
59: PETSC_INTERN PetscErrorCode MatGetOrdering_RowLength(Mat mat,MatOrderingType type,IS *irow,IS *icol)
60: {
62: PetscInt n,*permr,*lens,i;
63: const PetscInt *ia,*ja;
64: PetscBool done;
67: MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,&ia,&ja,&done);
68: if (!done) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get rows for matrix");
70: PetscMalloc2(n,&lens,n,&permr);
71: for (i=0; i<n; i++) {
72: lens[i] = ia[i+1] - ia[i];
73: permr[i] = i;
74: }
75: MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,&ia,&ja,&done);
77: PetscSortIntWithPermutation(n,lens,permr);
79: ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,irow);
80: ISCreateGeneral(PETSC_COMM_SELF,n,permr,PETSC_COPY_VALUES,icol);
81: PetscFree2(lens,permr);
82: return(0);
83: }
85: /*@C
86: MatOrderingRegister - Adds a new sparse matrix ordering to the matrix package.
88: Not Collective
90: Input Parameters:
91: + sname - name of ordering (for example MATORDERINGND)
92: - function - function pointer that creates the ordering
94: Level: developer
96: Sample usage:
97: .vb
98: MatOrderingRegister("my_order", MyOrder);
99: .ve
101: Then, your partitioner can be chosen with the procedural interface via
102: $ MatOrderingSetType(part,"my_order)
103: or at runtime via the option
104: $ -pc_factor_mat_ordering_type my_order
106: .seealso: MatOrderingRegisterAll()
107: @*/
108: PetscErrorCode MatOrderingRegister(const char sname[],PetscErrorCode (*function)(Mat,MatOrderingType,IS*,IS*))
109: {
113: MatInitializePackage();
114: PetscFunctionListAdd(&MatOrderingList,sname,function);
115: return(0);
116: }
118: #include <../src/mat/impls/aij/mpi/mpiaij.h>
119: /*@C
120: MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
121: improve numerical stability of LU factorization.
123: Collective on Mat
125: Input Parameters:
126: + mat - the matrix
127: - type - type of reordering, one of the following:
128: $ MATORDERINGNATURAL_OR_ND - Nested dissection unless matrix is SBAIJ then it is natural
129: $ MATORDERINGNATURAL - Natural
130: $ MATORDERINGND - Nested Dissection
131: $ MATORDERING1WD - One-way Dissection
132: $ MATORDERINGRCM - Reverse Cuthill-McKee
133: $ MATORDERINGQMD - Quotient Minimum Degree
134: $ MATORDERINGEXTERNAL - Use an ordering internal to the factorzation package and do not compute or use PETSc's
136: Output Parameters:
137: + rperm - row permutation indices
138: - cperm - column permutation indices
140: Options Database Key:
141: + -mat_view_ordering draw - plots matrix nonzero structure in new ordering
142: - -pc_factor_mat_ordering_type <nd,natural,..> - ordering to use with PCs based on factorization, LU, ILU, Cholesky, ICC
144: Level: intermediate
146: Notes:
147: This DOES NOT actually reorder the matrix; it merely returns two index sets
148: that define a reordering. This is usually not used directly, rather use the
149: options PCFactorSetMatOrderingType()
151: The user can define additional orderings; see MatOrderingRegister().
153: These are generally only implemented for sequential sparse matrices.
155: Some external packages that PETSc can use for direct factorization such as SuperLU do not accept orderings provided by
156: this call.
158: If MATORDERINGEXTERNAL is used then PETSc does not compute an ordering and utilizes one built into the factorization package
160: fill, reordering, natural, Nested Dissection,
161: One-way Dissection, Cholesky, Reverse Cuthill-McKee,
162: Quotient Minimum Degree
164: .seealso: MatOrderingRegister(), PCFactorSetMatOrderingType(), MatColoring, MatColoringCreate()
165: @*/
166: PetscErrorCode MatGetOrdering(Mat mat,MatOrderingType type,IS *rperm,IS *cperm)
167: {
169: PetscInt mmat,nmat,mis;
170: PetscErrorCode (*r)(Mat,MatOrderingType,IS*,IS*);
171: PetscBool flg,ismpiaij;
177: if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
178: if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
179: if (!type) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Ordering type cannot be null");
181: PetscStrcmp(type,MATORDERINGEXTERNAL,&flg);
182: if (flg) {
183: *rperm = NULL;
184: *cperm = NULL;
185: return(0);
186: }
188: /* This code is terrible. MatGetOrdering() multiple dispatch should use matrix and this code should move to impls/aij/mpi. */
189: PetscObjectTypeCompare((PetscObject)mat,MATMPIAIJ,&ismpiaij);
190: if (ismpiaij) { /* Reorder using diagonal block */
191: Mat Ad,Ao;
192: const PetscInt *colmap;
193: IS lrowperm,lcolperm;
194: PetscInt i,rstart,rend,*idx;
195: const PetscInt *lidx;
197: MatMPIAIJGetSeqAIJ(mat,&Ad,&Ao,&colmap);
198: MatGetOrdering(Ad,type,&lrowperm,&lcolperm);
199: MatGetOwnershipRange(mat,&rstart,&rend);
200: /* Remap row index set to global space */
201: ISGetIndices(lrowperm,&lidx);
202: PetscMalloc1(rend-rstart,&idx);
203: for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
204: ISRestoreIndices(lrowperm,&lidx);
205: ISDestroy(&lrowperm);
206: ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,rperm);
207: ISSetPermutation(*rperm);
208: /* Remap column index set to global space */
209: ISGetIndices(lcolperm,&lidx);
210: PetscMalloc1(rend-rstart,&idx);
211: for (i=0; i+rstart<rend; i++) idx[i] = rstart + lidx[i];
212: ISRestoreIndices(lcolperm,&lidx);
213: ISDestroy(&lcolperm);
214: ISCreateGeneral(PetscObjectComm((PetscObject)mat),rend-rstart,idx,PETSC_OWN_POINTER,cperm);
215: ISSetPermutation(*cperm);
216: return(0);
217: }
219: if (!mat->rmap->N) { /* matrix has zero rows */
220: ISCreateStride(PETSC_COMM_SELF,0,0,1,cperm);
221: ISCreateStride(PETSC_COMM_SELF,0,0,1,rperm);
222: ISSetIdentity(*cperm);
223: ISSetIdentity(*rperm);
224: return(0);
225: }
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: MatOrderingRegisterAll();
231: PetscFunctionListFind(MatOrderingList,type,&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: PetscOptionsHasName(((PetscObject)mat)->options,((PetscObject)mat)->prefix,"-mat_view_ordering",&flg);
244: if (flg) {
245: Mat tmat;
246: MatPermute(mat,*rperm,*cperm,&tmat);
247: MatViewFromOptions(tmat,(PetscObject)mat,"-mat_view_ordering");
248: MatDestroy(&tmat);
249: }
250: return(0);
251: }
253: PetscErrorCode MatGetOrderingList(PetscFunctionList *list)
254: {
256: *list = MatOrderingList;
257: return(0);
258: }