Actual source code: dapreallocate.c

petsc-3.7.4 2016-10-02
Report Typos and Errors
  1: #include <petsc/private/dmdaimpl.h>   /*I      "petscdmda.h"   I*/
  2: #include <petsc/private/isimpl.h>
  3: #include <petscsf.h>

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

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

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

 34: /*@
 35:   DMDASetPreallocationCenterDimension - Determine the topology used to determine adjacency

 37:   Input Parameters:
 38: + dm - The DM object
 39: - preallocCenterDim - The dimension of points which connect adjacent entries

 41:   Level: developer

 43:   Notes:
 44: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim
 45: $     FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1
 46: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0

 48: .seealso: DMCreateMatrix(), DMDAPreallocateOperator()
 49: @*/
 50: PetscErrorCode DMDASetPreallocationCenterDimension(DM dm, PetscInt preallocCenterDim)
 51: {
 52:   DM_DA *mesh = (DM_DA*) dm->data;

 56:   mesh->preallocCenterDim = preallocCenterDim;
 57:   return(0);
 58: }

 62: /*@
 63:   DMDAGetPreallocationCenterDimension - Return the topology used to determine adjacency

 65:   Input Parameter:
 66: . dm - The DM object

 68:   Output Parameter:
 69: . preallocCenterDim - The dimension of points which connect adjacent entries

 71:   Level: developer

 73:   Notes:
 74: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim
 75: $     FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1
 76: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0

 78: .seealso: DMCreateMatrix(), DMDAPreallocateOperator(), DMDASetPreallocationCenterDimension()
 79: @*/
 80: PetscErrorCode DMDAGetPreallocationCenterDimension(DM dm, PetscInt *preallocCenterDim)
 81: {
 82:   DM_DA *mesh = (DM_DA*) dm->data;

 87:   *preallocCenterDim = mesh->preallocCenterDim;
 88:   return(0);
 89: }

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

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

145:   PetscSectionGetChart(section, &pStart, &pEnd);
146:   PetscSectionGetStorageSize(section, &numDof);
147:   PetscSectionCreate(comm, &leafSectionAdj);
148:   PetscSectionSetChart(leafSectionAdj, 0, numDof);
149:   PetscSectionCreate(comm, &rootSectionAdj);
150:   PetscSectionSetChart(rootSectionAdj, 0, numDof);
151:   /*   Fill in the ghost dofs on the interface */
152:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);
153:   maxConeSize    = 6;
154:   maxSupportSize = 6;

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

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

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

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

