Actual source code: dapreallocate.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  1:  #include <petsc/private/dmdaimpl.h>
  2:  #include <petsc/private/isimpl.h>
  3:  #include <petscsf.h>

  5: static PetscErrorCode DMDAGetAdjacency_Internal(DM dm, PetscInt p, PetscBool useClosure, const PetscInt *tmpClosure, PetscInt *adjSize, PetscInt adj[])
  6: {
  7:   const PetscInt *star  = tmpClosure;
  8:   PetscInt        numAdj = 0, maxAdjSize = *adjSize, starSize, s;
  9:   PetscErrorCode  ierr;

 12:   DMDAGetTransitiveClosure(dm, p, useClosure, &starSize, (PetscInt**) &star);
 13:   for (s = 2; s < starSize*2; s += 2) {
 14:     const PetscInt *closure = NULL;
 15:     PetscInt        closureSize, c, q;

 17:     DMDAGetTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);
 18:     for (c = 0; c < closureSize*2; c += 2) {
 19:       for (q = 0; q < numAdj || ((void)(adj[numAdj++] = closure[c]),0); ++q) {
 20:         if (closure[c] == adj[q]) break;
 21:       }
 22:       if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
 23:     }
 24:     DMDARestoreTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);
 25:   }
 26:   *adjSize = numAdj;
 27:   return(0);
 28: }

 30: /*@
 31:   DMDASetPreallocationCenterDimension - Determine the topology used to determine adjacency

 33:   Input Parameters:
 34: + dm - The DM object
 35: - preallocCenterDim - The dimension of points which connect adjacent entries

 37:   Level: developer

 39:   Notes:
 40: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim
 41: $     FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1
 42: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0

 44: .seealso: DMCreateMatrix(), DMDAPreallocateOperator()
 45: @*/
 46: PetscErrorCode DMDASetPreallocationCenterDimension(DM dm, PetscInt preallocCenterDim)
 47: {
 48:   DM_DA *mesh = (DM_DA*) dm->data;

 52:   mesh->preallocCenterDim = preallocCenterDim;
 53:   return(0);
 54: }

 56: /*@
 57:   DMDAGetPreallocationCenterDimension - Return the topology used to determine adjacency

 59:   Input Parameter:
 60: . dm - The DM object

 62:   Output Parameter:
 63: . preallocCenterDim - The dimension of points which connect adjacent entries

 65:   Level: developer

 67:   Notes:
 68: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim
 69: $     FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1
 70: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0

 72: .seealso: DMCreateMatrix(), DMDAPreallocateOperator(), DMDASetPreallocationCenterDimension()
 73: @*/
 74: PetscErrorCode DMDAGetPreallocationCenterDimension(DM dm, PetscInt *preallocCenterDim)
 75: {
 76:   DM_DA *mesh = (DM_DA*) dm->data;

 81:   *preallocCenterDim = mesh->preallocCenterDim;
 82:   return(0);
 83: }

 85: PetscErrorCode DMDAPreallocateOperator(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
 86: {
 87:   MPI_Comm           comm;
 88:   MatType            mtype;
 89:   PetscSF            sf, sfDof, sfAdj;
 90:   PetscSection       leafSectionAdj, rootSectionAdj, sectionAdj;
 91:   PetscInt           nleaves, l, p;
 92:   const PetscInt    *leaves;
 93:   const PetscSFNode *remotes;
 94:   PetscInt           dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
 95:   PetscInt          *tmpClosure, *tmpAdj, *adj, *rootAdj, *cols, *remoteOffsets;
 96:   PetscInt           depth, centerDim, maxConeSize, maxSupportSize, maxClosureSize, maxAdjSize, adjSize;
 97:   PetscLayout        rLayout;
 98:   PetscInt           locRows, rStart, rEnd, r;
 99:   PetscMPIInt        size;
100:   PetscBool          useClosure, debug = PETSC_FALSE, isSymBlock, isSymSeqBlock, isSymMPIBlock;
101:   PetscErrorCode     ierr;

104:   PetscObjectGetComm((PetscObject)dm,&comm);
105:   PetscOptionsGetBool(NULL,NULL, "-dm_view_preallocation", &debug, NULL);
106:   MPI_Comm_size(comm, &size);
107:   DMDAGetInfo(dm, &dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
108:   depth = dim;
109:   DMGetPointSF(dm, &sf);
110:   /* Create dof SF based on point SF */
111:   if (debug) {
112:     PetscPrintf(comm, "Input Section for Preallocation:\n");
113:     PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);
114:     PetscPrintf(comm, "Input Global Section for Preallocation:\n");
115:     PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);
116:     PetscPrintf(comm, "Input SF for Preallocation:\n");
117:     PetscSFView(sf, NULL);
118:   }
119:   PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);
120:   PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);
121:   PetscFree(remoteOffsets);
122:   if (debug) {
123:     PetscPrintf(comm, "Dof SF for Preallocation:\n");
124:     PetscSFView(sfDof, NULL);
125:   }
126:   /* Create section for dof adjacency (dof ==> # adj dof) */
127:   /*   FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim   */
128:   /*   FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1 */
129:   /*   FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0     */
130:   DMDAGetPreallocationCenterDimension(dm, &centerDim);
131:   if (centerDim == dim) {
132:     useClosure = PETSC_FALSE;
133:   } else if (centerDim == 0) {
134:     useClosure = PETSC_TRUE;
135:   } else SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Do not support preallocation with center points of dimension %d", centerDim);

