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