182:     if ((p < pStart) || (p >= pEnd)) continue;
183:     PetscSectionGetDof(section, p, &dof);
184:     PetscSectionGetOffset(section, p, &off);
185:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
186:     for (q = 0; q < numAdj; ++q) {
187:       const PetscInt padj = tmpAdj[q];
188:       PetscInt ndof, ncdof;

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

216:     PetscSectionGetDof(section, p, &dof);
217:     PetscSectionGetOffset(section, p, &off);
218:     if (!dof) continue;
219:     PetscSectionGetDof(rootSectionAdj, off, &adof);
220:     if (adof <= 0) continue;
221:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
222:     for (q = 0; q < numAdj; ++q) {
223:       const PetscInt padj = tmpAdj[q];
224:       PetscInt ndof, ncdof;

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

257:     if ((p < pStart) || (p >= pEnd)) continue;
258:     PetscSectionGetDof(section, p, &dof);
259:     PetscSectionGetOffset(section, p, &off);
260:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
261:     for (d = off; d < off+dof; ++d) {
262:       PetscInt aoff, i = 0;

264:       PetscSectionGetOffset(leafSectionAdj, d, &aoff);
265:       for (q = 0; q < numAdj; ++q) {
266:         const PetscInt padj = tmpAdj[q];
267:         PetscInt ndof, ncdof, ngoff, nd;

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

310:     PetscSectionGetDof(section, p, &dof);
311:     PetscSectionGetOffset(section, p, &off);
312:     if (!dof) continue;
313:     PetscSectionGetDof(rootSectionAdj, off, &adof);
314:     if (adof <= 0) continue;
315:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
316:     for (d = off; d < off+dof; ++d) {
317:       PetscInt adof, aoff, i;

319:       PetscSectionGetDof(rootSectionAdj, d, &adof);
320:       PetscSectionGetOffset(rootSectionAdj, d, &aoff);
321:       i    = adof-1;
322:       for (q = 0; q < numAdj; ++q) {
323:         const PetscInt padj = tmpAdj[q];
324:         PetscInt ndof, ncdof, ngoff, nd;

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

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

382:     PetscSectionGetDof(section, p, &dof);
383:     PetscSectionGetConstraintDof(section, p, &cdof);
384:     PetscSectionGetOffset(section, p, &off);
385:     PetscSectionGetOffset(sectionGlobal, p, &goff);
386:     for (d = 0; d < dof-cdof; ++d) {
387:       PetscInt ldof, rdof;

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

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

428:     PetscSectionGetDof(section, p, &dof);
429:     PetscSectionGetConstraintDof(section, p, &cdof);
430:     PetscSectionGetOffset(section, p, &off);
431:     PetscSectionGetOffset(sectionGlobal, p, &goff);
432:     for (d = 0; d < dof-cdof; ++d) {
433:       PetscInt ldof, rdof;

435:       PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
436:       PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
437:       if (ldof > 0) {
438:         /* We do not own this point */
439:       } else if (rdof > 0) {
440:         PetscInt aoff, roff;

442:         PetscSectionGetOffset(sectionAdj, goff+d, &aoff);
443:         PetscSectionGetOffset(rootSectionAdj, off+d, &roff);
444:         PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));
445:       } else {
446:         found = PETSC_FALSE;
447:       }
448:     }
449:     if (found) continue;
450:     DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
451:     for (d = goff; d < goff+dof-cdof; ++d) {
452:       PetscInt adof, aoff, i = 0;

454:       PetscSectionGetDof(sectionAdj, d, &adof);
455:       PetscSectionGetOffset(sectionAdj, d, &aoff);
456:       for (q = 0; q < numAdj; ++q) {
457:         const PetscInt  padj = tmpAdj[q];
458:         PetscInt        ndof, ncdof, ngoff, nd;
459:         const PetscInt *ncind;

461:         /* Adjacent points may not be in the section chart */
462:         if ((padj < pStart) || (padj >= pEnd)) continue;
463:         PetscSectionGetDof(section, padj, &ndof);
464:         PetscSectionGetConstraintDof(section, padj, &ncdof);
465:         PetscSectionGetConstraintIndices(section, padj, &ncind);
466:         PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
467:         for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
468:           cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
469:         }
470:       }
471:       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);
472:     }
473:   }
474:   PetscSectionDestroy(&leafSectionAdj);
475:   PetscSectionDestroy(&rootSectionAdj);
476:   PetscFree(rootAdj);
477:   PetscFree2(tmpClosure, tmpAdj);
478:   /* Debugging */
479:   if (debug) {
480:     IS tmp;
481:     PetscPrintf(comm, "Column indices\n");
482:     ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);
483:     ISView(tmp, NULL);
484:     ISDestroy(&tmp);
485:   }
486:   /* Create allocation vectors from adjacency graph */
487:   MatGetLocalSize(A, &locRows, NULL);
488:   PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);
489:   PetscLayoutSetLocalSize(rLayout, locRows);
490:   PetscLayoutSetBlockSize(rLayout, 1);
491:   PetscLayoutSetUp(rLayout);
492:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
493:   PetscLayoutDestroy(&rLayout);
494:   /* Only loop over blocks of rows */
495:   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);
496:   for (r = rStart/bs; r < rEnd/bs; ++r) {
497:     const PetscInt row = r*bs;
498:     PetscInt       numCols, cStart, c;

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

534:     for (r = rStart; r < rEnd; ++r) {
535:       PetscInt len;

537:       PetscSectionGetDof(sectionAdj, r, &len);
538:       maxRowLen = PetscMax(maxRowLen, len);
539:     }
540:     PetscCalloc1(maxRowLen, &values);
541:     for (r = rStart; r < rEnd; ++r) {
542:       PetscInt numCols, cStart;

544:       PetscSectionGetDof(sectionAdj, r, &numCols);
545:       PetscSectionGetOffset(sectionAdj, r, &cStart);
546:       MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
547:     }
548:     PetscFree(values);
549:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
550:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
551:   }
552:   PetscSectionDestroy(&sectionAdj);
553:   PetscFree(cols);
554:   return(0);
555: }