137:   PetscSectionGetChart(section, &pStart, &pEnd);
138:   PetscSectionGetStorageSize(section, &numDof);
139:   PetscSectionCreate(comm, &leafSectionAdj);
140:   PetscSectionSetChart(leafSectionAdj, 0, numDof);
141:   PetscSectionCreate(comm, &rootSectionAdj);
142:   PetscSectionSetChart(rootSectionAdj, 0, numDof);
143:   /*   Fill in the ghost dofs on the interface */
144:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);
145:   maxConeSize    = 6;
146:   maxSupportSize = 6;

148:   maxClosureSize = 2*PetscMax(PetscPowInt(maxConeSize,depth+1),PetscPowInt(maxSupportSize,depth+1));
149:   maxAdjSize     = PetscPowInt(maxConeSize,depth+1) * PetscPowInt(maxSupportSize,depth+1) + 1;

151:   PetscMalloc2(maxClosureSize,&tmpClosure,maxAdjSize,&tmpAdj);

153:   /*
154:    ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
155:     1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
156:        Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
157:     2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
158:        Create sfAdj connecting rootSectionAdj and leafSectionAdj
159:     3. Visit unowned points on interface, write adjacencies to adj
160:        Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
161:     4. Visit owned points on interface, write adjacencies to rootAdj
162:        Remove redundancy in rootAdj
163:    ** The last two traversals use transitive closure
164:     5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
165:        Allocate memory addressed by sectionAdj (cols)
166:     6. Visit all owned points in the subdomain, insert dof adjacencies into cols
167:    ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
168:   */

