Actual source code: dapreallocate.c

petsc-3.6.1 2015-08-06
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, "-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:   if (debug) {
130:     PetscPrintf(comm, "Dof SF for Preallocation:\n");
131:     PetscSFView(sfDof, NULL);
132:   }
133:   /* Create section for dof adjacency (dof ==> # adj dof) */
134:   /*   FEM:   Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim   */
135:   /*   FVM:   Two points p and q are adjacent if q \in star(cone(p)),    preallocCenterDim = dim-1 */
136:   /*   FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0     */
137:   DMDAGetPreallocationCenterDimension(dm, &centerDim);
138:   if (centerDim == dim) {
139:     useClosure = PETSC_FALSE;
140:   } else if (centerDim == 0) {
141:     useClosure = PETSC_TRUE;
142:   } else SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Do not support preallocation with center points of dimension %d", centerDim);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

532:     for (r = rStart; r < rEnd; ++r) {
533:       PetscInt len;

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

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