1 2 /* 3 Defines the basic matrix operations for the ADJ adjacency list matrix data-structure. 4 */ 5 #include <../src/mat/impls/adj/mpi/mpiadj.h> /*I "petscmat.h" I*/ 6 #include <petscsf.h> 7 8 /* 9 The interface should be easy to use for both MatCreateSubMatrix (parallel sub-matrix) and MatCreateSubMatrices (sequential sub-matrices) 10 */ 11 static PetscErrorCode MatCreateSubMatrix_MPIAdj_data(Mat adj, IS irows, IS icols, PetscInt **sadj_xadj, PetscInt **sadj_adjncy, PetscInt **sadj_values) 12 { 13 PetscInt nlrows_is, icols_n, i, j, nroots, nleaves, rlocalindex, *ncols_send, *ncols_recv; 14 PetscInt nlrows_mat, *adjncy_recv, Ncols_recv, Ncols_send, *xadj_recv, *values_recv; 15 PetscInt *ncols_recv_offsets, loc, rnclos, *sadjncy, *sxadj, *svalues; 16 const PetscInt *irows_indices, *icols_indices, *xadj, *adjncy; 17 PetscMPIInt owner; 18 Mat_MPIAdj *a = (Mat_MPIAdj *)adj->data; 19 PetscLayout rmap; 20 MPI_Comm comm; 21 PetscSF sf; 22 PetscSFNode *iremote; 23 PetscBool done; 24 25 PetscFunctionBegin; 26 PetscCall(PetscObjectGetComm((PetscObject)adj, &comm)); 27 PetscCall(MatGetLayouts(adj, &rmap, NULL)); 28 PetscCall(ISGetLocalSize(irows, &nlrows_is)); 29 PetscCall(ISGetIndices(irows, &irows_indices)); 30 PetscCall(PetscMalloc1(nlrows_is, &iremote)); 31 /* construct sf graph*/ 32 nleaves = nlrows_is; 33 for (i = 0; i < nlrows_is; i++) { 34 owner = -1; 35 rlocalindex = -1; 36 PetscCall(PetscLayoutFindOwnerIndex(rmap, irows_indices[i], &owner, &rlocalindex)); 37 iremote[i].rank = owner; 38 iremote[i].index = rlocalindex; 39 } 40 PetscCall(MatGetRowIJ(adj, 0, PETSC_FALSE, PETSC_FALSE, &nlrows_mat, &xadj, &adjncy, &done)); 41 PetscCall(PetscCalloc4(nlrows_mat, &ncols_send, nlrows_is, &xadj_recv, nlrows_is + 1, &ncols_recv_offsets, nlrows_is, &ncols_recv)); 42 nroots = nlrows_mat; 43 for (i = 0; i < nlrows_mat; i++) ncols_send[i] = xadj[i + 1] - xadj[i]; 44 PetscCall(PetscSFCreate(comm, &sf)); 45 PetscCall(PetscSFSetGraph(sf, nroots, nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 46 PetscCall(PetscSFSetType(sf, PETSCSFBASIC)); 47 PetscCall(PetscSFSetFromOptions(sf)); 48 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, ncols_send, ncols_recv, MPI_REPLACE)); 49 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, ncols_send, ncols_recv, MPI_REPLACE)); 50 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, xadj, xadj_recv, MPI_REPLACE)); 51 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, xadj, xadj_recv, MPI_REPLACE)); 52 PetscCall(PetscSFDestroy(&sf)); 53 Ncols_recv = 0; 54 for (i = 0; i < nlrows_is; i++) { 55 Ncols_recv += ncols_recv[i]; 56 ncols_recv_offsets[i + 1] = ncols_recv[i] + ncols_recv_offsets[i]; 57 } 58 Ncols_send = 0; 59 for (i = 0; i < nlrows_mat; i++) Ncols_send += ncols_send[i]; 60 PetscCall(PetscCalloc1(Ncols_recv, &iremote)); 61 PetscCall(PetscCalloc1(Ncols_recv, &adjncy_recv)); 62 nleaves = Ncols_recv; 63 Ncols_recv = 0; 64 for (i = 0; i < nlrows_is; i++) { 65 PetscCall(PetscLayoutFindOwner(rmap, irows_indices[i], &owner)); 66 for (j = 0; j < ncols_recv[i]; j++) { 67 iremote[Ncols_recv].rank = owner; 68 iremote[Ncols_recv++].index = xadj_recv[i] + j; 69 } 70 } 71 PetscCall(ISRestoreIndices(irows, &irows_indices)); 72 /*if we need to deal with edge weights ???*/ 73 if (a->useedgeweights) PetscCall(PetscCalloc1(Ncols_recv, &values_recv)); 74 nroots = Ncols_send; 75 PetscCall(PetscSFCreate(comm, &sf)); 76 PetscCall(PetscSFSetGraph(sf, nroots, nleaves, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 77 PetscCall(PetscSFSetType(sf, PETSCSFBASIC)); 78 PetscCall(PetscSFSetFromOptions(sf)); 79 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, adjncy, adjncy_recv, MPI_REPLACE)); 80 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, adjncy, adjncy_recv, MPI_REPLACE)); 81 if (a->useedgeweights) { 82 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, a->values, values_recv, MPI_REPLACE)); 83 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, a->values, values_recv, MPI_REPLACE)); 84 } 85 PetscCall(PetscSFDestroy(&sf)); 86 PetscCall(MatRestoreRowIJ(adj, 0, PETSC_FALSE, PETSC_FALSE, &nlrows_mat, &xadj, &adjncy, &done)); 87 PetscCall(ISGetLocalSize(icols, &icols_n)); 88 PetscCall(ISGetIndices(icols, &icols_indices)); 89 rnclos = 0; 90 for (i = 0; i < nlrows_is; i++) { 91 for (j = ncols_recv_offsets[i]; j < ncols_recv_offsets[i + 1]; j++) { 92 PetscCall(PetscFindInt(adjncy_recv[j], icols_n, icols_indices, &loc)); 93 if (loc < 0) { 94 adjncy_recv[j] = -1; 95 if (a->useedgeweights) values_recv[j] = -1; 96 ncols_recv[i]--; 97 } else { 98 rnclos++; 99 } 100 } 101 } 102 PetscCall(ISRestoreIndices(icols, &icols_indices)); 103 PetscCall(PetscCalloc1(rnclos, &sadjncy)); 104 if (a->useedgeweights) PetscCall(PetscCalloc1(rnclos, &svalues)); 105 PetscCall(PetscCalloc1(nlrows_is + 1, &sxadj)); 106 rnclos = 0; 107 for (i = 0; i < nlrows_is; i++) { 108 for (j = ncols_recv_offsets[i]; j < ncols_recv_offsets[i + 1]; j++) { 109 if (adjncy_recv[j] < 0) continue; 110 sadjncy[rnclos] = adjncy_recv[j]; 111 if (a->useedgeweights) svalues[rnclos] = values_recv[j]; 112 rnclos++; 113 } 114 } 115 for (i = 0; i < nlrows_is; i++) sxadj[i + 1] = sxadj[i] + ncols_recv[i]; 116 if (sadj_xadj) { 117 *sadj_xadj = sxadj; 118 } else PetscCall(PetscFree(sxadj)); 119 if (sadj_adjncy) { 120 *sadj_adjncy = sadjncy; 121 } else PetscCall(PetscFree(sadjncy)); 122 if (sadj_values) { 123 if (a->useedgeweights) *sadj_values = svalues; 124 else *sadj_values = NULL; 125 } else { 126 if (a->useedgeweights) PetscCall(PetscFree(svalues)); 127 } 128 PetscCall(PetscFree4(ncols_send, xadj_recv, ncols_recv_offsets, ncols_recv)); 129 PetscCall(PetscFree(adjncy_recv)); 130 if (a->useedgeweights) PetscCall(PetscFree(values_recv)); 131 PetscFunctionReturn(PETSC_SUCCESS); 132 } 133 134 static PetscErrorCode MatCreateSubMatrices_MPIAdj_Private(Mat mat, PetscInt n, const IS irow[], const IS icol[], PetscBool subcomm, MatReuse scall, Mat *submat[]) 135 { 136 PetscInt i, irow_n, icol_n, *sxadj, *sadjncy, *svalues; 137 PetscInt *indices, nindx, j, k, loc; 138 PetscMPIInt issame; 139 const PetscInt *irow_indices, *icol_indices; 140 MPI_Comm scomm_row, scomm_col, scomm_mat; 141 142 PetscFunctionBegin; 143 nindx = 0; 144 /* 145 * Estimate a maximum number for allocating memory 146 */ 147 for (i = 0; i < n; i++) { 148 PetscCall(ISGetLocalSize(irow[i], &irow_n)); 149 PetscCall(ISGetLocalSize(icol[i], &icol_n)); 150 nindx = nindx > (irow_n + icol_n) ? nindx : (irow_n + icol_n); 151 } 152 PetscCall(PetscMalloc1(nindx, &indices)); 153 /* construct a submat */ 154 // if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscMalloc1(n,submat)); 155 156 for (i = 0; i < n; i++) { 157 if (subcomm) { 158 PetscCall(PetscObjectGetComm((PetscObject)irow[i], &scomm_row)); 159 PetscCall(PetscObjectGetComm((PetscObject)icol[i], &scomm_col)); 160 PetscCallMPI(MPI_Comm_compare(scomm_row, scomm_col, &issame)); 161 PetscCheck(issame == MPI_IDENT, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "row index set must have the same comm as the col index set"); 162 PetscCallMPI(MPI_Comm_compare(scomm_row, PETSC_COMM_SELF, &issame)); 163 PetscCheck(issame != MPI_IDENT, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, " can not use PETSC_COMM_SELF as comm when extracting a parallel submatrix"); 164 } else { 165 scomm_row = PETSC_COMM_SELF; 166 } 167 /*get sub-matrix data*/ 168 sxadj = NULL; 169 sadjncy = NULL; 170 svalues = NULL; 171 PetscCall(MatCreateSubMatrix_MPIAdj_data(mat, irow[i], icol[i], &sxadj, &sadjncy, &svalues)); 172 PetscCall(ISGetLocalSize(irow[i], &irow_n)); 173 PetscCall(ISGetLocalSize(icol[i], &icol_n)); 174 PetscCall(ISGetIndices(irow[i], &irow_indices)); 175 PetscCall(PetscArraycpy(indices, irow_indices, irow_n)); 176 PetscCall(ISRestoreIndices(irow[i], &irow_indices)); 177 PetscCall(ISGetIndices(icol[i], &icol_indices)); 178 PetscCall(PetscArraycpy(indices + irow_n, icol_indices, icol_n)); 179 PetscCall(ISRestoreIndices(icol[i], &icol_indices)); 180 nindx = irow_n + icol_n; 181 PetscCall(PetscSortRemoveDupsInt(&nindx, indices)); 182 /* renumber columns */ 183 for (j = 0; j < irow_n; j++) { 184 for (k = sxadj[j]; k < sxadj[j + 1]; k++) { 185 PetscCall(PetscFindInt(sadjncy[k], nindx, indices, &loc)); 186 PetscCheck(loc >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "can not find col %" PetscInt_FMT, sadjncy[k]); 187 sadjncy[k] = loc; 188 } 189 } 190 if (scall == MAT_INITIAL_MATRIX) { 191 PetscCall(MatCreateMPIAdj(scomm_row, irow_n, icol_n, sxadj, sadjncy, svalues, submat[i])); 192 } else { 193 Mat sadj = *(submat[i]); 194 Mat_MPIAdj *sa = (Mat_MPIAdj *)((sadj)->data); 195 PetscCall(PetscObjectGetComm((PetscObject)sadj, &scomm_mat)); 196 PetscCallMPI(MPI_Comm_compare(scomm_row, scomm_mat, &issame)); 197 PetscCheck(issame == MPI_IDENT, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "submatrix must have the same comm as the col index set"); 198 PetscCall(PetscArraycpy(sa->i, sxadj, irow_n + 1)); 199 PetscCall(PetscArraycpy(sa->j, sadjncy, sxadj[irow_n])); 200 if (svalues) PetscCall(PetscArraycpy(sa->values, svalues, sxadj[irow_n])); 201 PetscCall(PetscFree(sxadj)); 202 PetscCall(PetscFree(sadjncy)); 203 PetscCall(PetscFree(svalues)); 204 } 205 } 206 PetscCall(PetscFree(indices)); 207 PetscFunctionReturn(PETSC_SUCCESS); 208 } 209 210 static PetscErrorCode MatCreateSubMatricesMPI_MPIAdj(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[]) 211 { 212 /*get sub-matrices across a sub communicator */ 213 PetscFunctionBegin; 214 PetscCall(MatCreateSubMatrices_MPIAdj_Private(mat, n, irow, icol, PETSC_TRUE, scall, submat)); 215 PetscFunctionReturn(PETSC_SUCCESS); 216 } 217 218 static PetscErrorCode MatCreateSubMatrices_MPIAdj(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[]) 219 { 220 PetscFunctionBegin; 221 /*get sub-matrices based on PETSC_COMM_SELF */ 222 PetscCall(MatCreateSubMatrices_MPIAdj_Private(mat, n, irow, icol, PETSC_FALSE, scall, submat)); 223 PetscFunctionReturn(PETSC_SUCCESS); 224 } 225 226 static PetscErrorCode MatView_MPIAdj_ASCII(Mat A, PetscViewer viewer) 227 { 228 Mat_MPIAdj *a = (Mat_MPIAdj *)A->data; 229 PetscInt i, j, m = A->rmap->n; 230 const char *name; 231 PetscViewerFormat format; 232 233 PetscFunctionBegin; 234 PetscCall(PetscObjectGetName((PetscObject)A, &name)); 235 PetscCall(PetscViewerGetFormat(viewer, &format)); 236 if (format == PETSC_VIEWER_ASCII_INFO) { 237 PetscFunctionReturn(PETSC_SUCCESS); 238 } else { 239 PetscCheck(format != PETSC_VIEWER_ASCII_MATLAB, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MATLAB format not supported"); 240 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 241 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 242 for (i = 0; i < m; i++) { 243 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "row %" PetscInt_FMT ":", i + A->rmap->rstart)); 244 for (j = a->i[i]; j < a->i[i + 1]; j++) { 245 if (a->values) { 246 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%" PetscInt_FMT ", %" PetscInt_FMT ") ", a->j[j], a->values[j])); 247 } else { 248 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT " ", a->j[j])); 249 } 250 } 251 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n")); 252 } 253 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 254 PetscCall(PetscViewerFlush(viewer)); 255 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 256 } 257 PetscFunctionReturn(PETSC_SUCCESS); 258 } 259 260 static PetscErrorCode MatView_MPIAdj(Mat A, PetscViewer viewer) 261 { 262 PetscBool iascii; 263 264 PetscFunctionBegin; 265 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 266 if (iascii) PetscCall(MatView_MPIAdj_ASCII(A, viewer)); 267 PetscFunctionReturn(PETSC_SUCCESS); 268 } 269 270 static PetscErrorCode MatDestroy_MPIAdj(Mat mat) 271 { 272 Mat_MPIAdj *a = (Mat_MPIAdj *)mat->data; 273 274 PetscFunctionBegin; 275 #if defined(PETSC_USE_LOG) 276 PetscCall(PetscLogObjectState((PetscObject)mat, "Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT, mat->rmap->n, mat->cmap->n, a->nz)); 277 #endif 278 PetscCall(PetscFree(a->diag)); 279 if (a->freeaij) { 280 if (a->freeaijwithfree) { 281 if (a->i) free(a->i); 282 if (a->j) free(a->j); 283 } else { 284 PetscCall(PetscFree(a->i)); 285 PetscCall(PetscFree(a->j)); 286 PetscCall(PetscFree(a->values)); 287 } 288 } 289 PetscCall(PetscFree(a->rowvalues)); 290 PetscCall(PetscFree(mat->data)); 291 PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL)); 292 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPIAdjSetPreallocation_C", NULL)); 293 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPIAdjCreateNonemptySubcommMat_C", NULL)); 294 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPIAdjToSeq_C", NULL)); 295 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPIAdjToSeqRankZero_C", NULL)); 296 PetscFunctionReturn(PETSC_SUCCESS); 297 } 298 299 static PetscErrorCode MatSetOption_MPIAdj(Mat A, MatOption op, PetscBool flg) 300 { 301 Mat_MPIAdj *a = (Mat_MPIAdj *)A->data; 302 303 PetscFunctionBegin; 304 switch (op) { 305 case MAT_SYMMETRIC: 306 case MAT_STRUCTURALLY_SYMMETRIC: 307 case MAT_HERMITIAN: 308 case MAT_SPD: 309 a->symmetric = flg; 310 break; 311 case MAT_SYMMETRY_ETERNAL: 312 case MAT_STRUCTURAL_SYMMETRY_ETERNAL: 313 case MAT_SPD_ETERNAL: 314 break; 315 default: 316 PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op])); 317 break; 318 } 319 PetscFunctionReturn(PETSC_SUCCESS); 320 } 321 322 static PetscErrorCode MatGetRow_MPIAdj(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v) 323 { 324 Mat_MPIAdj *a = (Mat_MPIAdj *)A->data; 325 326 PetscFunctionBegin; 327 row -= A->rmap->rstart; 328 PetscCheck(row >= 0 && row < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row out of range"); 329 *nz = a->i[row + 1] - a->i[row]; 330 if (v) { 331 PetscInt j; 332 if (a->rowvalues_alloc < *nz) { 333 PetscCall(PetscFree(a->rowvalues)); 334 a->rowvalues_alloc = PetscMax(a->rowvalues_alloc * 2, *nz); 335 PetscCall(PetscMalloc1(a->rowvalues_alloc, &a->rowvalues)); 336 } 337 for (j = 0; j < *nz; j++) a->rowvalues[j] = a->values ? a->values[a->i[row] + j] : 1.0; 338 *v = (*nz) ? a->rowvalues : NULL; 339 } 340 if (idx) *idx = (*nz) ? a->j + a->i[row] : NULL; 341 PetscFunctionReturn(PETSC_SUCCESS); 342 } 343 344 static PetscErrorCode MatRestoreRow_MPIAdj(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v) 345 { 346 PetscFunctionBegin; 347 PetscFunctionReturn(PETSC_SUCCESS); 348 } 349 350 static PetscErrorCode MatEqual_MPIAdj(Mat A, Mat B, PetscBool *flg) 351 { 352 Mat_MPIAdj *a = (Mat_MPIAdj *)A->data, *b = (Mat_MPIAdj *)B->data; 353 PetscBool flag; 354 355 PetscFunctionBegin; 356 /* If the matrix dimensions are not equal,or no of nonzeros */ 357 if ((A->rmap->n != B->rmap->n) || (a->nz != b->nz)) flag = PETSC_FALSE; 358 359 /* if the a->i are the same */ 360 PetscCall(PetscArraycmp(a->i, b->i, A->rmap->n + 1, &flag)); 361 362 /* if a->j are the same */ 363 PetscCall(PetscMemcmp(a->j, b->j, (a->nz) * sizeof(PetscInt), &flag)); 364 365 PetscCall(MPIU_Allreduce(&flag, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 366 PetscFunctionReturn(PETSC_SUCCESS); 367 } 368 369 static PetscErrorCode MatGetRowIJ_MPIAdj(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *m, const PetscInt *inia[], const PetscInt *inja[], PetscBool *done) 370 { 371 PetscInt i; 372 Mat_MPIAdj *a = (Mat_MPIAdj *)A->data; 373 PetscInt **ia = (PetscInt **)inia, **ja = (PetscInt **)inja; 374 375 PetscFunctionBegin; 376 *m = A->rmap->n; 377 *ia = a->i; 378 *ja = a->j; 379 *done = PETSC_TRUE; 380 if (oshift) { 381 for (i = 0; i < (*ia)[*m]; i++) (*ja)[i]++; 382 for (i = 0; i <= (*m); i++) (*ia)[i]++; 383 } 384 PetscFunctionReturn(PETSC_SUCCESS); 385 } 386 387 static PetscErrorCode MatRestoreRowIJ_MPIAdj(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *m, const PetscInt *inia[], const PetscInt *inja[], PetscBool *done) 388 { 389 PetscInt i; 390 Mat_MPIAdj *a = (Mat_MPIAdj *)A->data; 391 PetscInt **ia = (PetscInt **)inia, **ja = (PetscInt **)inja; 392 393 PetscFunctionBegin; 394 PetscCheck(!ia || a->i == *ia, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "ia passed back is not one obtained with MatGetRowIJ()"); 395 PetscCheck(!ja || a->j == *ja, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "ja passed back is not one obtained with MatGetRowIJ()"); 396 if (oshift) { 397 PetscCheck(ia, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "If oshift then you must passed in inia[] argument"); 398 PetscCheck(ja, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "If oshift then you must passed in inja[] argument"); 399 for (i = 0; i <= (*m); i++) (*ia)[i]--; 400 for (i = 0; i < (*ia)[*m]; i++) (*ja)[i]--; 401 } 402 PetscFunctionReturn(PETSC_SUCCESS); 403 } 404 405 PetscErrorCode MatConvertFrom_MPIAdj(Mat A, MatType type, MatReuse reuse, Mat *newmat) 406 { 407 Mat B; 408 PetscInt i, m, N, nzeros = 0, *ia, *ja, len, rstart, cnt, j, *a; 409 const PetscInt *rj; 410 const PetscScalar *ra; 411 MPI_Comm comm; 412 413 PetscFunctionBegin; 414 PetscCall(MatGetSize(A, NULL, &N)); 415 PetscCall(MatGetLocalSize(A, &m, NULL)); 416 PetscCall(MatGetOwnershipRange(A, &rstart, NULL)); 417 418 /* count the number of nonzeros per row */ 419 for (i = 0; i < m; i++) { 420 PetscCall(MatGetRow(A, i + rstart, &len, &rj, NULL)); 421 for (j = 0; j < len; j++) { 422 if (rj[j] == i + rstart) { 423 len--; 424 break; 425 } /* don't count diagonal */ 426 } 427 nzeros += len; 428 PetscCall(MatRestoreRow(A, i + rstart, &len, &rj, NULL)); 429 } 430 431 /* malloc space for nonzeros */ 432 PetscCall(PetscMalloc1(nzeros + 1, &a)); 433 PetscCall(PetscMalloc1(N + 1, &ia)); 434 PetscCall(PetscMalloc1(nzeros + 1, &ja)); 435 436 nzeros = 0; 437 ia[0] = 0; 438 for (i = 0; i < m; i++) { 439 PetscCall(MatGetRow(A, i + rstart, &len, &rj, &ra)); 440 cnt = 0; 441 for (j = 0; j < len; j++) { 442 if (rj[j] != i + rstart) { /* if not diagonal */ 443 a[nzeros + cnt] = (PetscInt)PetscAbsScalar(ra[j]); 444 ja[nzeros + cnt++] = rj[j]; 445 } 446 } 447 PetscCall(MatRestoreRow(A, i + rstart, &len, &rj, &ra)); 448 nzeros += cnt; 449 ia[i + 1] = nzeros; 450 } 451 452 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 453 PetscCall(MatCreate(comm, &B)); 454 PetscCall(MatSetSizes(B, m, PETSC_DETERMINE, PETSC_DETERMINE, N)); 455 PetscCall(MatSetType(B, type)); 456 PetscCall(MatMPIAdjSetPreallocation(B, ia, ja, a)); 457 458 if (reuse == MAT_INPLACE_MATRIX) { 459 PetscCall(MatHeaderReplace(A, &B)); 460 } else { 461 *newmat = B; 462 } 463 PetscFunctionReturn(PETSC_SUCCESS); 464 } 465 466 PetscErrorCode MatSetValues_MPIAdj(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode im) 467 { 468 Mat_MPIAdj *adj = (Mat_MPIAdj *)A->data; 469 PetscInt rStart, rEnd, cStart, cEnd; 470 471 PetscFunctionBegin; 472 PetscCheck(!adj->i, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is already assembled, cannot change its values"); 473 PetscCall(MatGetOwnershipRange(A, &rStart, &rEnd)); 474 PetscCall(MatGetOwnershipRangeColumn(A, &cStart, &cEnd)); 475 if (!adj->ht) { 476 PetscCall(PetscHSetIJCreate(&adj->ht)); 477 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)A), 1, &A->stash)); 478 PetscCall(PetscLayoutSetUp(A->rmap)); 479 PetscCall(PetscLayoutSetUp(A->cmap)); 480 } 481 for (PetscInt r = 0; r < m; ++r) { 482 PetscHashIJKey key; 483 484 key.i = rows[r]; 485 if (key.i < 0) continue; 486 if ((key.i < rStart) || (key.i >= rEnd)) { 487 PetscCall(MatStashValuesRow_Private(&A->stash, key.i, n, cols, values, PETSC_FALSE)); 488 } else { 489 for (PetscInt c = 0; c < n; ++c) { 490 key.j = cols[c]; 491 if (key.j < 0 || key.i == key.j) continue; 492 PetscCall(PetscHSetIJAdd(adj->ht, key)); 493 } 494 } 495 } 496 PetscFunctionReturn(PETSC_SUCCESS); 497 } 498 499 PetscErrorCode MatAssemblyBegin_MPIAdj(Mat A, MatAssemblyType type) 500 { 501 PetscInt nstash, reallocs; 502 Mat_MPIAdj *adj = (Mat_MPIAdj *)A->data; 503 504 PetscFunctionBegin; 505 if (!adj->ht) { 506 PetscCall(PetscHSetIJCreate(&adj->ht)); 507 PetscCall(PetscLayoutSetUp(A->rmap)); 508 PetscCall(PetscLayoutSetUp(A->cmap)); 509 } 510 PetscCall(MatStashScatterBegin_Private(A, &A->stash, A->rmap->range)); 511 PetscCall(MatStashGetInfo_Private(&A->stash, &nstash, &reallocs)); 512 PetscCall(PetscInfo(A, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs)); 513 PetscFunctionReturn(PETSC_SUCCESS); 514 } 515 516 PetscErrorCode MatAssemblyEnd_MPIAdj(Mat A, MatAssemblyType type) 517 { 518 PetscScalar *val; 519 PetscInt *row, *col, m, rstart, *rowstarts; 520 PetscInt i, j, ncols, flg, nz; 521 PetscMPIInt n; 522 Mat_MPIAdj *adj = (Mat_MPIAdj *)A->data; 523 PetscHashIter hi; 524 PetscHashIJKey key; 525 PetscHSetIJ ht = adj->ht; 526 527 PetscFunctionBegin; 528 while (1) { 529 PetscCall(MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg)); 530 if (!flg) break; 531 532 for (i = 0; i < n;) { 533 /* Identify the consecutive vals belonging to the same row */ 534 for (j = i, rstart = row[j]; j < n; j++) { 535 if (row[j] != rstart) break; 536 } 537 if (j < n) ncols = j - i; 538 else ncols = n - i; 539 /* Set all these values with a single function call */ 540 PetscCall(MatSetValues_MPIAdj(A, 1, row + i, ncols, col + i, val + i, INSERT_VALUES)); 541 i = j; 542 } 543 } 544 PetscCall(MatStashScatterEnd_Private(&A->stash)); 545 PetscCall(MatStashDestroy_Private(&A->stash)); 546 547 PetscCall(MatGetLocalSize(A, &m, NULL)); 548 PetscCall(MatGetOwnershipRange(A, &rstart, NULL)); 549 PetscCall(PetscCalloc1(m + 1, &rowstarts)); 550 PetscHashIterBegin(ht, hi); 551 for (; !PetscHashIterAtEnd(ht, hi);) { 552 PetscHashIterGetKey(ht, hi, key); 553 rowstarts[key.i - rstart + 1]++; 554 PetscHashIterNext(ht, hi); 555 } 556 for (i = 1; i < m + 1; i++) rowstarts[i] = rowstarts[i - 1] + rowstarts[i]; 557 558 PetscCall(PetscHSetIJGetSize(ht, &nz)); 559 PetscCall(PetscMalloc1(nz, &col)); 560 PetscHashIterBegin(ht, hi); 561 for (; !PetscHashIterAtEnd(ht, hi);) { 562 PetscHashIterGetKey(ht, hi, key); 563 col[rowstarts[key.i - rstart]++] = key.j; 564 PetscHashIterNext(ht, hi); 565 } 566 PetscCall(PetscHSetIJDestroy(&ht)); 567 for (i = m; i > 0; i--) rowstarts[i] = rowstarts[i - 1]; 568 rowstarts[0] = 0; 569 570 for (PetscInt i = 0; i < m; i++) PetscCall(PetscSortInt(rowstarts[i + 1] - rowstarts[i], &col[rowstarts[i]])); 571 572 adj->i = rowstarts; 573 adj->j = col; 574 adj->nz = rowstarts[m]; 575 adj->freeaij = PETSC_TRUE; 576 PetscFunctionReturn(PETSC_SUCCESS); 577 } 578 579 static struct _MatOps MatOps_Values = {MatSetValues_MPIAdj, 580 MatGetRow_MPIAdj, 581 MatRestoreRow_MPIAdj, 582 NULL, 583 /* 4*/ NULL, 584 NULL, 585 NULL, 586 NULL, 587 NULL, 588 NULL, 589 /*10*/ NULL, 590 NULL, 591 NULL, 592 NULL, 593 NULL, 594 /*15*/ NULL, 595 MatEqual_MPIAdj, 596 NULL, 597 NULL, 598 NULL, 599 /*20*/ MatAssemblyBegin_MPIAdj, 600 MatAssemblyEnd_MPIAdj, 601 MatSetOption_MPIAdj, 602 NULL, 603 /*24*/ NULL, 604 NULL, 605 NULL, 606 NULL, 607 NULL, 608 /*29*/ NULL, 609 NULL, 610 NULL, 611 NULL, 612 NULL, 613 /*34*/ NULL, 614 NULL, 615 NULL, 616 NULL, 617 NULL, 618 /*39*/ NULL, 619 MatCreateSubMatrices_MPIAdj, 620 NULL, 621 NULL, 622 NULL, 623 /*44*/ NULL, 624 NULL, 625 MatShift_Basic, 626 NULL, 627 NULL, 628 /*49*/ NULL, 629 MatGetRowIJ_MPIAdj, 630 MatRestoreRowIJ_MPIAdj, 631 NULL, 632 NULL, 633 /*54*/ NULL, 634 NULL, 635 NULL, 636 NULL, 637 NULL, 638 /*59*/ NULL, 639 MatDestroy_MPIAdj, 640 MatView_MPIAdj, 641 MatConvertFrom_MPIAdj, 642 NULL, 643 /*64*/ NULL, 644 NULL, 645 NULL, 646 NULL, 647 NULL, 648 /*69*/ NULL, 649 NULL, 650 NULL, 651 NULL, 652 NULL, 653 /*74*/ NULL, 654 NULL, 655 NULL, 656 NULL, 657 NULL, 658 /*79*/ NULL, 659 NULL, 660 NULL, 661 NULL, 662 NULL, 663 /*84*/ NULL, 664 NULL, 665 NULL, 666 NULL, 667 NULL, 668 /*89*/ NULL, 669 NULL, 670 NULL, 671 NULL, 672 NULL, 673 /*94*/ NULL, 674 NULL, 675 NULL, 676 NULL, 677 NULL, 678 /*99*/ NULL, 679 NULL, 680 NULL, 681 NULL, 682 NULL, 683 /*104*/ NULL, 684 NULL, 685 NULL, 686 NULL, 687 NULL, 688 /*109*/ NULL, 689 NULL, 690 NULL, 691 NULL, 692 NULL, 693 /*114*/ NULL, 694 NULL, 695 NULL, 696 NULL, 697 NULL, 698 /*119*/ NULL, 699 NULL, 700 NULL, 701 NULL, 702 NULL, 703 /*124*/ NULL, 704 NULL, 705 NULL, 706 NULL, 707 MatCreateSubMatricesMPI_MPIAdj, 708 /*129*/ NULL, 709 NULL, 710 NULL, 711 NULL, 712 NULL, 713 /*134*/ NULL, 714 NULL, 715 NULL, 716 NULL, 717 NULL, 718 /*139*/ NULL, 719 NULL, 720 NULL, 721 NULL, 722 NULL, 723 /*144*/ NULL, 724 NULL, 725 NULL, 726 NULL, 727 NULL, 728 NULL, 729 /*150*/ NULL, 730 NULL}; 731 732 static PetscErrorCode MatMPIAdjSetPreallocation_MPIAdj(Mat B, PetscInt *i, PetscInt *j, PetscInt *values) 733 { 734 Mat_MPIAdj *b = (Mat_MPIAdj *)B->data; 735 PetscBool useedgeweights; 736 737 PetscFunctionBegin; 738 PetscCall(PetscLayoutSetUp(B->rmap)); 739 PetscCall(PetscLayoutSetUp(B->cmap)); 740 if (values) useedgeweights = PETSC_TRUE; 741 else useedgeweights = PETSC_FALSE; 742 /* Make everybody knows if they are using edge weights or not */ 743 PetscCall(MPIU_Allreduce((int *)&useedgeweights, (int *)&b->useedgeweights, 1, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)B))); 744 745 if (PetscDefined(USE_DEBUG)) { 746 PetscInt ii; 747 748 PetscCheck(i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "First i[] index must be zero, instead it is %" PetscInt_FMT, i[0]); 749 for (ii = 1; ii < B->rmap->n; ii++) { 750 PetscCheck(i[ii] >= 0 && i[ii] >= i[ii - 1], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i[%" PetscInt_FMT "]=%" PetscInt_FMT " index is out of range: i[%" PetscInt_FMT "]=%" PetscInt_FMT, ii, i[ii], ii - 1, i[ii - 1]); 751 } 752 for (ii = 0; ii < i[B->rmap->n]; ii++) PetscCheck(j[ii] >= 0 && j[ii] < B->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column index %" PetscInt_FMT " out of range %" PetscInt_FMT, ii, j[ii]); 753 } 754 b->j = j; 755 b->i = i; 756 b->values = values; 757 758 b->nz = i[B->rmap->n]; 759 b->diag = NULL; 760 b->symmetric = PETSC_FALSE; 761 b->freeaij = PETSC_TRUE; 762 763 B->ops->assemblybegin = NULL; 764 B->ops->assemblyend = NULL; 765 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 766 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 767 PetscCall(MatStashDestroy_Private(&B->stash)); 768 PetscFunctionReturn(PETSC_SUCCESS); 769 } 770 771 static PetscErrorCode MatMPIAdjCreateNonemptySubcommMat_MPIAdj(Mat A, Mat *B) 772 { 773 Mat_MPIAdj *a = (Mat_MPIAdj *)A->data; 774 const PetscInt *ranges; 775 MPI_Comm acomm, bcomm; 776 MPI_Group agroup, bgroup; 777 PetscMPIInt i, rank, size, nranks, *ranks; 778 779 PetscFunctionBegin; 780 *B = NULL; 781 PetscCall(PetscObjectGetComm((PetscObject)A, &acomm)); 782 PetscCallMPI(MPI_Comm_size(acomm, &size)); 783 PetscCallMPI(MPI_Comm_size(acomm, &rank)); 784 PetscCall(MatGetOwnershipRanges(A, &ranges)); 785 for (i = 0, nranks = 0; i < size; i++) { 786 if (ranges[i + 1] - ranges[i] > 0) nranks++; 787 } 788 if (nranks == size) { /* All ranks have a positive number of rows, so we do not need to create a subcomm; */ 789 PetscCall(PetscObjectReference((PetscObject)A)); 790 *B = A; 791 PetscFunctionReturn(PETSC_SUCCESS); 792 } 793 794 PetscCall(PetscMalloc1(nranks, &ranks)); 795 for (i = 0, nranks = 0; i < size; i++) { 796 if (ranges[i + 1] - ranges[i] > 0) ranks[nranks++] = i; 797 } 798 PetscCallMPI(MPI_Comm_group(acomm, &agroup)); 799 PetscCallMPI(MPI_Group_incl(agroup, nranks, ranks, &bgroup)); 800 PetscCall(PetscFree(ranks)); 801 PetscCallMPI(MPI_Comm_create(acomm, bgroup, &bcomm)); 802 PetscCallMPI(MPI_Group_free(&agroup)); 803 PetscCallMPI(MPI_Group_free(&bgroup)); 804 if (bcomm != MPI_COMM_NULL) { 805 PetscInt m, N; 806 Mat_MPIAdj *b; 807 PetscCall(MatGetLocalSize(A, &m, NULL)); 808 PetscCall(MatGetSize(A, NULL, &N)); 809 PetscCall(MatCreateMPIAdj(bcomm, m, N, a->i, a->j, a->values, B)); 810 b = (Mat_MPIAdj *)(*B)->data; 811 b->freeaij = PETSC_FALSE; 812 PetscCallMPI(MPI_Comm_free(&bcomm)); 813 } 814 PetscFunctionReturn(PETSC_SUCCESS); 815 } 816 817 PetscErrorCode MatMPIAdjToSeq_MPIAdj(Mat A, Mat *B) 818 { 819 PetscInt M, N, *II, *J, NZ, nz, m, nzstart, i; 820 PetscInt *Values = NULL; 821 Mat_MPIAdj *adj = (Mat_MPIAdj *)A->data; 822 PetscMPIInt mnz, mm, *allnz, *allm, size, *dispnz, *dispm; 823 824 PetscFunctionBegin; 825 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 826 PetscCall(MatGetSize(A, &M, &N)); 827 PetscCall(MatGetLocalSize(A, &m, NULL)); 828 nz = adj->nz; 829 PetscCheck(adj->i[m] == nz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nz %" PetscInt_FMT " not correct i[m] %" PetscInt_FMT, nz, adj->i[m]); 830 PetscCall(MPIU_Allreduce(&nz, &NZ, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A))); 831 832 PetscCall(PetscMPIIntCast(nz, &mnz)); 833 PetscCall(PetscMalloc2(size, &allnz, size, &dispnz)); 834 PetscCallMPI(MPI_Allgather(&mnz, 1, MPI_INT, allnz, 1, MPI_INT, PetscObjectComm((PetscObject)A))); 835 dispnz[0] = 0; 836 for (i = 1; i < size; i++) dispnz[i] = dispnz[i - 1] + allnz[i - 1]; 837 if (adj->values) { 838 PetscCall(PetscMalloc1(NZ, &Values)); 839 PetscCallMPI(MPI_Allgatherv(adj->values, mnz, MPIU_INT, Values, allnz, dispnz, MPIU_INT, PetscObjectComm((PetscObject)A))); 840 } 841 PetscCall(PetscMalloc1(NZ, &J)); 842 PetscCallMPI(MPI_Allgatherv(adj->j, mnz, MPIU_INT, J, allnz, dispnz, MPIU_INT, PetscObjectComm((PetscObject)A))); 843 PetscCall(PetscFree2(allnz, dispnz)); 844 PetscCallMPI(MPI_Scan(&nz, &nzstart, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A))); 845 nzstart -= nz; 846 /* shift the i[] values so they will be correct after being received */ 847 for (i = 0; i < m; i++) adj->i[i] += nzstart; 848 PetscCall(PetscMalloc1(M + 1, &II)); 849 PetscCall(PetscMPIIntCast(m, &mm)); 850 PetscCall(PetscMalloc2(size, &allm, size, &dispm)); 851 PetscCallMPI(MPI_Allgather(&mm, 1, MPI_INT, allm, 1, MPI_INT, PetscObjectComm((PetscObject)A))); 852 dispm[0] = 0; 853 for (i = 1; i < size; i++) dispm[i] = dispm[i - 1] + allm[i - 1]; 854 PetscCallMPI(MPI_Allgatherv(adj->i, mm, MPIU_INT, II, allm, dispm, MPIU_INT, PetscObjectComm((PetscObject)A))); 855 PetscCall(PetscFree2(allm, dispm)); 856 II[M] = NZ; 857 /* shift the i[] values back */ 858 for (i = 0; i < m; i++) adj->i[i] -= nzstart; 859 PetscCall(MatCreateMPIAdj(PETSC_COMM_SELF, M, N, II, J, Values, B)); 860 PetscFunctionReturn(PETSC_SUCCESS); 861 } 862 863 PetscErrorCode MatMPIAdjToSeqRankZero_MPIAdj(Mat A, Mat *B) 864 { 865 PetscInt M, N, *II, *J, NZ, nz, m, nzstart, i; 866 PetscInt *Values = NULL; 867 Mat_MPIAdj *adj = (Mat_MPIAdj *)A->data; 868 PetscMPIInt mnz, mm, *allnz = NULL, *allm, size, *dispnz, *dispm, rank; 869 870 PetscFunctionBegin; 871 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 872 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank)); 873 PetscCall(MatGetSize(A, &M, &N)); 874 PetscCall(MatGetLocalSize(A, &m, NULL)); 875 nz = adj->nz; 876 PetscCheck(adj->i[m] == nz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nz %" PetscInt_FMT " not correct i[m] %" PetscInt_FMT, nz, adj->i[m]); 877 PetscCall(MPIU_Allreduce(&nz, &NZ, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A))); 878 879 PetscCall(PetscMPIIntCast(nz, &mnz)); 880 if (!rank) PetscCall(PetscMalloc2(size, &allnz, size, &dispnz)); 881 PetscCallMPI(MPI_Gather(&mnz, 1, MPI_INT, allnz, 1, MPI_INT, 0, PetscObjectComm((PetscObject)A))); 882 if (!rank) { 883 dispnz[0] = 0; 884 for (i = 1; i < size; i++) dispnz[i] = dispnz[i - 1] + allnz[i - 1]; 885 if (adj->values) { 886 PetscCall(PetscMalloc1(NZ, &Values)); 887 PetscCallMPI(MPI_Gatherv(adj->values, mnz, MPIU_INT, Values, allnz, dispnz, MPIU_INT, 0, PetscObjectComm((PetscObject)A))); 888 } 889 PetscCall(PetscMalloc1(NZ, &J)); 890 PetscCallMPI(MPI_Gatherv(adj->j, mnz, MPIU_INT, J, allnz, dispnz, MPIU_INT, 0, PetscObjectComm((PetscObject)A))); 891 PetscCall(PetscFree2(allnz, dispnz)); 892 } else { 893 if (adj->values) PetscCallMPI(MPI_Gatherv(adj->values, mnz, MPIU_INT, NULL, NULL, NULL, MPIU_INT, 0, PetscObjectComm((PetscObject)A))); 894 PetscCallMPI(MPI_Gatherv(adj->j, mnz, MPIU_INT, NULL, NULL, NULL, MPIU_INT, 0, PetscObjectComm((PetscObject)A))); 895 } 896 PetscCallMPI(MPI_Scan(&nz, &nzstart, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A))); 897 nzstart -= nz; 898 /* shift the i[] values so they will be correct after being received */ 899 for (i = 0; i < m; i++) adj->i[i] += nzstart; 900 PetscCall(PetscMPIIntCast(m, &mm)); 901 if (!rank) { 902 PetscCall(PetscMalloc1(M + 1, &II)); 903 PetscCall(PetscMalloc2(size, &allm, size, &dispm)); 904 PetscCallMPI(MPI_Gather(&mm, 1, MPI_INT, allm, 1, MPI_INT, 0, PetscObjectComm((PetscObject)A))); 905 dispm[0] = 0; 906 for (i = 1; i < size; i++) dispm[i] = dispm[i - 1] + allm[i - 1]; 907 PetscCallMPI(MPI_Gatherv(adj->i, mm, MPIU_INT, II, allm, dispm, MPIU_INT, 0, PetscObjectComm((PetscObject)A))); 908 PetscCall(PetscFree2(allm, dispm)); 909 II[M] = NZ; 910 } else { 911 PetscCallMPI(MPI_Gather(&mm, 1, MPI_INT, NULL, 1, MPI_INT, 0, PetscObjectComm((PetscObject)A))); 912 PetscCallMPI(MPI_Gatherv(adj->i, mm, MPIU_INT, NULL, NULL, NULL, MPIU_INT, 0, PetscObjectComm((PetscObject)A))); 913 } 914 /* shift the i[] values back */ 915 for (i = 0; i < m; i++) adj->i[i] -= nzstart; 916 if (!rank) PetscCall(MatCreateMPIAdj(PETSC_COMM_SELF, M, N, II, J, Values, B)); 917 PetscFunctionReturn(PETSC_SUCCESS); 918 } 919 920 /*@ 921 MatMPIAdjCreateNonemptySubcommMat - create the same `MATMPIADJ` matrix on a subcommunicator containing only processes owning a positive number of rows 922 923 Collective 924 925 Input Parameter: 926 . A - original `MATMPIADJ` matrix 927 928 Output Parameter: 929 . B - matrix on subcommunicator, `NULL` on ranks that owned zero rows of `A` 930 931 Level: developer 932 933 Note: 934 This function is mostly useful for internal use by mesh partitioning packages that require that every process owns at least one row. 935 936 The matrix `B` should be destroyed with `MatDestroy()`. The arrays are not copied, so `B` should be destroyed before `A` is destroyed. 937 938 .seealso: [](ch_matrices), `Mat`, `MATMPIADJ`, `MatCreateMPIAdj()` 939 @*/ 940 PetscErrorCode MatMPIAdjCreateNonemptySubcommMat(Mat A, Mat *B) 941 { 942 PetscFunctionBegin; 943 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 944 PetscUseMethod(A, "MatMPIAdjCreateNonemptySubcommMat_C", (Mat, Mat *), (A, B)); 945 PetscFunctionReturn(PETSC_SUCCESS); 946 } 947 948 /*MC 949 MATMPIADJ - MATMPIADJ = "mpiadj" - A matrix type to be used for distributed adjacency matrices, 950 intended for use constructing orderings and partitionings. 951 952 Level: beginner 953 954 Note: 955 You can provide values to the matrix using `MatMPIAdjSetPreallocation()`, `MatCreateMPIAdj()`, or 956 by calling `MatSetValues()` and `MatAssemblyBegin()` followed by `MatAssemblyEnd()` 957 958 .seealso: [](ch_matrices), `Mat`, `MatCreateMPIAdj()`, `MatMPIAdjSetPreallocation()`, `MatSetValues()` 959 M*/ 960 PETSC_EXTERN PetscErrorCode MatCreate_MPIAdj(Mat B) 961 { 962 Mat_MPIAdj *b; 963 964 PetscFunctionBegin; 965 PetscCall(PetscNew(&b)); 966 B->data = (void *)b; 967 B->ops[0] = MatOps_Values; 968 B->assembled = PETSC_FALSE; 969 B->preallocated = PETSC_TRUE; /* so that MatSetValues() may be used */ 970 971 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPIAdjSetPreallocation_C", MatMPIAdjSetPreallocation_MPIAdj)); 972 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPIAdjCreateNonemptySubcommMat_C", MatMPIAdjCreateNonemptySubcommMat_MPIAdj)); 973 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPIAdjToSeq_C", MatMPIAdjToSeq_MPIAdj)); 974 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMPIAdjToSeqRankZero_C", MatMPIAdjToSeqRankZero_MPIAdj)); 975 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATMPIADJ)); 976 PetscFunctionReturn(PETSC_SUCCESS); 977 } 978 979 /*@C 980 MatMPIAdjToSeq - Converts an parallel `MATMPIADJ` matrix to complete `MATMPIADJ` on each process (needed by sequential partitioner) 981 982 Logically Collective 983 984 Input Parameter: 985 . A - the matrix 986 987 Output Parameter: 988 . B - the same matrix on all processes 989 990 Level: intermediate 991 992 .seealso: [](ch_matrices), `Mat`, `MATMPIADJ`, `MatCreate()`, `MatCreateMPIAdj()`, `MatSetValues()`, `MatMPIAdjToSeqRankZero()` 993 @*/ 994 PetscErrorCode MatMPIAdjToSeq(Mat A, Mat *B) 995 { 996 PetscFunctionBegin; 997 PetscUseMethod(A, "MatMPIAdjToSeq_C", (Mat, Mat *), (A, B)); 998 PetscFunctionReturn(PETSC_SUCCESS); 999 } 1000 1001 /*@C 1002 MatMPIAdjToSeqRankZero - Converts an parallel `MATMPIADJ` matrix to complete `MATMPIADJ` on rank zero (needed by sequential partitioner) 1003 1004 Logically Collective 1005 1006 Input Parameter: 1007 . A - the matrix 1008 1009 Output Parameter: 1010 . B - the same matrix on rank zero, not set on other ranks 1011 1012 Level: intermediate 1013 1014 Note: 1015 This routine has the advantage on systems with multiple ranks per node since only one copy of the matrix 1016 is stored on the first node, instead of the number of ranks copies. This can allow partitioning much larger 1017 parallel graph sequentially. 1018 1019 .seealso: [](ch_matrices), `Mat`, `MATMPIADJ`, `MatCreate()`, `MatCreateMPIAdj()`, `MatSetValues()`, `MatMPIAdjToSeq()` 1020 @*/ 1021 PetscErrorCode MatMPIAdjToSeqRankZero(Mat A, Mat *B) 1022 { 1023 PetscFunctionBegin; 1024 PetscUseMethod(A, "MatMPIAdjToSeqRankZero_C", (Mat, Mat *), (A, B)); 1025 PetscFunctionReturn(PETSC_SUCCESS); 1026 } 1027 1028 /*@C 1029 MatMPIAdjSetPreallocation - Sets the array used for storing the matrix elements 1030 1031 Logically Collective 1032 1033 Input Parameters: 1034 + A - the matrix 1035 . i - the indices into `j` for the start of each row 1036 . j - the column indices for each row (sorted for each row). 1037 The indices in `i` and `j` start with zero (NOT with one). 1038 - values - [use `NULL` if not provided] edge weights 1039 1040 Level: intermediate 1041 1042 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateMPIAdj()`, `MatSetValues()`, `MATMPIADJ`, `MatCreateMPIAdj()` 1043 @*/ 1044 PetscErrorCode MatMPIAdjSetPreallocation(Mat B, PetscInt *i, PetscInt *j, PetscInt *values) 1045 { 1046 PetscFunctionBegin; 1047 PetscTryMethod(B, "MatMPIAdjSetPreallocation_C", (Mat, PetscInt *, PetscInt *, PetscInt *), (B, i, j, values)); 1048 PetscFunctionReturn(PETSC_SUCCESS); 1049 } 1050 1051 /*@C 1052 MatCreateMPIAdj - Creates a sparse matrix representing an adjacency list. 1053 The matrix does not have numerical values associated with it, but is 1054 intended for ordering (to reduce bandwidth etc) and partitioning. 1055 1056 Collective 1057 1058 Input Parameters: 1059 + comm - MPI communicator 1060 . m - number of local rows 1061 . N - number of global columns 1062 . i - the indices into `j` for the start of each row 1063 . j - the column indices for each row (sorted for each row). 1064 The indices in `i` and `j` start with zero (NOT with one). 1065 - values -[use `NULL` if not provided] edge weights 1066 1067 Output Parameter: 1068 . A - the matrix 1069 1070 Level: intermediate 1071 1072 Notes: 1073 You must NOT free the `i`, `values` and `j` arrays yourself. PETSc will free them 1074 when the matrix is destroyed; you must allocate them with `PetscMalloc()`. If you 1075 call from Fortran you need not create the arrays with `PetscMalloc()`. 1076 1077 You should not include the matrix diagonals. 1078 1079 If you already have a matrix, you can create its adjacency matrix by a call 1080 to `MatConvert()`, specifying a type of `MATMPIADJ`. 1081 1082 Possible values for `MatSetOption()` - `MAT_STRUCTURALLY_SYMMETRIC` 1083 1084 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatConvert()`, `MatGetOrdering()`, `MATMPIADJ`, `MatMPIAdjSetPreallocation()` 1085 @*/ 1086 PetscErrorCode MatCreateMPIAdj(MPI_Comm comm, PetscInt m, PetscInt N, PetscInt *i, PetscInt *j, PetscInt *values, Mat *A) 1087 { 1088 PetscFunctionBegin; 1089 PetscCall(MatCreate(comm, A)); 1090 PetscCall(MatSetSizes(*A, m, PETSC_DETERMINE, PETSC_DETERMINE, N)); 1091 PetscCall(MatSetType(*A, MATMPIADJ)); 1092 PetscCall(MatMPIAdjSetPreallocation(*A, i, j, values)); 1093 PetscFunctionReturn(PETSC_SUCCESS); 1094 } 1095