170:   for (l = 0; l < nleaves; ++l) {
171:     PetscInt dof, off, d, q;
172:     PetscInt p = leaves[l], numAdj = maxAdjSize;

174:     if ((p < pStart) || (p >= pEnd)) continue;
175:     PetscSectionGetDof(section, p, &dof);
176:     PetscSectionGetOffset(section, p, &off);
177:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
178:     for (q = 0; q < numAdj; ++q) {
179:       const PetscInt padj = tmpAdj[q];
180:       PetscInt ndof, ncdof;

182:       if ((padj < pStart) || (padj >= pEnd)) continue;
183:       PetscSectionGetDof(section, padj, &ndof);
184:       PetscSectionGetConstraintDof(section, padj, &ncdof);
185:       for (d = off; d < off+dof; ++d) {
186:         PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);
187:       }
188:     }
189:   }
190:   PetscSectionSetUp(leafSectionAdj);
191:   if (debug) {
192:     PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");
193:     PetscSectionView(leafSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
194:   }
195:   /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
196:   if (size > 1) {
197:     PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
198:     PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
199:   }
200:   if (debug) {
201:     PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");
202:     PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
203:   }
204:   /* Add in local adjacency sizes for owned dofs on interface (roots) */
205:   for (p = pStart; p < pEnd; ++p) {
206:     PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;

208:     PetscSectionGetDof(section, p, &dof);
209:     PetscSectionGetOffset(section, p, &off);
210:     if (!dof) continue;
211:     PetscSectionGetDof(rootSectionAdj, off, &adof);
212:     if (adof <= 0) continue;
213:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
214:     for (q = 0; q < numAdj; ++q) {
215:       const PetscInt padj = tmpAdj[q];
216:       PetscInt ndof, ncdof;

218:       if ((padj < pStart) || (padj >= pEnd)) continue;
219:       PetscSectionGetDof(section, padj, &ndof);
220:       PetscSectionGetConstraintDof(section, padj, &ncdof);
221:       for (d = off; d < off+dof; ++d) {
222:         PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);
223:       }
224:     }
225:   }
226:   PetscSectionSetUp(rootSectionAdj);
227:   if (debug) {
228:     PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");
229:     PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
230:   }
231:   /* Create adj SF based on dof SF */
232:   PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);
233:   PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);
234:   PetscFree(remoteOffsets);
235:   if (debug) {
236:     PetscPrintf(comm, "Adjacency SF for Preallocation:\n");
237:     PetscSFView(sfAdj, NULL);
238:   }
239:   PetscSFDestroy(&sfDof);
240:   /* Create leaf adjacency */
241:   PetscSectionSetUp(leafSectionAdj);
242:   PetscSectionGetStorageSize(leafSectionAdj, &adjSize);
243:   PetscMalloc1(adjSize, &adj);
244:   PetscMemzero(adj, adjSize * sizeof(PetscInt));
245:   for (l = 0; l < nleaves; ++l) {
246:     PetscInt dof, off, d, q;
247:     PetscInt p = leaves[l], numAdj = maxAdjSize;

249:     if ((p < pStart) || (p >= pEnd)) continue;
250:     PetscSectionGetDof(section, p, &dof);
251:     PetscSectionGetOffset(section, p, &off);
252:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
253:     for (d = off; d < off+dof; ++d) {
254:       PetscInt aoff, i = 0;

256:       PetscSectionGetOffset(leafSectionAdj, d, &aoff);
257:       for (q = 0; q < numAdj; ++q) {
258:         const PetscInt padj = tmpAdj[q];
259:         PetscInt ndof, ncdof, ngoff, nd;

261:         if ((padj < pStart) || (padj >= pEnd)) continue;
262:         PetscSectionGetDof(section, padj, &ndof);
263:         PetscSectionGetConstraintDof(section, padj, &ncdof);
264:         PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
265:         for (nd = 0; nd < ndof-ncdof; ++nd) {
266:           adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
267:           ++i;
268:         }
269:       }
270:     }
271:   }
272:   /* Debugging */
273:   if (debug) {
274:     IS tmp;
275:     PetscPrintf(comm, "Leaf adjacency indices\n");
276:     ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);
277:     ISView(tmp, NULL);
278:     ISDestroy(&tmp);
279:   }
280:   /* Gather adjacenct indices to root */
281:   PetscSectionGetStorageSize(rootSectionAdj, &adjSize);
282:   PetscMalloc1(adjSize, &rootAdj);
283:   for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
284:   if (size > 1) {
285:     PetscSFGatherBegin(sfAdj, MPIU_INT, adj, rootAdj);
286:     PetscSFGatherEnd(sfAdj, MPIU_INT, adj, rootAdj);
287:   }
288:   PetscSFDestroy(&sfAdj);
289:   PetscFree(adj);
290:   /* Debugging */
291:   if (debug) {
292:     IS tmp;
293:     PetscPrintf(comm, "Root adjacency indices after gather\n");
294:     ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
295:     ISView(tmp, NULL);
296:     ISDestroy(&tmp);
297:   }
298:   /* Add in local adjacency indices for owned dofs on interface (roots) */
299:   for (p = pStart; p < pEnd; ++p) {
300:     PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;

302:     PetscSectionGetDof(section, p, &dof);
303:     PetscSectionGetOffset(section, p, &off);
304:     if (!dof) continue;
305:     PetscSectionGetDof(rootSectionAdj, off, &adof);
306:     if (adof <= 0) continue;
307:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
308:     for (d = off; d < off+dof; ++d) {
309:       PetscInt adof, aoff, i;

311:       PetscSectionGetDof(rootSectionAdj, d, &adof);
312:       PetscSectionGetOffset(rootSectionAdj, d, &aoff);
313:       i    = adof-1;
314:       for (q = 0; q < numAdj; ++q) {
315:         const PetscInt padj = tmpAdj[q];
316:         PetscInt ndof, ncdof, ngoff, nd;

318:         if ((padj < pStart) || (padj >= pEnd)) continue;
319:         PetscSectionGetDof(section, padj, &ndof);
320:         PetscSectionGetConstraintDof(section, padj, &ncdof);
321:         PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
322:         for (nd = 0; nd < ndof-ncdof; ++nd) {
323:           rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
324:           --i;
325:         }
326:       }
327:     }
328:   }
329:   /* Debugging */
330:   if (debug) {
331:     IS tmp;
332:     PetscPrintf(comm, "Root adjacency indices\n");
333:     ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
334:     ISView(tmp, NULL);
335:     ISDestroy(&tmp);
336:   }
337:   /* Compress indices */
338:   PetscSectionSetUp(rootSectionAdj);
339:   for (p = pStart; p < pEnd; ++p) {
340:     PetscInt dof, cdof, off, d;
341:     PetscInt adof, aoff;

343:     PetscSectionGetDof(section, p, &dof);
344:     PetscSectionGetConstraintDof(section, p, &cdof);
345:     PetscSectionGetOffset(section, p, &off);
346:     if (!dof) continue;
347:     PetscSectionGetDof(rootSectionAdj, off, &adof);
348:     if (adof <= 0) continue;
349:     for (d = off; d < off+dof-cdof; ++d) {
350:       PetscSectionGetDof(rootSectionAdj, d, &adof);
351:       PetscSectionGetOffset(rootSectionAdj, d, &aoff);
352:       PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);
353:       PetscSectionSetDof(rootSectionAdj, d, adof);
354:     }
355:   }
356:   /* Debugging */
357:   if (debug) {
358:     IS tmp;
359:     PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");
360:     PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
361:     PetscPrintf(comm, "Root adjacency indices after compression\n");
362:     ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
363:     ISView(tmp, NULL);
364:     ISDestroy(&tmp);
365:   }
366:   /* Build adjacency section: Maps global indices to sets of adjacent global indices */
367:   PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);
368:   PetscSectionCreate(comm, &sectionAdj);
369:   PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);
370:   for (p = pStart; p < pEnd; ++p) {
371:     PetscInt  numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
372:     PetscBool found  = PETSC_TRUE;

374:     PetscSectionGetDof(section, p, &dof);
375:     PetscSectionGetConstraintDof(section, p, &cdof);
376:     PetscSectionGetOffset(section, p, &off);
377:     PetscSectionGetOffset(sectionGlobal, p, &goff);
378:     for (d = 0; d < dof-cdof; ++d) {
379:       PetscInt ldof, rdof;

381:       PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
382:       PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
383:       if (ldof > 0) {
384:         /* We do not own this point */
385:       } else if (rdof > 0) {
386:         PetscSectionSetDof(sectionAdj, goff+d, rdof);
387:       } else {
388:         found = PETSC_FALSE;
389:       }
390:     }
391:     if (found) continue;
392:     PetscSectionGetDof(section, p, &dof);
393:     PetscSectionGetOffset(sectionGlobal, p, &goff);
394:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
395:     for (q = 0; q < numAdj; ++q) {
396:       const PetscInt padj = tmpAdj[q];
397:       PetscInt ndof, ncdof, noff;

399:       if ((padj < pStart) || (padj >= pEnd)) continue;
400:       PetscSectionGetDof(section, padj, &ndof);
401:       PetscSectionGetConstraintDof(section, padj, &ncdof);
402:       PetscSectionGetOffset(section, padj, &noff);
403:       for (d = goff; d < goff+dof-cdof; ++d) {
404:         PetscSectionAddDof(sectionAdj, d, ndof-ncdof);
405:       }
406:     }
407:   }
408:   PetscSectionSetUp(sectionAdj);
409:   if (debug) {
410:     PetscPrintf(comm, "Adjacency Section for Preallocation:\n");
411:     PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);
412:   }
413:   /* Get adjacent indices */
414:   PetscSectionGetStorageSize(sectionAdj, &numCols);
415:   PetscMalloc1(numCols, &cols);
416:   for (p = pStart; p < pEnd; ++p) {
417:     PetscInt  numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
418:     PetscBool found  = PETSC_TRUE;

420:     PetscSectionGetDof(section, p, &dof);
421:     PetscSectionGetConstraintDof(section, p, &cdof);
422:     PetscSectionGetOffset(section, p, &off);
423:     PetscSectionGetOffset(sectionGlobal, p, &goff);
424:     for (d = 0; d < dof-cdof; ++d) {
425:       PetscInt ldof, rdof;

427:       PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
428:       PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
429:       if (ldof > 0) {
430:         /* We do not own this point */
431:       } else if (rdof > 0) {
432:         PetscInt aoff, roff;

434:         PetscSectionGetOffset(sectionAdj, goff+d, &aoff);
435:         PetscSectionGetOffset(rootSectionAdj, off+d, &roff);
436:         PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));
437:       } else {
438:         found = PETSC_FALSE;
439:       }
440:     }
441:     if (found) continue;
442:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
443:     for (d = goff; d < goff+dof-cdof; ++d) {
444:       PetscInt adof, aoff, i = 0;

446:       PetscSectionGetDof(sectionAdj, d, &adof);
447:       PetscSectionGetOffset(sectionAdj, d, &aoff);
448:       for (q = 0; q < numAdj; ++q) {
449:         const PetscInt  padj = tmpAdj[q];
450:         PetscInt        ndof, ncdof, ngoff, nd;
451:         const PetscInt *ncind;

453:         /* Adjacent points may not be in the section chart */
454:         if ((padj < pStart) || (padj >= pEnd)) continue;
455:         PetscSectionGetDof(section, padj, &ndof);
456:         PetscSectionGetConstraintDof(section, padj, &ncdof);
457:         PetscSectionGetConstraintIndices(section, padj, &ncind);
458:         PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
459:         for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
460:           cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
461:         }
462:       }
463:       if (i != adof) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of entries %D != %D for dof %D (point %D)", i, adof, d, p);
464:     }
465:   }
466:   PetscSectionDestroy(&leafSectionAdj);
467:   PetscSectionDestroy(&rootSectionAdj);
468:   PetscFree(rootAdj);
469:   PetscFree2(tmpClosure, tmpAdj);
470:   /* Debugging */
471:   if (debug) {
472:     IS tmp;
473:     PetscPrintf(comm, "Column indices\n");
474:     ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);
475:     ISView(tmp, NULL);
476:     ISDestroy(&tmp);
477:   }
478:   /* Create allocation vectors from adjacency graph */
479:   MatGetLocalSize(A, &locRows, NULL);
480:   PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);
481:   PetscLayoutSetLocalSize(rLayout, locRows);
482:   PetscLayoutSetBlockSize(rLayout, 1);
483:   PetscLayoutSetUp(rLayout);
484:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
485:   PetscLayoutDestroy(&rLayout);
486:   /* Only loop over blocks of rows */
487:   if (rStart%bs || rEnd%bs) SETERRQ3(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid layout [%d, %d) for matrix, must be divisible by block size %d", rStart, rEnd, bs);
488:   for (r = rStart/bs; r < rEnd/bs; ++r) {
489:     const PetscInt row = r*bs;
490:     PetscInt       numCols, cStart, c;

492:     PetscSectionGetDof(sectionAdj, row, &numCols);
493:     PetscSectionGetOffset(sectionAdj, row, &cStart);
494:     for (c = cStart; c < cStart+numCols; ++c) {
495:       if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
496:         ++dnz[r-rStart];
497:         if (cols[c] >= row) ++dnzu[r-rStart];
498:       } else {
499:         ++onz[r-rStart];
500:         if (cols[c] >= row) ++onzu[r-rStart];
501:       }
502:     }
503:   }
504:   if (bs > 1) {
505:     for (r = 0; r < locRows/bs; ++r) {
506:       dnz[r]  /= bs;
507:       onz[r]  /= bs;
508:       dnzu[r] /= bs;
509:       onzu[r] /= bs;
510:     }
511:   }
512:   /* Set matrix pattern */
513:   MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);
514:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
515:   /* Check for symmetric storage */
516:   MatGetType(A, &mtype);
517:   PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
518:   PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
519:   PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
520:   if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);}
521:   /* Fill matrix with zeros */
522:   if (fillMatrix) {
523:     PetscScalar *values;
524:     PetscInt     maxRowLen = 0;

526:     for (r = rStart; r < rEnd; ++r) {
527:       PetscInt len;

529:       PetscSectionGetDof(sectionAdj, r, &len);
530:       maxRowLen = PetscMax(maxRowLen, len);
531:     }
532:     PetscCalloc1(maxRowLen, &values);
533:     for (r = rStart; r < rEnd; ++r) {
534:       PetscInt numCols, cStart;

536:       PetscSectionGetDof(sectionAdj, r, &numCols);
537:       PetscSectionGetOffset(sectionAdj, r, &cStart);
538:       MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
539:     }
540:     PetscFree(values);
541:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
542:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
543:   }
544:   PetscSectionDestroy(&sectionAdj);
545:   PetscFree(cols);
546:   return(0);
547: }