Dear petsc devs
After some profiling of the loading of a large gmsh mesh in dmplex
(which seemed to take forever), I found that it was spending 99% of its
time in IsGetIndices_Stride in the following calltree
DMPlexCreateGmshFromFile -> DMPlexGetFullJoin -> DMLabelGetStratumBounds
-> ISGetIndices -> ISGetIndices_Stride
where it seems to work out the bounds of each stratum by explicitly
computing all indices of a strided index set to then only take its first
and last entry, and this is done for each facet in the mesh a couple of
times.
With the attached patch I managed to reduce the load time of a mesh from
4 hours to 22 seconds! I'm not a 100% sure if the patch is correct in
all circumstances: I presume the assumption was that the index set is
sorted?
Cheers
Stephan
diff --git a/src/dm/label/dmlabel.c b/src/dm/label/dmlabel.c
index f3e784c..fb8d1b9 100644
--- a/src/dm/label/dmlabel.c
+++ b/src/dm/label/dmlabel.c
@@ -770,15 +770,14 @@ PetscErrorCode DMLabelGetStratumBounds(DMLabel label, PetscInt value, PetscInt *
if (start) {PetscValidPointer(start, 3); *start = 0;}
if (end) {PetscValidPointer(end, 4); *end = 0;}
for (v = 0; v < label->numStrata; ++v) {
- const PetscInt *points;
+ PetscInt min, max;
if (label->stratumValues[v] != value) continue;
ierr = DMLabelMakeValid_Private(label, v);CHKERRQ(ierr);
if (label->stratumSizes[v] <= 0) break;
- ierr = ISGetIndices(label->points[v],&points);CHKERRQ(ierr);
- if (start) *start = points[0];
- if (end) *end = points[label->stratumSizes[v]-1]+1;
- ierr = ISRestoreIndices(label->points[v],&points);CHKERRQ(ierr);
+ ierr = ISGetMinMax(label->points[v], &min, &max);CHKERRQ(ierr);
+ if (start) *start = min;
+ if (end) *end = max+1;
break;
}
PetscFunctionReturn(0);
