1. petsc
  2. PETSc
  3. petsc

Commits

Matt Knepley  committed 7f5b169

DMPlex: Proper preallocation for interpolators

  • Participants
  • Parent commits 75a6906
  • Branches master

Comments (0)

Files changed (2)

File src/dm/impls/plex/plex.c

View file
  • Ignore whitespace
   ierr = MatCreate(PetscObjectComm((PetscObject) dmCoarse), interpolation);CHKERRQ(ierr);
   ierr = MatSetSizes(*interpolation, m, n, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr);
   ierr = MatSetType(*interpolation, dmCoarse->mattype);CHKERRQ(ierr);
-  ierr = MatSetUp(*interpolation);CHKERRQ(ierr);
-  ierr = MatSetFromOptions(*interpolation);CHKERRQ(ierr);
-  ierr = MatSetOption(*interpolation, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHKERRQ(ierr);
   ierr = DMGetApplicationContext(dmFine, &ctx);CHKERRQ(ierr);
   ierr = DMPlexComputeInterpolatorFEM(dmCoarse, dmFine, *interpolation, ctx);CHKERRQ(ierr);
   /* Use naive scaling */

File src/dm/impls/plex/plexfem.c

View file
  • Ignore whitespace
     rTotDim += rNb*Nc;
   }
   ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr);
-  ierr = MatZeroEntries(In);CHKERRQ(ierr);
   ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr);
   ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr);
   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
     ierr = PetscFree(points);CHKERRQ(ierr);
   }
   if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);}
+  /* Preallocate matrix */
+  {
+    PetscHashJK ht;
+    PetscLayout rLayout;
+    PetscInt   *dnz, *onz, *cellCIndices, *cellFIndices;
+    PetscInt    locRows, rStart, rEnd, cell, r;
+
+    ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr);
+    ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr);
+    ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
+    ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
+    ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
+    ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
+    ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
+    ierr = PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr);
+    ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr);
+    for (cell = cStart; cell < cEnd; ++cell) {
+      ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr);
+      for (r = 0; r < rTotDim; ++r) {
+        PetscHashJKKey  key;
+        PetscHashJKIter missing, iter;
+
+        key.j = cellFIndices[r];
+        if (key.j < 0) continue;
+        for (c = 0; c < cTotDim; ++c) {
+          key.k = cellCIndices[c];
+          if (key.k < 0) continue;
+          ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr);
+          if (missing) {
+            ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr);
+            if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
+            else                                     ++onz[key.j-rStart];
+          }
+        }
+      }
+    }
+    ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr);
+    ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr);
+    ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
+    ierr = PetscFree4(dnz,onz,cellCIndices,cellFIndices);CHKERRQ(ierr);
+  }
+  /* Fill matrix */
+  ierr = MatZeroEntries(In);CHKERRQ(ierr);
   for (c = cStart; c < cEnd; ++c) {
     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr);
   }