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