1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/isimpl.h> 3 #include <petscsf.h> 4 #include <petscds.h> 5 6 /* get adjacencies due to point-to-point constraints that can't be found with DMPlexGetAdjacency() */ 7 static PetscErrorCode DMPlexComputeAnchorAdjacencies(DM dm, PetscBool useCone, PetscBool useClosure, PetscSection *anchorSectionAdj, PetscInt *anchorAdj[]) 8 { 9 PetscInt pStart, pEnd; 10 PetscSection section, sectionGlobal, adjSec, aSec; 11 IS aIS; 12 13 PetscFunctionBegin; 14 PetscCall(DMGetLocalSection(dm, §ion)); 15 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 16 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &adjSec)); 17 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 18 PetscCall(PetscSectionSetChart(adjSec, pStart, pEnd)); 19 20 PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS)); 21 if (aSec) { 22 const PetscInt *anchors; 23 PetscInt p, q, a, aSize, *offsets, aStart, aEnd, *inverse, iSize, *adj, adjSize; 24 PetscInt *tmpAdjP = NULL, *tmpAdjQ = NULL; 25 PetscSection inverseSec; 26 27 /* invert the constraint-to-anchor map */ 28 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)aSec), &inverseSec)); 29 PetscCall(PetscSectionSetChart(inverseSec, pStart, pEnd)); 30 PetscCall(ISGetLocalSize(aIS, &aSize)); 31 PetscCall(ISGetIndices(aIS, &anchors)); 32 33 for (p = 0; p < aSize; p++) { 34 PetscInt a = anchors[p]; 35 36 PetscCall(PetscSectionAddDof(inverseSec, a, 1)); 37 } 38 PetscCall(PetscSectionSetUp(inverseSec)); 39 PetscCall(PetscSectionGetStorageSize(inverseSec, &iSize)); 40 PetscCall(PetscMalloc1(iSize, &inverse)); 41 PetscCall(PetscCalloc1(pEnd - pStart, &offsets)); 42 PetscCall(PetscSectionGetChart(aSec, &aStart, &aEnd)); 43 for (p = aStart; p < aEnd; p++) { 44 PetscInt dof, off; 45 46 PetscCall(PetscSectionGetDof(aSec, p, &dof)); 47 PetscCall(PetscSectionGetOffset(aSec, p, &off)); 48 49 for (q = 0; q < dof; q++) { 50 PetscInt iOff; 51 52 a = anchors[off + q]; 53 PetscCall(PetscSectionGetOffset(inverseSec, a, &iOff)); 54 inverse[iOff + offsets[a - pStart]++] = p; 55 } 56 } 57 PetscCall(ISRestoreIndices(aIS, &anchors)); 58 PetscCall(PetscFree(offsets)); 59 60 /* construct anchorAdj and adjSec 61 * 62 * loop over anchors: 63 * construct anchor adjacency 64 * loop over constrained: 65 * construct constrained adjacency 66 * if not in anchor adjacency, add to dofs 67 * setup adjSec, allocate anchorAdj 68 * loop over anchors: 69 * construct anchor adjacency 70 * loop over constrained: 71 * construct constrained adjacency 72 * if not in anchor adjacency 73 * if not already in list, put in list 74 * sort, unique, reduce dof count 75 * optional: compactify 76 */ 77 for (p = pStart; p < pEnd; p++) { 78 PetscInt iDof, iOff, i, r, s, numAdjP = PETSC_DETERMINE; 79 80 PetscCall(PetscSectionGetDof(inverseSec, p, &iDof)); 81 if (!iDof) continue; 82 PetscCall(PetscSectionGetOffset(inverseSec, p, &iOff)); 83 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, PETSC_TRUE, &numAdjP, &tmpAdjP)); 84 for (i = 0; i < iDof; i++) { 85 PetscInt iNew = 0, qAdj, qAdjDof, qAdjCDof, numAdjQ = PETSC_DETERMINE; 86 87 q = inverse[iOff + i]; 88 PetscCall(DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ)); 89 for (r = 0; r < numAdjQ; r++) { 90 qAdj = tmpAdjQ[r]; 91 if ((qAdj < pStart) || (qAdj >= pEnd)) continue; 92 for (s = 0; s < numAdjP; s++) { 93 if (qAdj == tmpAdjP[s]) break; 94 } 95 if (s < numAdjP) continue; 96 PetscCall(PetscSectionGetDof(section, qAdj, &qAdjDof)); 97 PetscCall(PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof)); 98 iNew += qAdjDof - qAdjCDof; 99 } 100 PetscCall(PetscSectionAddDof(adjSec, p, iNew)); 101 } 102 } 103 104 PetscCall(PetscSectionSetUp(adjSec)); 105 PetscCall(PetscSectionGetStorageSize(adjSec, &adjSize)); 106 PetscCall(PetscMalloc1(adjSize, &adj)); 107 108 for (p = pStart; p < pEnd; p++) { 109 PetscInt iDof, iOff, i, r, s, aOff, aOffOrig, aDof, numAdjP = PETSC_DETERMINE; 110 111 PetscCall(PetscSectionGetDof(inverseSec, p, &iDof)); 112 if (!iDof) continue; 113 PetscCall(PetscSectionGetOffset(inverseSec, p, &iOff)); 114 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, PETSC_TRUE, &numAdjP, &tmpAdjP)); 115 PetscCall(PetscSectionGetDof(adjSec, p, &aDof)); 116 PetscCall(PetscSectionGetOffset(adjSec, p, &aOff)); 117 aOffOrig = aOff; 118 for (i = 0; i < iDof; i++) { 119 PetscInt qAdj, qAdjDof, qAdjCDof, qAdjOff, nd, numAdjQ = PETSC_DETERMINE; 120 121 q = inverse[iOff + i]; 122 PetscCall(DMPlexGetAdjacency_Internal(dm, q, useCone, useClosure, PETSC_TRUE, &numAdjQ, &tmpAdjQ)); 123 for (r = 0; r < numAdjQ; r++) { 124 qAdj = tmpAdjQ[r]; 125 if ((qAdj < pStart) || (qAdj >= pEnd)) continue; 126 for (s = 0; s < numAdjP; s++) { 127 if (qAdj == tmpAdjP[s]) break; 128 } 129 if (s < numAdjP) continue; 130 PetscCall(PetscSectionGetDof(section, qAdj, &qAdjDof)); 131 PetscCall(PetscSectionGetConstraintDof(section, qAdj, &qAdjCDof)); 132 PetscCall(PetscSectionGetOffset(sectionGlobal, qAdj, &qAdjOff)); 133 for (nd = 0; nd < qAdjDof - qAdjCDof; ++nd) adj[aOff++] = (qAdjOff < 0 ? -(qAdjOff + 1) : qAdjOff) + nd; 134 } 135 } 136 PetscCall(PetscSortRemoveDupsInt(&aDof, PetscSafePointerPlusOffset(adj, aOffOrig))); 137 PetscCall(PetscSectionSetDof(adjSec, p, aDof)); 138 } 139 *anchorAdj = adj; 140 141 /* clean up */ 142 PetscCall(PetscSectionDestroy(&inverseSec)); 143 PetscCall(PetscFree(inverse)); 144 PetscCall(PetscFree(tmpAdjP)); 145 PetscCall(PetscFree(tmpAdjQ)); 146 } else { 147 *anchorAdj = NULL; 148 PetscCall(PetscSectionSetUp(adjSec)); 149 } 150 *anchorSectionAdj = adjSec; 151 PetscFunctionReturn(PETSC_SUCCESS); 152 } 153 154 // Based off of `PetscIntViewNumColumns()` 155 static PetscErrorCode PetscIntViewPairs(PetscInt N, PetscInt Ncol, const PetscInt idx1[], const PetscInt idx2[], PetscViewer viewer) 156 { 157 PetscMPIInt rank, size; 158 PetscInt j, i, n = N / Ncol, p = N % Ncol; 159 PetscBool isascii; 160 MPI_Comm comm; 161 162 PetscFunctionBegin; 163 if (!viewer) viewer = PETSC_VIEWER_STDOUT_SELF; 164 if (N) PetscAssertPointer(idx1, 3); 165 if (N) PetscAssertPointer(idx2, 4); 166 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 5); 167 PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm)); 168 PetscCallMPI(MPI_Comm_size(comm, &size)); 169 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 170 171 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 172 if (isascii) { 173 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 174 for (i = 0; i < n; i++) { 175 if (size > 1) { 176 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT ":", rank, Ncol * i)); 177 } else { 178 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT ":", Ncol * i)); 179 } 180 for (j = 0; j < Ncol; j++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%" PetscInt_FMT ", %" PetscInt_FMT ")", idx1[i * Ncol + j], idx2[i * Ncol + j])); 181 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n")); 182 } 183 if (p) { 184 if (size > 1) { 185 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] %" PetscInt_FMT ":", rank, Ncol * n)); 186 } else { 187 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%" PetscInt_FMT ":", Ncol * n)); 188 } 189 for (i = 0; i < p; i++) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " (%" PetscInt_FMT ", %" PetscInt_FMT ")", idx1[Ncol * n + i], idx2[Ncol * n + i])); 190 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n")); 191 } 192 PetscCall(PetscViewerFlush(viewer)); 193 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 194 } else { 195 const char *tname; 196 PetscCall(PetscObjectGetName((PetscObject)viewer, &tname)); 197 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot handle that PetscViewer of type %s", tname); 198 } 199 PetscFunctionReturn(PETSC_SUCCESS); 200 } 201 202 // Determine if any of the local adjacencies match a leaf and root of the pointSF. 203 // When using isoperiodic boundary conditions, it is possible for periodic (leaf) and donor (root) pairs to be on the same rank. 204 // This check is done to ensure the adjacency in these cases is only counted for one of the mesh points rather than both. 205 static inline PetscErrorCode AdjancencyContainsLeafRootPair(PetscInt num_pairs, const PetscInt leaves[], const PetscInt roots[], PetscInt numAdj, const PetscInt tmpAdj[], PetscInt *num_leaves_found, PetscInt leaves_found[]) 206 { 207 PetscInt root_index = -1, leaf_, num_roots_found = 0; 208 209 PetscFunctionBeginUser; 210 *num_leaves_found = 0; 211 for (PetscInt q = 0; q < numAdj; q++) { 212 const PetscInt padj = tmpAdj[q]; 213 PetscCall(PetscFindInt(padj, num_pairs, roots, &root_index)); 214 if (root_index >= 0) { 215 leaves_found[num_roots_found] = root_index; // Initially use leaves_found to store pair indices 216 num_roots_found++; 217 break; 218 } 219 } 220 if (num_roots_found == 0) PetscFunctionReturn(PETSC_SUCCESS); 221 222 for (PetscInt i = 0; i < num_roots_found; i++) { 223 leaf_ = leaves[root_index]; 224 for (PetscInt q = 0; q < numAdj; q++) { 225 if (tmpAdj[q] == leaf_) { 226 leaves_found[*num_leaves_found] = leaf_; 227 (*num_leaves_found)++; 228 continue; 229 } 230 } 231 } 232 233 PetscCall(PetscIntSortSemiOrdered(*num_leaves_found, leaves_found)); 234 PetscFunctionReturn(PETSC_SUCCESS); 235 } 236 237 static PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx) 238 { 239 MPI_Comm comm; 240 PetscMPIInt myrank; 241 PetscBool doCommLocal, doComm, debug = PETSC_FALSE; 242 PetscSF sf, sfAdj; 243 PetscSection section, sectionGlobal, leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj; 244 PetscInt nroots, nleaves, l, p, r; 245 const PetscInt *leaves; 246 const PetscSFNode *remotes; 247 PetscInt dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols; 248 PetscInt *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets, *rootsMyRankPair = NULL, *leavesMyRankPair = NULL; 249 PetscInt adjSize, numMyRankPair = 0; 250 251 PetscFunctionBegin; 252 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 253 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL)); 254 PetscCallMPI(MPI_Comm_rank(comm, &myrank)); 255 PetscCall(DMGetDimension(dm, &dim)); 256 PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf)); 257 PetscCall(DMGetLocalSection(dm, §ion)); 258 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 259 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 260 doCommLocal = nroots >= 0 ? PETSC_TRUE : PETSC_FALSE; 261 PetscCallMPI(MPIU_Allreduce(&doCommLocal, &doComm, 1, MPI_C_BOOL, MPI_LAND, comm)); 262 /* Create section for dof adjacency (dof ==> # adj dof) */ 263 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 264 PetscCall(PetscSectionGetStorageSize(section, &numDof)); 265 PetscCall(PetscSectionCreate(comm, &leafSectionAdj)); 266 PetscCall(PetscSectionSetChart(leafSectionAdj, 0, numDof)); 267 PetscCall(PetscSectionCreate(comm, &rootSectionAdj)); 268 PetscCall(PetscSectionSetChart(rootSectionAdj, 0, numDof)); 269 PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes)); 270 271 // Store leaf-root pairs if remote.rank is the current rank 272 if (nleaves >= 0) PetscCall(PetscMalloc2(nleaves, &rootsMyRankPair, nleaves, &leavesMyRankPair)); 273 for (PetscInt l = 0; l < nleaves; l++) { 274 if (remotes[l].rank == myrank) { 275 rootsMyRankPair[numMyRankPair] = remotes[l].index; 276 leavesMyRankPair[numMyRankPair] = leaves[l]; 277 numMyRankPair++; 278 } 279 } 280 PetscCall(PetscIntSortSemiOrdered(numMyRankPair, rootsMyRankPair)); 281 PetscCall(PetscIntSortSemiOrdered(numMyRankPair, leavesMyRankPair)); 282 if (debug) { 283 PetscCall(PetscPrintf(comm, "Root/leaf pairs on the same rank:\n")); 284 PetscCall(PetscIntViewPairs(numMyRankPair, 5, rootsMyRankPair, leavesMyRankPair, NULL)); 285 } 286 /* 287 section - maps points to (# dofs, local dofs) 288 sectionGlobal - maps points to (# dofs, global dofs) 289 leafSectionAdj - maps unowned local dofs to # adj dofs 290 rootSectionAdj - maps owned local dofs to # adj dofs 291 adj - adj global dofs indexed by leafSectionAdj 292 rootAdj - adj global dofs indexed by rootSectionAdj 293 sf - describes shared points across procs 294 sfDof - describes shared dofs across procs 295 sfAdj - describes shared adjacent dofs across procs 296 ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point. 297 (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj 298 (This is done in DMPlexComputeAnchorAdjacencies()) 299 1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj 300 Reduce those counts to rootSectionAdj (now redundantly counting some interface points) 301 2. Visit owned points on interface, count adjacencies placing in rootSectionAdj 302 Create sfAdj connecting rootSectionAdj and leafSectionAdj 303 3. Visit unowned points on interface, write adjacencies to adj 304 Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies) 305 4. Visit owned points on interface, write adjacencies to rootAdj 306 Remove redundancy in rootAdj 307 ** The last two traversals use transitive closure 308 5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj) 309 Allocate memory addressed by sectionAdj (cols) 310 6. Visit all owned points in the subdomain, insert dof adjacencies into cols 311 ** Knowing all the column adjacencies, check ownership and sum into dnz and onz 312 */ 313 PetscCall(DMPlexComputeAnchorAdjacencies(dm, useCone, useClosure, &anchorSectionAdj, &anchorAdj)); 314 for (l = 0; l < nleaves; ++l) { 315 PetscInt dof, off, d, q, anDof; 316 PetscInt p = leaves[l], numAdj = PETSC_DETERMINE; 317 318 if ((p < pStart) || (p >= pEnd)) continue; 319 PetscCall(PetscSectionGetDof(section, p, &dof)); 320 PetscCall(PetscSectionGetOffset(section, p, &off)); 321 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 322 for (q = 0; q < numAdj; ++q) { 323 const PetscInt padj = tmpAdj[q]; 324 PetscInt ndof, ncdof; 325 326 if ((padj < pStart) || (padj >= pEnd)) continue; 327 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 328 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 329 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, ndof - ncdof)); 330 } 331 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 332 if (anDof) { 333 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, anDof)); 334 } 335 } 336 PetscCall(PetscSectionSetUp(leafSectionAdj)); 337 if (debug) { 338 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n")); 339 PetscCall(PetscSectionView(leafSectionAdj, NULL)); 340 } 341 /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */ 342 if (doComm) { 343 PetscCall(PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM)); 344 PetscCall(PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM)); 345 PetscCall(PetscSectionInvalidateMaxDof_Internal(rootSectionAdj)); 346 } 347 if (debug) { 348 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots:\n")); 349 PetscCall(PetscSectionView(rootSectionAdj, NULL)); 350 } 351 /* Add in local adjacency sizes for owned dofs on interface (roots) */ 352 for (p = pStart; p < pEnd; ++p) { 353 PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof; 354 355 PetscCall(PetscSectionGetDof(section, p, &dof)); 356 PetscCall(PetscSectionGetOffset(section, p, &off)); 357 if (!dof) continue; 358 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 359 if (adof <= 0) continue; 360 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 361 for (q = 0; q < numAdj; ++q) { 362 const PetscInt padj = tmpAdj[q]; 363 PetscInt ndof, ncdof; 364 365 if ((padj < pStart) || (padj >= pEnd)) continue; 366 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 367 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 368 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, ndof - ncdof)); 369 } 370 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 371 if (anDof) { 372 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, anDof)); 373 } 374 } 375 PetscCall(PetscSectionSetUp(rootSectionAdj)); 376 if (debug) { 377 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after local additions:\n")); 378 PetscCall(PetscSectionView(rootSectionAdj, NULL)); 379 } 380 /* Create adj SF based on dof SF */ 381 PetscCall(PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets)); 382 PetscCall(PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj)); 383 PetscCall(PetscFree(remoteOffsets)); 384 if (debug) { 385 PetscCall(PetscPrintf(comm, "Adjacency SF for Preallocation:\n")); 386 PetscCall(PetscSFView(sfAdj, NULL)); 387 } 388 /* Create leaf adjacency */ 389 PetscCall(PetscSectionSetUp(leafSectionAdj)); 390 PetscCall(PetscSectionGetStorageSize(leafSectionAdj, &adjSize)); 391 PetscCall(PetscCalloc1(adjSize, &adj)); 392 for (l = 0; l < nleaves; ++l) { 393 PetscInt dof, off, d, q, anDof, anOff; 394 PetscInt p = leaves[l], numAdj = PETSC_DETERMINE; 395 396 if ((p < pStart) || (p >= pEnd)) continue; 397 PetscCall(PetscSectionGetDof(section, p, &dof)); 398 PetscCall(PetscSectionGetOffset(section, p, &off)); 399 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 400 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 401 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 402 for (d = off; d < off + dof; ++d) { 403 PetscInt aoff, i = 0; 404 405 PetscCall(PetscSectionGetOffset(leafSectionAdj, d, &aoff)); 406 for (q = 0; q < numAdj; ++q) { 407 const PetscInt padj = tmpAdj[q]; 408 PetscInt ndof, ncdof, ngoff, nd; 409 410 if ((padj < pStart) || (padj >= pEnd)) continue; 411 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 412 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 413 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 414 for (nd = 0; nd < ndof - ncdof; ++nd) { 415 adj[aoff + i] = (ngoff < 0 ? -(ngoff + 1) : ngoff) + nd; 416 ++i; 417 } 418 } 419 for (q = 0; q < anDof; q++) { 420 adj[aoff + i] = anchorAdj[anOff + q]; 421 ++i; 422 } 423 } 424 } 425 /* Debugging */ 426 if (debug) { 427 PetscCall(PetscPrintf(comm, "Leaf adjacency indices\n")); 428 PetscCall(PetscSectionArrayView(leafSectionAdj, adj, PETSC_INT, NULL)); 429 } 430 /* Gather adjacent indices to root */ 431 PetscCall(PetscSectionGetStorageSize(rootSectionAdj, &adjSize)); 432 PetscCall(PetscMalloc1(adjSize, &rootAdj)); 433 for (r = 0; r < adjSize; ++r) rootAdj[r] = -1; 434 if (doComm) { 435 const PetscInt *indegree; 436 PetscInt *remoteadj, radjsize = 0; 437 438 PetscCall(PetscSFComputeDegreeBegin(sfAdj, &indegree)); 439 PetscCall(PetscSFComputeDegreeEnd(sfAdj, &indegree)); 440 for (p = 0; p < adjSize; ++p) radjsize += indegree[p]; 441 PetscCall(PetscMalloc1(radjsize, &remoteadj)); 442 PetscCall(PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj)); 443 PetscCall(PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj)); 444 for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p - 1])) { 445 PetscInt s; 446 for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l + s] = remoteadj[r]; 447 } 448 PetscCheck(r == radjsize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, r, radjsize); 449 PetscCheck(l == adjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, l, adjSize); 450 PetscCall(PetscFree(remoteadj)); 451 } 452 PetscCall(PetscSFDestroy(&sfAdj)); 453 PetscCall(PetscFree(adj)); 454 /* Debugging */ 455 if (debug) { 456 PetscCall(PetscPrintf(comm, "Root adjacency indices after gather\n")); 457 PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL)); 458 } 459 /* Add in local adjacency indices for owned dofs on interface (roots) */ 460 for (p = pStart; p < pEnd; ++p) { 461 PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff; 462 463 PetscCall(PetscSectionGetDof(section, p, &dof)); 464 PetscCall(PetscSectionGetOffset(section, p, &off)); 465 if (!dof) continue; 466 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 467 if (adof <= 0) continue; 468 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 469 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 470 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 471 for (d = off; d < off + dof; ++d) { 472 PetscInt adof, aoff, i; 473 474 PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof)); 475 PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff)); 476 i = adof - 1; 477 for (q = 0; q < anDof; q++) { 478 rootAdj[aoff + i] = anchorAdj[anOff + q]; 479 --i; 480 } 481 for (q = 0; q < numAdj; ++q) { 482 const PetscInt padj = tmpAdj[q]; 483 PetscInt ndof, ncdof, ngoff, nd; 484 485 if ((padj < pStart) || (padj >= pEnd)) continue; 486 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 487 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 488 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 489 for (nd = 0; nd < ndof - ncdof; ++nd) { 490 rootAdj[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd; 491 --i; 492 } 493 } 494 } 495 } 496 /* Debugging */ 497 if (debug) { 498 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots before compression:\n")); 499 PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL)); 500 } 501 /* Compress indices */ 502 PetscCall(PetscSectionSetUp(rootSectionAdj)); 503 for (p = pStart; p < pEnd; ++p) { 504 PetscInt dof, cdof, off, d; 505 PetscInt adof, aoff; 506 507 PetscCall(PetscSectionGetDof(section, p, &dof)); 508 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 509 PetscCall(PetscSectionGetOffset(section, p, &off)); 510 if (!dof) continue; 511 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 512 if (adof <= 0) continue; 513 for (d = off; d < off + dof - cdof; ++d) { 514 PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof)); 515 PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff)); 516 PetscCall(PetscSortRemoveDupsInt(&adof, &rootAdj[aoff])); 517 PetscCall(PetscSectionSetDof(rootSectionAdj, d, adof)); 518 } 519 } 520 /* Debugging */ 521 if (debug) { 522 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after compression:\n")); 523 PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL)); 524 } 525 /* Build adjacency section: Maps global indices to sets of adjacent global indices */ 526 PetscCall(PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd)); 527 PetscCall(PetscSectionCreate(comm, §ionAdj)); 528 PetscCall(PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd)); 529 530 PetscInt *exclude_leaves, num_exclude_leaves, max_adjacency_size = 0; 531 PetscCall(DMPlexGetMaxAdjacencySize_Internal(dm, useAnchors, &max_adjacency_size)); 532 PetscCall(PetscMalloc1(max_adjacency_size, &exclude_leaves)); 533 534 for (p = pStart; p < pEnd; ++p) { 535 PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof; 536 PetscBool found = PETSC_TRUE; 537 538 PetscCall(PetscSectionGetDof(section, p, &dof)); 539 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 540 PetscCall(PetscSectionGetOffset(section, p, &off)); 541 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 542 for (d = 0; d < dof - cdof; ++d) { 543 PetscInt ldof, rdof; 544 545 PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof)); 546 PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof)); 547 if (ldof > 0) { 548 /* We do not own this point */ 549 } else if (rdof > 0) { 550 PetscCall(PetscSectionSetDof(sectionAdj, goff + d, rdof)); 551 } else { 552 found = PETSC_FALSE; 553 } 554 } 555 if (found) continue; 556 PetscCall(PetscSectionGetDof(section, p, &dof)); 557 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 558 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 559 PetscCall(AdjancencyContainsLeafRootPair(numMyRankPair, leavesMyRankPair, rootsMyRankPair, numAdj, tmpAdj, &num_exclude_leaves, exclude_leaves)); 560 for (q = 0; q < numAdj; ++q) { 561 const PetscInt padj = tmpAdj[q]; 562 PetscInt ndof, ncdof, noff, count; 563 564 if ((padj < pStart) || (padj >= pEnd)) continue; 565 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 566 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 567 PetscCall(PetscSectionGetOffset(section, padj, &noff)); 568 // If leaf-root pair are both on this rank, only count root 569 PetscCall(PetscFindInt(padj, num_exclude_leaves, exclude_leaves, &count)); 570 if (count >= 0) continue; 571 for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, ndof - ncdof)); 572 } 573 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 574 if (anDof) { 575 for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, anDof)); 576 } 577 } 578 PetscCall(PetscSectionSetUp(sectionAdj)); 579 if (debug) { 580 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation:\n")); 581 PetscCall(PetscSectionView(sectionAdj, NULL)); 582 } 583 /* Get adjacent indices */ 584 PetscCall(PetscSectionGetStorageSize(sectionAdj, &numCols)); 585 PetscCall(PetscMalloc1(numCols, &cols)); 586 for (p = pStart; p < pEnd; ++p) { 587 PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff; 588 PetscBool found = PETSC_TRUE; 589 590 PetscCall(PetscSectionGetDof(section, p, &dof)); 591 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 592 PetscCall(PetscSectionGetOffset(section, p, &off)); 593 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 594 for (d = 0; d < dof - cdof; ++d) { 595 PetscInt ldof, rdof; 596 597 PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof)); 598 PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof)); 599 if (ldof > 0) { 600 /* We do not own this point */ 601 } else if (rdof > 0) { 602 PetscInt aoff, roff; 603 604 PetscCall(PetscSectionGetOffset(sectionAdj, goff + d, &aoff)); 605 PetscCall(PetscSectionGetOffset(rootSectionAdj, off + d, &roff)); 606 PetscCall(PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof)); 607 } else { 608 found = PETSC_FALSE; 609 } 610 } 611 if (found) continue; 612 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 613 PetscCall(AdjancencyContainsLeafRootPair(numMyRankPair, leavesMyRankPair, rootsMyRankPair, numAdj, tmpAdj, &num_exclude_leaves, exclude_leaves)); 614 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 615 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 616 for (d = goff; d < goff + dof - cdof; ++d) { 617 PetscInt adof, aoff, i = 0; 618 619 PetscCall(PetscSectionGetDof(sectionAdj, d, &adof)); 620 PetscCall(PetscSectionGetOffset(sectionAdj, d, &aoff)); 621 for (q = 0; q < numAdj; ++q) { 622 const PetscInt padj = tmpAdj[q], *ncind; 623 PetscInt ndof, ncdof, ngoff, nd, count; 624 625 /* Adjacent points may not be in the section chart */ 626 if ((padj < pStart) || (padj >= pEnd)) continue; 627 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 628 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 629 PetscCall(PetscSectionGetConstraintIndices(section, padj, &ncind)); 630 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 631 // If leaf-root pair are both on this rank, only count root 632 PetscCall(PetscFindInt(padj, num_exclude_leaves, exclude_leaves, &count)); 633 if (count >= 0) continue; 634 for (nd = 0; nd < ndof - ncdof; ++nd, ++i) cols[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd; 635 } 636 for (q = 0; q < anDof; q++, i++) cols[aoff + i] = anchorAdj[anOff + q]; 637 PetscCheck(i == adof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of entries %" PetscInt_FMT " != %" PetscInt_FMT " for dof %" PetscInt_FMT " (point %" PetscInt_FMT ")", i, adof, d, p); 638 } 639 } 640 PetscCall(PetscSectionDestroy(&anchorSectionAdj)); 641 PetscCall(PetscSectionDestroy(&leafSectionAdj)); 642 PetscCall(PetscSectionDestroy(&rootSectionAdj)); 643 PetscCall(PetscFree(exclude_leaves)); 644 PetscCall(PetscFree2(rootsMyRankPair, leavesMyRankPair)); 645 PetscCall(PetscFree(anchorAdj)); 646 PetscCall(PetscFree(rootAdj)); 647 PetscCall(PetscFree(tmpAdj)); 648 /* Debugging */ 649 if (debug) { 650 PetscCall(PetscPrintf(comm, "Column indices\n")); 651 PetscCall(PetscSectionArrayView(sectionAdj, cols, PETSC_INT, NULL)); 652 } 653 654 *sA = sectionAdj; 655 *colIdx = cols; 656 PetscFunctionReturn(PETSC_SUCCESS); 657 } 658 659 static PetscErrorCode DMPlexUpdateAllocation_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[]) 660 { 661 PetscSection section; 662 PetscInt rStart, rEnd, r, pStart, pEnd, p; 663 664 PetscFunctionBegin; 665 /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */ 666 PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd)); 667 PetscCheck((rStart % bs) == 0 && (rEnd % bs) == 0, PetscObjectComm((PetscObject)rLayout), PETSC_ERR_ARG_WRONG, "Invalid layout [%" PetscInt_FMT ", %" PetscInt_FMT ") for matrix, must be divisible by block size %" PetscInt_FMT, rStart, rEnd, bs); 668 if (f >= 0 && bs == 1) { 669 PetscCall(DMGetLocalSection(dm, §ion)); 670 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 671 for (p = pStart; p < pEnd; ++p) { 672 PetscInt rS, rE; 673 674 PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE)); 675 for (r = rS; r < rE; ++r) { 676 PetscInt numCols, cStart, c; 677 678 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 679 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 680 for (c = cStart; c < cStart + numCols; ++c) { 681 if ((cols[c] >= rStart) && (cols[c] < rEnd)) { 682 ++dnz[r - rStart]; 683 if (cols[c] >= r) ++dnzu[r - rStart]; 684 } else { 685 ++onz[r - rStart]; 686 if (cols[c] >= r) ++onzu[r - rStart]; 687 } 688 } 689 } 690 } 691 } else { 692 /* Only loop over blocks of rows */ 693 for (r = rStart / bs; r < rEnd / bs; ++r) { 694 const PetscInt row = r * bs; 695 PetscInt numCols, cStart, c; 696 697 PetscCall(PetscSectionGetDof(sectionAdj, row, &numCols)); 698 PetscCall(PetscSectionGetOffset(sectionAdj, row, &cStart)); 699 for (c = cStart; c < cStart + numCols; ++c) { 700 if ((cols[c] >= rStart) && (cols[c] < rEnd)) { 701 ++dnz[r - rStart / bs]; 702 if (cols[c] >= row) ++dnzu[r - rStart / bs]; 703 } else { 704 ++onz[r - rStart / bs]; 705 if (cols[c] >= row) ++onzu[r - rStart / bs]; 706 } 707 } 708 } 709 for (r = 0; r < (rEnd - rStart) / bs; ++r) { 710 dnz[r] /= bs; 711 onz[r] /= bs; 712 dnzu[r] /= bs; 713 onzu[r] /= bs; 714 } 715 } 716 PetscFunctionReturn(PETSC_SUCCESS); 717 } 718 719 static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A) 720 { 721 PetscSection section; 722 PetscScalar *values; 723 PetscInt rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0; 724 725 PetscFunctionBegin; 726 PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd)); 727 for (r = rStart; r < rEnd; ++r) { 728 PetscCall(PetscSectionGetDof(sectionAdj, r, &len)); 729 maxRowLen = PetscMax(maxRowLen, len); 730 } 731 PetscCall(PetscCalloc1(maxRowLen, &values)); 732 if (f >= 0 && bs == 1) { 733 PetscCall(DMGetLocalSection(dm, §ion)); 734 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 735 for (p = pStart; p < pEnd; ++p) { 736 PetscInt rS, rE; 737 738 PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE)); 739 for (r = rS; r < rE; ++r) { 740 PetscInt numCols, cStart; 741 742 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 743 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 744 PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES)); 745 } 746 } 747 } else { 748 for (r = rStart; r < rEnd; ++r) { 749 PetscInt numCols, cStart; 750 751 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 752 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 753 PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES)); 754 } 755 } 756 PetscCall(PetscFree(values)); 757 PetscFunctionReturn(PETSC_SUCCESS); 758 } 759 760 /*@ 761 DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the `DM`, 762 the `PetscDS` it contains, and the default `PetscSection`. 763 764 Collective 765 766 Input Parameters: 767 + dm - The `DMPLEX` 768 . bs - The matrix blocksize 769 . dnz - An array to hold the number of nonzeros in the diagonal block 770 . onz - An array to hold the number of nonzeros in the off-diagonal block 771 . dnzu - An array to hold the number of nonzeros in the upper triangle of the diagonal block 772 . onzu - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block 773 - fillMatrix - If `PETSC_TRUE`, fill the matrix with zeros 774 775 Output Parameter: 776 . A - The preallocated matrix 777 778 Level: advanced 779 780 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreateMatrix()` 781 @*/ 782 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix) 783 { 784 MPI_Comm comm; 785 MatType mtype; 786 PetscSF sf, sfDof; 787 PetscSection section; 788 PetscInt *remoteOffsets; 789 PetscSection sectionAdj[4] = {NULL, NULL, NULL, NULL}; 790 PetscInt *cols[4] = {NULL, NULL, NULL, NULL}; 791 PetscBool useCone, useClosure; 792 PetscInt Nf, f, idx, locRows; 793 PetscLayout rLayout; 794 PetscBool isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE; 795 796 PetscFunctionBegin; 797 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 798 PetscValidHeaderSpecific(A, MAT_CLASSID, 7); 799 if (dnz) PetscAssertPointer(dnz, 3); 800 if (onz) PetscAssertPointer(onz, 4); 801 if (dnzu) PetscAssertPointer(dnzu, 5); 802 if (onzu) PetscAssertPointer(onzu, 6); 803 PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf)); 804 PetscCall(DMGetLocalSection(dm, §ion)); 805 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL)); 806 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 807 PetscCall(PetscLogEventBegin(DMPLEX_Preallocate, dm, 0, 0, 0)); 808 /* Create dof SF based on point SF */ 809 if (debug) { 810 PetscSection section, sectionGlobal; 811 PetscSF sf; 812 813 PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf)); 814 PetscCall(DMGetLocalSection(dm, §ion)); 815 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 816 PetscCall(PetscPrintf(comm, "Input Section for Preallocation:\n")); 817 PetscCall(PetscSectionView(section, NULL)); 818 PetscCall(PetscPrintf(comm, "Input Global Section for Preallocation:\n")); 819 PetscCall(PetscSectionView(sectionGlobal, NULL)); 820 PetscCall(PetscPrintf(comm, "Input SF for Preallocation:\n")); 821 PetscCall(PetscSFView(sf, NULL)); 822 } 823 PetscCall(PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets)); 824 PetscCall(PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof)); 825 PetscCall(PetscFree(remoteOffsets)); 826 if (debug) { 827 PetscCall(PetscPrintf(comm, "Dof SF for Preallocation:\n")); 828 PetscCall(PetscSFView(sfDof, NULL)); 829 } 830 /* Create allocation vectors from adjacency graph */ 831 PetscCall(MatGetLocalSize(A, &locRows, NULL)); 832 PetscCall(PetscLayoutCreate(comm, &rLayout)); 833 PetscCall(PetscLayoutSetLocalSize(rLayout, locRows)); 834 PetscCall(PetscLayoutSetBlockSize(rLayout, 1)); 835 PetscCall(PetscLayoutSetUp(rLayout)); 836 /* There are 4 types of adjacency */ 837 PetscCall(PetscSectionGetNumFields(section, &Nf)); 838 if (Nf < 1 || bs > 1) { 839 PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 840 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 841 PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx])); 842 PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu)); 843 } else { 844 for (f = 0; f < Nf; ++f) { 845 PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure)); 846 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 847 if (!sectionAdj[idx]) PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx])); 848 PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu)); 849 } 850 } 851 PetscCall(PetscSFDestroy(&sfDof)); 852 /* Set matrix pattern */ 853 PetscCall(MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu)); 854 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 855 /* Check for symmetric storage */ 856 PetscCall(MatGetType(A, &mtype)); 857 PetscCall(PetscStrcmp(mtype, MATSBAIJ, &isSymBlock)); 858 PetscCall(PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock)); 859 PetscCall(PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock)); 860 if (isSymBlock || isSymSeqBlock || isSymMPIBlock) PetscCall(MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE)); 861 /* Fill matrix with zeros */ 862 if (fillMatrix) { 863 if (Nf < 1 || bs > 1) { 864 PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 865 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 866 PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A)); 867 } else { 868 for (f = 0; f < Nf; ++f) { 869 PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure)); 870 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 871 PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A)); 872 } 873 } 874 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 875 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 876 } 877 PetscCall(PetscLayoutDestroy(&rLayout)); 878 for (idx = 0; idx < 4; ++idx) { 879 PetscCall(PetscSectionDestroy(§ionAdj[idx])); 880 PetscCall(PetscFree(cols[idx])); 881 } 882 PetscCall(PetscLogEventEnd(DMPLEX_Preallocate, dm, 0, 0, 0)); 883 PetscFunctionReturn(PETSC_SUCCESS); 884 } 885 886 #if 0 887 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix) 888 { 889 PetscInt *tmpClosure,*tmpAdj,*visits; 890 PetscInt c,cStart,cEnd,pStart,pEnd; 891 892 PetscFunctionBegin; 893 PetscCall(DMGetDimension(dm, &dim)); 894 PetscCall(DMPlexGetDepth(dm, &depth)); 895 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize)); 896 897 maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)); 898 899 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 900 npoints = pEnd - pStart; 901 902 PetscCall(PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits)); 903 PetscCall(PetscArrayzero(lvisits,pEnd-pStart)); 904 PetscCall(PetscArrayzero(visits,pEnd-pStart)); 905 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 906 for (c=cStart; c<cEnd; c++) { 907 PetscInt *support = tmpClosure; 908 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support)); 909 for (p=0; p<supportSize; p++) lvisits[support[p]]++; 910 } 911 PetscCall(PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM)); 912 PetscCall(PetscSFReduceEnd (sf,MPIU_INT,lvisits,visits,MPI_SUM)); 913 PetscCall(PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits,MPI_REPLACE)); 914 PetscCall(PetscSFBcastEnd (sf,MPIU_INT,visits,lvisits)); 915 916 PetscCall(PetscSFGetRootRanks()); 917 918 PetscCall(PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner)); 919 for (c=cStart; c<cEnd; c++) { 920 PetscCall(PetscArrayzero(cellmat,maxClosureSize*maxClosureSize)); 921 /* 922 Depth-first walk of transitive closure. 923 At each leaf frame f of transitive closure that we see, add 1/visits[f] to each pair (p,q) not marked as done in cellmat. 924 This contribution is added to dnz if owning ranks of p and q match, to onz otherwise. 925 */ 926 } 927 928 PetscCall(PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM)); 929 PetscCall(PetscSFReduceEnd (sf,MPIU_INT,lonz,onz,MPI_SUM)); 930 PetscFunctionReturn(PETSC_SUCCESS); 931 } 932 #endif 933