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(PetscSection myRankPairSection, 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, num_leaves; 208 209 PetscFunctionBeginUser; 210 PetscCall(PetscSectionGetChart(myRankPairSection, NULL, &num_roots)); 211 PetscCall(PetscSectionGetStorageSize(myRankPairSection, &num_leaves)); 212 *num_leaves_found = 0; 213 for (PetscInt q = 0; q < numAdj; q++) { 214 const PetscInt padj = tmpAdj[q]; 215 PetscCall(PetscFindInt(padj, num_roots, roots, &root_index)); 216 if (root_index >= 0) { 217 PetscInt dof, offset; 218 219 PetscCall(PetscSectionGetDof(myRankPairSection, root_index, &dof)); 220 PetscCall(PetscSectionGetOffset(myRankPairSection, root_index, &offset)); 221 222 for (PetscInt l = 0; l < dof; l++) { 223 leaf_ = leaves[offset + l]; 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 break; 229 } 230 } 231 } 232 } 233 } 234 235 PetscCall(PetscIntSortSemiOrdered(*num_leaves_found, leaves_found)); 236 PetscFunctionReturn(PETSC_SUCCESS); 237 } 238 239 static PetscErrorCode DMPlexCreateAdjacencySection_Static(DM dm, PetscInt bs, PetscSF sfDof, PetscBool useCone, PetscBool useClosure, PetscBool useAnchors, PetscSection *sA, PetscInt **colIdx) 240 { 241 MPI_Comm comm; 242 PetscMPIInt myrank; 243 PetscBool doCommLocal, doComm, debug = PETSC_FALSE; 244 PetscSF sf, sfAdj; 245 PetscSection section, sectionGlobal, leafSectionAdj, rootSectionAdj, sectionAdj, anchorSectionAdj, myRankPairSection; 246 PetscInt nroots, nleaves, l, p, r; 247 const PetscInt *leaves; 248 const PetscSFNode *remotes; 249 PetscInt dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols, adjSize; 250 PetscInt *tmpAdj = NULL, *adj, *rootAdj, *anchorAdj = NULL, *cols, *remoteOffsets, *myRankPairRoots = NULL, *myRankPairLeaves = NULL; 251 252 PetscFunctionBegin; 253 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 254 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL)); 255 PetscCallMPI(MPI_Comm_rank(comm, &myrank)); 256 PetscCall(DMGetDimension(dm, &dim)); 257 PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf)); 258 PetscCall(DMGetLocalSection(dm, §ion)); 259 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 260 PetscCall(PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL)); 261 doCommLocal = nroots >= 0 ? PETSC_TRUE : PETSC_FALSE; 262 PetscCallMPI(MPIU_Allreduce(&doCommLocal, &doComm, 1, MPI_C_BOOL, MPI_LAND, comm)); 263 /* Create section for dof adjacency (dof ==> # adj dof) */ 264 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 265 PetscCall(PetscSectionGetStorageSize(section, &numDof)); 266 PetscCall(PetscSectionCreate(comm, &leafSectionAdj)); 267 PetscCall(PetscSectionSetChart(leafSectionAdj, 0, numDof)); 268 PetscCall(PetscSectionCreate(comm, &rootSectionAdj)); 269 PetscCall(PetscSectionSetChart(rootSectionAdj, 0, numDof)); 270 PetscCall(PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes)); 271 272 // Store leaf-root pairs if the leaf and root are both located on current rank. 273 // The point in myRankPairSection is an index into myRankPairRoots. 274 // The section then defines the number of pairs for that root and the offset into myRankPairLeaves to access them. 275 { 276 PetscSegBuffer seg_roots, seg_leaves; 277 PetscCount buffer_size; 278 PetscInt *roots_with_dups, num_non_dups, num_pairs = 0; 279 280 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 32, &seg_roots)); 281 PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 32, &seg_leaves)); 282 for (PetscInt l = 0; l < nleaves; l++) { 283 if (remotes[l].rank == myrank) { 284 PetscInt *slot; 285 PetscCall(PetscSegBufferGetInts(seg_roots, 1, &slot)); 286 *slot = remotes[l].index; 287 PetscCall(PetscSegBufferGetInts(seg_leaves, 1, &slot)); 288 *slot = leaves[l]; 289 } 290 } 291 PetscCall(PetscSegBufferGetSize(seg_roots, &buffer_size)); 292 PetscCall(PetscIntCast(buffer_size, &num_pairs)); 293 PetscCall(PetscSegBufferExtractAlloc(seg_roots, &roots_with_dups)); 294 PetscCall(PetscSegBufferExtractAlloc(seg_leaves, &myRankPairLeaves)); 295 PetscCall(PetscSegBufferDestroy(&seg_roots)); 296 PetscCall(PetscSegBufferDestroy(&seg_leaves)); 297 298 PetscCall(PetscIntSortSemiOrderedWithArray(num_pairs, roots_with_dups, myRankPairLeaves)); 299 if (debug) { 300 PetscCall(PetscPrintf(comm, "Root/leaf pairs on the same rank:\n")); 301 PetscCall(PetscIntViewPairs(num_pairs, 5, roots_with_dups, myRankPairLeaves, NULL)); 302 } 303 PetscCall(PetscMalloc1(num_pairs, &myRankPairRoots)); 304 PetscCall(PetscArraycpy(myRankPairRoots, roots_with_dups, num_pairs)); 305 306 num_non_dups = num_pairs; 307 PetscCall(PetscSortedRemoveDupsInt(&num_non_dups, myRankPairRoots)); 308 PetscCall(PetscSectionCreate(comm, &myRankPairSection)); 309 PetscCall(PetscSectionSetChart(myRankPairSection, 0, num_non_dups)); 310 for (PetscInt p = 0; p < num_pairs; p++) { 311 PetscInt root = roots_with_dups[p], index; 312 PetscCall(PetscFindInt(root, num_non_dups, myRankPairRoots, &index)); 313 PetscAssert(index >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Root not found after removal of duplicates"); 314 PetscCall(PetscSectionAddDof(myRankPairSection, index, 1)); 315 } 316 PetscCall(PetscSectionSetUp(myRankPairSection)); 317 318 if (debug) { 319 PetscCall(PetscPrintf(comm, "Root/leaf pair section on the same rank:\n")); 320 PetscCall(PetscIntView(num_non_dups, myRankPairRoots, NULL)); 321 PetscCall(PetscSectionArrayView(myRankPairSection, myRankPairLeaves, PETSC_INT, NULL)); 322 } 323 PetscCall(PetscFree(roots_with_dups)); 324 } 325 326 /* 327 section - maps points to (# dofs, local dofs) 328 sectionGlobal - maps points to (# dofs, global dofs) 329 leafSectionAdj - maps unowned local dofs to # adj dofs 330 rootSectionAdj - maps owned local dofs to # adj dofs 331 adj - adj global dofs indexed by leafSectionAdj 332 rootAdj - adj global dofs indexed by rootSectionAdj 333 sf - describes shared points across procs 334 sfDof - describes shared dofs across procs 335 sfAdj - describes shared adjacent dofs across procs 336 ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point. 337 (0). If there are point-to-point constraints, add the adjacencies of constrained points to anchors in anchorAdj 338 (This is done in DMPlexComputeAnchorAdjacencies()) 339 1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj 340 Reduce those counts to rootSectionAdj (now redundantly counting some interface points) 341 2. Visit owned points on interface, count adjacencies placing in rootSectionAdj 342 Create sfAdj connecting rootSectionAdj and leafSectionAdj 343 3. Visit unowned points on interface, write adjacencies to adj 344 Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies) 345 4. Visit owned points on interface, write adjacencies to rootAdj 346 Remove redundancy in rootAdj 347 ** The last two traversals use transitive closure 348 5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj) 349 Allocate memory addressed by sectionAdj (cols) 350 6. Visit all owned points in the subdomain, insert dof adjacencies into cols 351 ** Knowing all the column adjacencies, check ownership and sum into dnz and onz 352 */ 353 PetscCall(DMPlexComputeAnchorAdjacencies(dm, useCone, useClosure, &anchorSectionAdj, &anchorAdj)); 354 for (l = 0; l < nleaves; ++l) { 355 PetscInt dof, off, d, q, anDof; 356 PetscInt p = leaves[l], numAdj = PETSC_DETERMINE; 357 358 if ((p < pStart) || (p >= pEnd)) continue; 359 PetscCall(PetscSectionGetDof(section, p, &dof)); 360 PetscCall(PetscSectionGetOffset(section, p, &off)); 361 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 362 for (q = 0; q < numAdj; ++q) { 363 const PetscInt padj = tmpAdj[q]; 364 PetscInt ndof, ncdof; 365 366 if ((padj < pStart) || (padj >= pEnd)) continue; 367 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 368 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 369 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, ndof - ncdof)); 370 } 371 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 372 if (anDof) { 373 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(leafSectionAdj, d, anDof)); 374 } 375 } 376 PetscCall(PetscSectionSetUp(leafSectionAdj)); 377 if (debug) { 378 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n")); 379 PetscCall(PetscSectionView(leafSectionAdj, NULL)); 380 } 381 /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */ 382 if (doComm) { 383 PetscCall(PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM)); 384 PetscCall(PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM)); 385 PetscCall(PetscSectionInvalidateMaxDof_Internal(rootSectionAdj)); 386 } 387 if (debug) { 388 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots:\n")); 389 PetscCall(PetscSectionView(rootSectionAdj, NULL)); 390 } 391 /* Add in local adjacency sizes for owned dofs on interface (roots) */ 392 for (p = pStart; p < pEnd; ++p) { 393 PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof; 394 395 PetscCall(PetscSectionGetDof(section, p, &dof)); 396 PetscCall(PetscSectionGetOffset(section, p, &off)); 397 if (!dof) continue; 398 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 399 if (adof <= 0) continue; 400 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 401 for (q = 0; q < numAdj; ++q) { 402 const PetscInt padj = tmpAdj[q]; 403 PetscInt ndof, ncdof; 404 405 if ((padj < pStart) || (padj >= pEnd)) continue; 406 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 407 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 408 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, ndof - ncdof)); 409 } 410 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 411 if (anDof) { 412 for (d = off; d < off + dof; ++d) PetscCall(PetscSectionAddDof(rootSectionAdj, d, anDof)); 413 } 414 } 415 PetscCall(PetscSectionSetUp(rootSectionAdj)); 416 if (debug) { 417 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after local additions:\n")); 418 PetscCall(PetscSectionView(rootSectionAdj, NULL)); 419 } 420 /* Create adj SF based on dof SF */ 421 PetscCall(PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets)); 422 PetscCall(PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj)); 423 PetscCall(PetscFree(remoteOffsets)); 424 if (debug) { 425 PetscCall(PetscPrintf(comm, "Adjacency SF for Preallocation:\n")); 426 PetscCall(PetscSFView(sfAdj, NULL)); 427 } 428 /* Create leaf adjacency */ 429 PetscCall(PetscSectionSetUp(leafSectionAdj)); 430 PetscCall(PetscSectionGetStorageSize(leafSectionAdj, &adjSize)); 431 PetscCall(PetscCalloc1(adjSize, &adj)); 432 for (l = 0; l < nleaves; ++l) { 433 PetscInt dof, off, d, q, anDof, anOff; 434 PetscInt p = leaves[l], numAdj = PETSC_DETERMINE; 435 436 if ((p < pStart) || (p >= pEnd)) continue; 437 PetscCall(PetscSectionGetDof(section, p, &dof)); 438 PetscCall(PetscSectionGetOffset(section, p, &off)); 439 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 440 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 441 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 442 for (d = off; d < off + dof; ++d) { 443 PetscInt aoff, i = 0; 444 445 PetscCall(PetscSectionGetOffset(leafSectionAdj, d, &aoff)); 446 for (q = 0; q < numAdj; ++q) { 447 const PetscInt padj = tmpAdj[q]; 448 PetscInt ndof, ncdof, ngoff, nd; 449 450 if ((padj < pStart) || (padj >= pEnd)) continue; 451 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 452 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 453 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 454 for (nd = 0; nd < ndof - ncdof; ++nd) { 455 adj[aoff + i] = (ngoff < 0 ? -(ngoff + 1) : ngoff) + nd; 456 ++i; 457 } 458 } 459 for (q = 0; q < anDof; q++) { 460 adj[aoff + i] = anchorAdj[anOff + q]; 461 ++i; 462 } 463 } 464 } 465 /* Debugging */ 466 if (debug) { 467 PetscCall(PetscPrintf(comm, "Leaf adjacency indices\n")); 468 PetscCall(PetscSectionArrayView(leafSectionAdj, adj, PETSC_INT, NULL)); 469 } 470 /* Gather adjacent indices to root */ 471 PetscCall(PetscSectionGetStorageSize(rootSectionAdj, &adjSize)); 472 PetscCall(PetscMalloc1(adjSize, &rootAdj)); 473 for (r = 0; r < adjSize; ++r) rootAdj[r] = -1; 474 if (doComm) { 475 const PetscInt *indegree; 476 PetscInt *remoteadj, radjsize = 0; 477 478 PetscCall(PetscSFComputeDegreeBegin(sfAdj, &indegree)); 479 PetscCall(PetscSFComputeDegreeEnd(sfAdj, &indegree)); 480 for (p = 0; p < adjSize; ++p) radjsize += indegree[p]; 481 PetscCall(PetscMalloc1(radjsize, &remoteadj)); 482 PetscCall(PetscSFGatherBegin(sfAdj, MPIU_INT, adj, remoteadj)); 483 PetscCall(PetscSFGatherEnd(sfAdj, MPIU_INT, adj, remoteadj)); 484 for (p = 0, l = 0, r = 0; p < adjSize; ++p, l = PetscMax(p, l + indegree[p - 1])) { 485 PetscInt s; 486 for (s = 0; s < indegree[p]; ++s, ++r) rootAdj[l + s] = remoteadj[r]; 487 } 488 PetscCheck(r == radjsize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, r, radjsize); 489 PetscCheck(l == adjSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistency in communication %" PetscInt_FMT " != %" PetscInt_FMT, l, adjSize); 490 PetscCall(PetscFree(remoteadj)); 491 } 492 PetscCall(PetscSFDestroy(&sfAdj)); 493 PetscCall(PetscFree(adj)); 494 /* Debugging */ 495 if (debug) { 496 PetscCall(PetscPrintf(comm, "Root adjacency indices after gather\n")); 497 PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL)); 498 } 499 /* Add in local adjacency indices for owned dofs on interface (roots) */ 500 for (p = pStart; p < pEnd; ++p) { 501 PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q, anDof, anOff; 502 503 PetscCall(PetscSectionGetDof(section, p, &dof)); 504 PetscCall(PetscSectionGetOffset(section, p, &off)); 505 if (!dof) continue; 506 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 507 if (adof <= 0) continue; 508 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 509 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 510 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 511 for (d = off; d < off + dof; ++d) { 512 PetscInt adof, aoff, i; 513 514 PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof)); 515 PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff)); 516 i = adof - 1; 517 for (q = 0; q < anDof; q++) { 518 rootAdj[aoff + i] = anchorAdj[anOff + q]; 519 --i; 520 } 521 for (q = 0; q < numAdj; ++q) { 522 const PetscInt padj = tmpAdj[q]; 523 PetscInt ndof, ncdof, ngoff, nd; 524 525 if ((padj < pStart) || (padj >= pEnd)) continue; 526 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 527 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 528 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 529 for (nd = 0; nd < ndof - ncdof; ++nd) { 530 rootAdj[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd; 531 --i; 532 } 533 } 534 } 535 } 536 /* Debugging */ 537 if (debug) { 538 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots before compression:\n")); 539 PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL)); 540 } 541 /* Compress indices */ 542 PetscCall(PetscSectionSetUp(rootSectionAdj)); 543 for (p = pStart; p < pEnd; ++p) { 544 PetscInt dof, cdof, off, d; 545 PetscInt adof, aoff; 546 547 PetscCall(PetscSectionGetDof(section, p, &dof)); 548 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 549 PetscCall(PetscSectionGetOffset(section, p, &off)); 550 if (!dof) continue; 551 PetscCall(PetscSectionGetDof(rootSectionAdj, off, &adof)); 552 if (adof <= 0) continue; 553 for (d = off; d < off + dof - cdof; ++d) { 554 PetscCall(PetscSectionGetDof(rootSectionAdj, d, &adof)); 555 PetscCall(PetscSectionGetOffset(rootSectionAdj, d, &aoff)); 556 PetscCall(PetscSortRemoveDupsInt(&adof, &rootAdj[aoff])); 557 PetscCall(PetscSectionSetDof(rootSectionAdj, d, adof)); 558 } 559 } 560 /* Debugging */ 561 if (debug) { 562 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation on Roots after compression:\n")); 563 PetscCall(PetscSectionArrayView(rootSectionAdj, rootAdj, PETSC_INT, NULL)); 564 } 565 /* Build adjacency section: Maps global indices to sets of adjacent global indices */ 566 PetscCall(PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd)); 567 PetscCall(PetscSectionCreate(comm, §ionAdj)); 568 PetscCall(PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd)); 569 570 PetscInt *exclude_leaves, num_exclude_leaves = 0, max_adjacency_size = 0; 571 PetscCall(DMPlexGetMaxAdjacencySize_Internal(dm, useAnchors, &max_adjacency_size)); 572 PetscCall(PetscMalloc1(max_adjacency_size, &exclude_leaves)); 573 574 for (p = pStart; p < pEnd; ++p) { 575 PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof; 576 PetscBool found = PETSC_TRUE; 577 578 PetscCall(PetscSectionGetDof(section, p, &dof)); 579 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 580 PetscCall(PetscSectionGetOffset(section, p, &off)); 581 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 582 for (d = 0; d < dof - cdof; ++d) { 583 PetscInt ldof, rdof; 584 585 PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof)); 586 PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof)); 587 if (ldof > 0) { 588 /* We do not own this point */ 589 } else if (rdof > 0) { 590 PetscCall(PetscSectionSetDof(sectionAdj, goff + d, rdof)); 591 } else { 592 found = PETSC_FALSE; 593 } 594 } 595 if (found) continue; 596 PetscCall(PetscSectionGetDof(section, p, &dof)); 597 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 598 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 599 PetscCall(AdjancencyContainsLeafRootPair(myRankPairSection, myRankPairLeaves, myRankPairRoots, numAdj, tmpAdj, &num_exclude_leaves, exclude_leaves)); 600 for (q = 0; q < numAdj; ++q) { 601 const PetscInt padj = tmpAdj[q]; 602 PetscInt ndof, ncdof, noff, count; 603 604 if ((padj < pStart) || (padj >= pEnd)) continue; 605 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 606 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 607 PetscCall(PetscSectionGetOffset(section, padj, &noff)); 608 // If leaf-root pair are both on this rank, only count root 609 PetscCall(PetscFindInt(padj, num_exclude_leaves, exclude_leaves, &count)); 610 if (count >= 0) continue; 611 for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, ndof - ncdof)); 612 } 613 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 614 if (anDof) { 615 for (d = goff; d < goff + dof - cdof; ++d) PetscCall(PetscSectionAddDof(sectionAdj, d, anDof)); 616 } 617 } 618 PetscCall(PetscSectionSetUp(sectionAdj)); 619 if (debug) { 620 PetscCall(PetscPrintf(comm, "Adjacency Section for Preallocation:\n")); 621 PetscCall(PetscSectionView(sectionAdj, NULL)); 622 } 623 /* Get adjacent indices */ 624 PetscCall(PetscSectionGetStorageSize(sectionAdj, &numCols)); 625 PetscCall(PetscMalloc1(numCols, &cols)); 626 for (p = pStart; p < pEnd; ++p) { 627 PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q, anDof, anOff; 628 PetscBool found = PETSC_TRUE; 629 630 PetscCall(PetscSectionGetDof(section, p, &dof)); 631 PetscCall(PetscSectionGetConstraintDof(section, p, &cdof)); 632 PetscCall(PetscSectionGetOffset(section, p, &off)); 633 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &goff)); 634 for (d = 0; d < dof - cdof; ++d) { 635 PetscInt ldof, rdof; 636 637 PetscCall(PetscSectionGetDof(leafSectionAdj, off + d, &ldof)); 638 PetscCall(PetscSectionGetDof(rootSectionAdj, off + d, &rdof)); 639 if (ldof > 0) { 640 /* We do not own this point */ 641 } else if (rdof > 0) { 642 PetscInt aoff, roff; 643 644 PetscCall(PetscSectionGetOffset(sectionAdj, goff + d, &aoff)); 645 PetscCall(PetscSectionGetOffset(rootSectionAdj, off + d, &roff)); 646 PetscCall(PetscArraycpy(&cols[aoff], &rootAdj[roff], rdof)); 647 } else { 648 found = PETSC_FALSE; 649 } 650 } 651 if (found) continue; 652 PetscCall(DMPlexGetAdjacency_Internal(dm, p, useCone, useClosure, useAnchors, &numAdj, &tmpAdj)); 653 PetscCall(AdjancencyContainsLeafRootPair(myRankPairSection, myRankPairLeaves, myRankPairRoots, numAdj, tmpAdj, &num_exclude_leaves, exclude_leaves)); 654 PetscCall(PetscSectionGetDof(anchorSectionAdj, p, &anDof)); 655 PetscCall(PetscSectionGetOffset(anchorSectionAdj, p, &anOff)); 656 for (d = goff; d < goff + dof - cdof; ++d) { 657 PetscInt adof, aoff, i = 0; 658 659 PetscCall(PetscSectionGetDof(sectionAdj, d, &adof)); 660 PetscCall(PetscSectionGetOffset(sectionAdj, d, &aoff)); 661 for (q = 0; q < numAdj; ++q) { 662 const PetscInt padj = tmpAdj[q], *ncind; 663 PetscInt ndof, ncdof, ngoff, nd, count; 664 665 /* Adjacent points may not be in the section chart */ 666 if ((padj < pStart) || (padj >= pEnd)) continue; 667 PetscCall(PetscSectionGetDof(section, padj, &ndof)); 668 PetscCall(PetscSectionGetConstraintDof(section, padj, &ncdof)); 669 PetscCall(PetscSectionGetConstraintIndices(section, padj, &ncind)); 670 PetscCall(PetscSectionGetOffset(sectionGlobal, padj, &ngoff)); 671 // If leaf-root pair are both on this rank, only count root 672 PetscCall(PetscFindInt(padj, num_exclude_leaves, exclude_leaves, &count)); 673 if (count >= 0) continue; 674 for (nd = 0; nd < ndof - ncdof; ++nd, ++i) cols[aoff + i] = ngoff < 0 ? -(ngoff + 1) + nd : ngoff + nd; 675 } 676 for (q = 0; q < anDof; q++, i++) cols[aoff + i] = anchorAdj[anOff + q]; 677 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); 678 } 679 } 680 PetscCall(PetscSectionDestroy(&anchorSectionAdj)); 681 PetscCall(PetscSectionDestroy(&leafSectionAdj)); 682 PetscCall(PetscSectionDestroy(&rootSectionAdj)); 683 PetscCall(PetscSectionDestroy(&myRankPairSection)); 684 PetscCall(PetscFree(exclude_leaves)); 685 PetscCall(PetscFree(myRankPairLeaves)); 686 PetscCall(PetscFree(myRankPairRoots)); 687 PetscCall(PetscFree(anchorAdj)); 688 PetscCall(PetscFree(rootAdj)); 689 PetscCall(PetscFree(tmpAdj)); 690 /* Debugging */ 691 if (debug) { 692 PetscCall(PetscPrintf(comm, "Column indices\n")); 693 PetscCall(PetscSectionArrayView(sectionAdj, cols, PETSC_INT, NULL)); 694 } 695 696 *sA = sectionAdj; 697 *colIdx = cols; 698 PetscFunctionReturn(PETSC_SUCCESS); 699 } 700 701 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[]) 702 { 703 PetscSection section; 704 PetscInt rStart, rEnd, r, pStart, pEnd, p; 705 706 PetscFunctionBegin; 707 /* This loop needs to change to a loop over points, then field dofs, which means we need to look both sections */ 708 PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd)); 709 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); 710 if (f >= 0 && bs == 1) { 711 PetscCall(DMGetLocalSection(dm, §ion)); 712 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 713 for (p = pStart; p < pEnd; ++p) { 714 PetscInt rS, rE; 715 716 PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE)); 717 for (r = rS; r < rE; ++r) { 718 PetscInt numCols, cStart, c; 719 720 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 721 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 722 for (c = cStart; c < cStart + numCols; ++c) { 723 if ((cols[c] >= rStart) && (cols[c] < rEnd)) { 724 ++dnz[r - rStart]; 725 if (cols[c] >= r) ++dnzu[r - rStart]; 726 } else { 727 ++onz[r - rStart]; 728 if (cols[c] >= r) ++onzu[r - rStart]; 729 } 730 } 731 } 732 } 733 } else { 734 /* Only loop over blocks of rows */ 735 for (r = rStart / bs; r < rEnd / bs; ++r) { 736 const PetscInt row = r * bs; 737 PetscInt numCols, cStart, c; 738 739 PetscCall(PetscSectionGetDof(sectionAdj, row, &numCols)); 740 PetscCall(PetscSectionGetOffset(sectionAdj, row, &cStart)); 741 for (c = cStart; c < cStart + numCols; ++c) { 742 if ((cols[c] >= rStart) && (cols[c] < rEnd)) { 743 ++dnz[r - rStart / bs]; 744 if (cols[c] >= row) ++dnzu[r - rStart / bs]; 745 } else { 746 ++onz[r - rStart / bs]; 747 if (cols[c] >= row) ++onzu[r - rStart / bs]; 748 } 749 } 750 } 751 for (r = 0; r < (rEnd - rStart) / bs; ++r) { 752 dnz[r] /= bs; 753 onz[r] /= bs; 754 dnzu[r] /= bs; 755 onzu[r] /= bs; 756 } 757 } 758 PetscFunctionReturn(PETSC_SUCCESS); 759 } 760 761 static PetscErrorCode DMPlexFillMatrix_Static(DM dm, PetscLayout rLayout, PetscInt bs, PetscInt f, PetscSection sectionAdj, const PetscInt cols[], Mat A) 762 { 763 PetscSection section; 764 PetscScalar *values; 765 PetscInt rStart, rEnd, r, pStart, pEnd, p, len, maxRowLen = 0; 766 767 PetscFunctionBegin; 768 PetscCall(PetscLayoutGetRange(rLayout, &rStart, &rEnd)); 769 for (r = rStart; r < rEnd; ++r) { 770 PetscCall(PetscSectionGetDof(sectionAdj, r, &len)); 771 maxRowLen = PetscMax(maxRowLen, len); 772 } 773 PetscCall(PetscCalloc1(maxRowLen, &values)); 774 if (f >= 0 && bs == 1) { 775 PetscCall(DMGetLocalSection(dm, §ion)); 776 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 777 for (p = pStart; p < pEnd; ++p) { 778 PetscInt rS, rE; 779 780 PetscCall(DMGetGlobalFieldOffset_Private(dm, p, f, &rS, &rE)); 781 for (r = rS; r < rE; ++r) { 782 PetscInt numCols, cStart; 783 784 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 785 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 786 PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES)); 787 } 788 } 789 } else { 790 for (r = rStart; r < rEnd; ++r) { 791 PetscInt numCols, cStart; 792 793 PetscCall(PetscSectionGetDof(sectionAdj, r, &numCols)); 794 PetscCall(PetscSectionGetOffset(sectionAdj, r, &cStart)); 795 PetscCall(MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES)); 796 } 797 } 798 PetscCall(PetscFree(values)); 799 PetscFunctionReturn(PETSC_SUCCESS); 800 } 801 802 /*@ 803 DMPlexPreallocateOperator - Calculate the matrix nonzero pattern based upon the information in the `DM`, 804 the `PetscDS` it contains, and the default `PetscSection`. 805 806 Collective 807 808 Input Parameters: 809 + dm - The `DMPLEX` 810 . bs - The matrix blocksize 811 . dnz - An array to hold the number of nonzeros in the diagonal block 812 . onz - An array to hold the number of nonzeros in the off-diagonal block 813 . dnzu - An array to hold the number of nonzeros in the upper triangle of the diagonal block 814 . onzu - An array to hold the number of nonzeros in the upper triangle of the off-diagonal block 815 - fillMatrix - If `PETSC_TRUE`, fill the matrix with zeros 816 817 Output Parameter: 818 . A - The preallocated matrix 819 820 Level: advanced 821 822 .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreateMatrix()` 823 @*/ 824 PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix) 825 { 826 MPI_Comm comm; 827 MatType mtype; 828 PetscSF sf, sfDof; 829 PetscSection section; 830 PetscInt *remoteOffsets; 831 PetscSection sectionAdj[4] = {NULL, NULL, NULL, NULL}; 832 PetscInt *cols[4] = {NULL, NULL, NULL, NULL}; 833 PetscBool useCone, useClosure; 834 PetscInt Nf, f, idx, locRows; 835 PetscLayout rLayout; 836 PetscBool isSymBlock, isSymSeqBlock, isSymMPIBlock, debug = PETSC_FALSE; 837 838 PetscFunctionBegin; 839 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 840 PetscValidHeaderSpecific(A, MAT_CLASSID, 7); 841 if (dnz) PetscAssertPointer(dnz, 3); 842 if (onz) PetscAssertPointer(onz, 4); 843 if (dnzu) PetscAssertPointer(dnzu, 5); 844 if (onzu) PetscAssertPointer(onzu, 6); 845 PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf)); 846 PetscCall(DMGetLocalSection(dm, §ion)); 847 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_view_preallocation", &debug, NULL)); 848 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 849 PetscCall(PetscLogEventBegin(DMPLEX_Preallocate, dm, 0, 0, 0)); 850 /* Create dof SF based on point SF */ 851 if (debug) { 852 PetscSection section, sectionGlobal; 853 PetscSF sf; 854 855 PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &sf)); 856 PetscCall(DMGetLocalSection(dm, §ion)); 857 PetscCall(DMGetGlobalSection(dm, §ionGlobal)); 858 PetscCall(PetscPrintf(comm, "Input Section for Preallocation:\n")); 859 PetscCall(PetscSectionView(section, NULL)); 860 PetscCall(PetscPrintf(comm, "Input Global Section for Preallocation:\n")); 861 PetscCall(PetscSectionView(sectionGlobal, NULL)); 862 PetscCall(PetscPrintf(comm, "Input SF for Preallocation:\n")); 863 PetscCall(PetscSFView(sf, NULL)); 864 } 865 PetscCall(PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets)); 866 PetscCall(PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof)); 867 PetscCall(PetscFree(remoteOffsets)); 868 if (debug) { 869 PetscCall(PetscPrintf(comm, "Dof SF for Preallocation:\n")); 870 PetscCall(PetscSFView(sfDof, NULL)); 871 } 872 /* Create allocation vectors from adjacency graph */ 873 PetscCall(MatGetLocalSize(A, &locRows, NULL)); 874 PetscCall(PetscLayoutCreate(comm, &rLayout)); 875 PetscCall(PetscLayoutSetLocalSize(rLayout, locRows)); 876 PetscCall(PetscLayoutSetBlockSize(rLayout, 1)); 877 PetscCall(PetscLayoutSetUp(rLayout)); 878 /* There are 4 types of adjacency */ 879 PetscCall(PetscSectionGetNumFields(section, &Nf)); 880 if (Nf < 1 || bs > 1) { 881 PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 882 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 883 PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx])); 884 PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu)); 885 } else { 886 for (f = 0; f < Nf; ++f) { 887 PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure)); 888 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 889 if (!sectionAdj[idx]) PetscCall(DMPlexCreateAdjacencySection_Static(dm, bs, sfDof, useCone, useClosure, PETSC_TRUE, §ionAdj[idx], &cols[idx])); 890 PetscCall(DMPlexUpdateAllocation_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], dnz, onz, dnzu, onzu)); 891 } 892 } 893 PetscCall(PetscSFDestroy(&sfDof)); 894 /* Set matrix pattern */ 895 PetscCall(MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu)); 896 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 897 /* Check for symmetric storage */ 898 PetscCall(MatGetType(A, &mtype)); 899 PetscCall(PetscStrcmp(mtype, MATSBAIJ, &isSymBlock)); 900 PetscCall(PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock)); 901 PetscCall(PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock)); 902 if (isSymBlock || isSymSeqBlock || isSymMPIBlock) PetscCall(MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE)); 903 /* Fill matrix with zeros */ 904 if (fillMatrix) { 905 if (Nf < 1 || bs > 1) { 906 PetscCall(DMGetBasicAdjacency(dm, &useCone, &useClosure)); 907 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 908 PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, -1, sectionAdj[idx], cols[idx], A)); 909 } else { 910 for (f = 0; f < Nf; ++f) { 911 PetscCall(DMGetAdjacency(dm, f, &useCone, &useClosure)); 912 idx = (useCone ? 1 : 0) + (useClosure ? 2 : 0); 913 PetscCall(DMPlexFillMatrix_Static(dm, rLayout, bs, f, sectionAdj[idx], cols[idx], A)); 914 } 915 } 916 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 917 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 918 } 919 PetscCall(PetscLayoutDestroy(&rLayout)); 920 for (idx = 0; idx < 4; ++idx) { 921 PetscCall(PetscSectionDestroy(§ionAdj[idx])); 922 PetscCall(PetscFree(cols[idx])); 923 } 924 PetscCall(PetscLogEventEnd(DMPLEX_Preallocate, dm, 0, 0, 0)); 925 PetscFunctionReturn(PETSC_SUCCESS); 926 } 927 928 #if 0 929 PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix) 930 { 931 PetscInt *tmpClosure,*tmpAdj,*visits; 932 PetscInt c,cStart,cEnd,pStart,pEnd; 933 934 PetscFunctionBegin; 935 PetscCall(DMGetDimension(dm, &dim)); 936 PetscCall(DMPlexGetDepth(dm, &depth)); 937 PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize)); 938 939 maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1)); 940 941 PetscCall(PetscSectionGetChart(section, &pStart, &pEnd)); 942 npoints = pEnd - pStart; 943 944 PetscCall(PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits)); 945 PetscCall(PetscArrayzero(lvisits,pEnd-pStart)); 946 PetscCall(PetscArrayzero(visits,pEnd-pStart)); 947 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 948 for (c=cStart; c<cEnd; c++) { 949 PetscInt *support = tmpClosure; 950 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support)); 951 for (p=0; p<supportSize; p++) lvisits[support[p]]++; 952 } 953 PetscCall(PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM)); 954 PetscCall(PetscSFReduceEnd (sf,MPIU_INT,lvisits,visits,MPI_SUM)); 955 PetscCall(PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits,MPI_REPLACE)); 956 PetscCall(PetscSFBcastEnd (sf,MPIU_INT,visits,lvisits)); 957 958 PetscCall(PetscSFGetRootRanks()); 959 960 PetscCall(PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner)); 961 for (c=cStart; c<cEnd; c++) { 962 PetscCall(PetscArrayzero(cellmat,maxClosureSize*maxClosureSize)); 963 /* 964 Depth-first walk of transitive closure. 965 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. 966 This contribution is added to dnz if owning ranks of p and q match, to onz otherwise. 967 */ 968 } 969 970 PetscCall(PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM)); 971 PetscCall(PetscSFReduceEnd (sf,MPIU_INT,lonz,onz,MPI_SUM)); 972 PetscFunctionReturn(PETSC_SUCCESS); 973 } 974 #endif 975