Matt Knepley avatar Matt Knepley committed 1fa8f14 Merge with conflicts

Merge branch 'knepley/feature-plex-reordering'

* knepley/feature-plex-reordering:
doc: update makefile with order.h -> ../../../include/petsc-private/matorderimpl.h rename
DMPlex: Fixed bug with DMPlexPermute - Fixed label permutation
DMPlex: Check whether default section exists to be reordered
DMPlex: Fixed compiler warning
DMPlex ex10: Fixed test output
DMPlex: Added missing declarations
DMPlex: Fixed misuse of PetscObjectSetName()
MatComputeBandwidth: use PetscInt instead of PetscMPIInt
Mat: Fixed declaration
Mat: Changed MatCalcBandwidth() to MatComputeBandwidth() - Added Fortran interface
DMPlex ex10: This tests reordering of meshes
DMPlex: Added DMPlexGetOrdering() and DMPlexPermute()
PetscSection: Added PetscSectionPermute()
DMLabel: Added DMLabelPermute()
DMPlex: Added DMPlexCreateDoublet()
DMPlex: Added DMPlexCopyLabels()
Mat: Added MatCalcBandwidth()
MatOrdering: Moved private header into petsc-private - Fix to pointer checks in sorder.c

Conflicts:
include/petscdmplex.h

Comments (0)

Files changed (43)

config/builder.py

                                                                  # 3D Hex P_2 scalar tests
                                                                  # 3D Hex P_2 vector tests
                                                                  ],
+                        'src/dm/impls/plex/examples/tests/ex10': [# Two cell tests
+                                                                  {'numProcs': 1, 'args': '-dim 2 -interpolate 1 -cell_simplex 1 -num_dof 1,0,0 -mat_view'},
+                                                                  {'numProcs': 1, 'args': '-dim 2 -interpolate 1 -cell_simplex 0 -num_dof 1,0,0 -mat_view'},
+                                                                  {'numProcs': 1, 'args': '-dim 3 -interpolate 1 -cell_simplex 1 -num_dof 1,0,0,0 -mat_view'},
+                                                                  {'numProcs': 1, 'args': '-dim 3 -interpolate 1 -cell_simplex 0 -num_dof 1,0,0,0 -mat_view'},
+                                                                  # Refined tests
+                                                                  {'numProcs': 1, 'args': '-dim 2 -interpolate 1 -cell_simplex 1 -refinement_limit 0.00625 -num_dof 1,0,0'},
+                                                                  {'numProcs': 1, 'args': '-dim 2 -interpolate 1 -cell_simplex 0 -refinement_uniform       -num_dof 1,0,0'},
+                                                                  {'numProcs': 1, 'args': '-dim 3 -interpolate 1 -cell_simplex 1 -refinement_limit 0.00625 -num_dof 1,0,0,0'},
+                                                                  {'numProcs': 1, 'args': '-dim 3 -interpolate 1 -cell_simplex 0 -refinement_uniform       -num_dof 1,0,0,0'},
+                                                                  # Parallel tests
+                                                                  ],
                         'src/dm/impls/plex/examples/tests/ex1f90': [{'numProcs': 1, 'args': ''}],
                         'src/dm/impls/plex/examples/tests/ex2f90': [{'numProcs': 1, 'args': ''}],
                         'src/dm/impls/plex/examples/tutorials/ex1': [{'numProcs': 1, 'args': ''},

include/petsc-private/matorderimpl.h

+#ifndef __MATORDERIMPL_H
+#define __MATORDERIMPL_H
+
+#include <petscmat.h>
+#include <petsc-private/petscimpl.h>
+
+/*
+   Defines the interface to the SparsePack routines, translated into C.
+*/
+PETSC_EXTERN PetscErrorCode SPARSEPACKgen1wd(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKgennd(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKgenrcm(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKgenqmd(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+
+PETSC_EXTERN PetscErrorCode SPARSEPACKqmdrch(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKqmdmrg(const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKqmdqt(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKqmdupd(const PetscInt*,const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKfnroot(PetscInt*,const PetscInt*,const PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKrootls(const PetscInt*,const PetscInt*,const PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKfn1wd(PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKrevrse(const PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKrootls(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKfndsep(PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKdegree(const PetscInt*,const PetscInt*,const PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscInt*);
+PETSC_EXTERN PetscErrorCode SPARSEPACKrcm(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
+
+#endif

include/petscdmplex.h

 PETSC_EXTERN PetscErrorCode DMLabelDestroyIndex(DMLabel);
 PETSC_EXTERN PetscErrorCode DMLabelHasPoint(DMLabel, PetscInt, PetscBool *);
 PETSC_EXTERN PetscErrorCode DMLabelFilter(DMLabel, PetscInt, PetscInt);
+PETSC_EXTERN PetscErrorCode DMLabelPermute(DMLabel, IS, DMLabel *);
 
 PETSC_EXTERN PetscErrorCode DMPlexCreateLabel(DM, const char []);
 PETSC_EXTERN PetscErrorCode DMPlexGetLabelValue(DM, const char[], PetscInt, PetscInt *);
 PETSC_EXTERN PetscErrorCode DMPlexInvertCell(PetscInt, PetscInt, int []);
 PETSC_EXTERN PetscErrorCode DMPlexInterpolate(DM, DM *);
 PETSC_EXTERN PetscErrorCode DMPlexCopyCoordinates(DM, DM);
+PETSC_EXTERN PetscErrorCode DMPlexCopyLabels(DM, DM);
 PETSC_EXTERN PetscErrorCode DMPlexDistribute(DM, const char[], PetscInt, PetscSF*, DM*);
 PETSC_EXTERN PetscErrorCode DMPlexDistributeField(DM,PetscSF,PetscSection,Vec,PetscSection,Vec);
 PETSC_EXTERN PetscErrorCode DMPlexDistributeData(DM,PetscSF,PetscSection,MPI_Datatype,void*,PetscSection,void**);
+PETSC_EXTERN PetscErrorCode DMPlexGetOrdering(DM, MatOrderingType, IS *);
+PETSC_EXTERN PetscErrorCode DMPlexPermute(DM, IS, DM *);
 PETSC_EXTERN PetscErrorCode DMPlexLoad(PetscViewer, DM);
 PETSC_EXTERN PetscErrorCode DMPlexGetSubpointMap(DM, DMLabel*);
 PETSC_EXTERN PetscErrorCode DMPlexSetSubpointMap(DM, DMLabel);
 PETSC_EXTERN PetscErrorCode DMPlexCreateSubpointIS(DM, IS *);
 
+PETSC_EXTERN PetscErrorCode DMPlexCreateDoublet(MPI_Comm, PetscInt, PetscBool, PetscBool, PetscBool, PetscReal, DM *);
 PETSC_EXTERN PetscErrorCode DMPlexCreateCubeBoundary(DM, const PetscReal [], const PetscReal [], const PetscInt []);
 PETSC_EXTERN PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm, PetscInt, PetscBool, DM *);
 PETSC_EXTERN PetscErrorCode DMPlexCreateHexBoxMesh(MPI_Comm,PetscInt,const PetscInt[],DM *);

include/petscis.h

 PETSC_EXTERN PetscErrorCode PetscSectionCreateSubmeshSection(PetscSection, IS, PetscSection *);
 PETSC_EXTERN PetscErrorCode PetscSectionGetPointLayout(MPI_Comm, PetscSection, PetscLayout *);
 PETSC_EXTERN PetscErrorCode PetscSectionGetValueLayout(MPI_Comm, PetscSection, PetscLayout *);
+PETSC_EXTERN PetscErrorCode PetscSectionPermute(PetscSection, IS, PetscSection *);
 
 PETSC_EXTERN PetscErrorCode PetscSectionSetClosureIndex(PetscSection, PetscObject, PetscSection, IS);
 PETSC_EXTERN PetscErrorCode PetscSectionGetClosureIndex(PetscSection, PetscObject, PetscSection *, IS *);

include/petscmat.h

 PETSC_EXTERN PetscErrorCode MatNestSetSubMat(Mat,PetscInt,PetscInt,Mat);
 
 PETSC_EXTERN PetscErrorCode MatChop(Mat,PetscReal);
+PETSC_EXTERN PetscErrorCode MatComputeBandwidth(Mat,PetscReal,PetscInt*);
 
 #endif

src/dm/impls/plex/examples/tests/ex10.c

+static char help[] = "Test for mesh reordering\n\n";
+
+#include <petscdmplex.h>
+
+typedef struct {
+  PetscInt  dim;               /* The topological mesh dimension */
+  PetscBool cellSimplex;       /* Flag for simplices */
+  PetscBool interpolate;       /* Flag for mesh interpolation */
+  PetscBool refinementUniform; /* Uniformly refine the mesh */
+  PetscReal refinementLimit;   /* Maximum volume of a refined cell */
+  PetscInt  numFields;         /* The number of section fields */
+  PetscInt *numComponents;     /* The number of field components */
+  PetscInt *numDof;            /* The dof signature for the section */
+} AppCtx;
+
+#undef __FUNCT__
+#define __FUNCT__ "ProcessOptions"
+PetscErrorCode ProcessOptions(AppCtx *options)
+{
+  PetscInt       len;
+  PetscBool      flg;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  options->dim               = 2;
+  options->cellSimplex       = PETSC_TRUE;
+  options->interpolate       = PETSC_FALSE;
+  options->refinementUniform = PETSC_FALSE;
+  options->refinementLimit   = 0.0;
+  options->numFields         = 1;
+  options->numComponents     = NULL;
+  options->numDof            = NULL;
+
+  ierr = PetscOptionsBegin(PETSC_COMM_SELF, "", "Meshing Problem Options", "DMPLEX");CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex10.c", options->dim, &options->dim, NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-cell_simplex", "Flag for simplices", "ex10.c", options->cellSimplex, &options->cellSimplex, NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-interpolate", "Flag for mesh interpolation", "ex10.c", options->interpolate, &options->interpolate, NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-refinement_uniform", "Uniformly refine the mesh", "ex10.c", options->refinementUniform, &options->refinementUniform, NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-refinement_limit", "The maximum volume of a refined cell", "ex10.c", options->refinementLimit, &options->refinementLimit, NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-num_fields", "The number of section fields", "ex10.c", options->numFields, &options->numFields, NULL);CHKERRQ(ierr);
+  if (options->numFields) {
+    len  = options->numFields;
+    ierr = PetscMalloc(len * sizeof(PetscInt), &options->numComponents);CHKERRQ(ierr);
+    ierr = PetscOptionsIntArray("-num_components", "The number of components per field", "ex10.c", options->numComponents, &len, &flg);CHKERRQ(ierr);
+    if (flg && (len != options->numFields)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of components array is %d should be %d", len, options->numFields);
+  }
+  len  = (options->dim+1) * PetscMax(1, options->numFields);
+  ierr = PetscMalloc(len * sizeof(PetscInt), &options->numDof);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray("-num_dof", "The dof signature for the section", "ex10.c", options->numDof, &len, &flg);CHKERRQ(ierr);
+  if (flg && (len != (options->dim+1) * PetscMax(1, options->numFields))) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Length of dof array is %d should be %d", len, (options->dim+1) * PetscMax(1, options->numFields));
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+};
+
+#undef __FUNCT__
+#define __FUNCT__ "CleanupContext"
+PetscErrorCode CleanupContext(AppCtx *user)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = PetscFree(user->numComponents);CHKERRQ(ierr);
+  ierr = PetscFree(user->numDof);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "TestReordering"
+PetscErrorCode TestReordering(DM dm, AppCtx *user)
+{
+  DM              pdm;
+  IS              perm;
+  Mat             A, pA;
+  PetscInt        bw, pbw;
+  MatOrderingType order = MATORDERINGRCM;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetOrdering(dm, order, &perm);CHKERRQ(ierr);
+  ierr = DMPlexPermute(dm, perm, &pdm);CHKERRQ(ierr);
+  ierr = DMSetFromOptions(pdm);CHKERRQ(ierr);
+  ierr = ISDestroy(&perm);CHKERRQ(ierr);
+  ierr = DMCreateMatrix(dm, MATAIJ, &A);CHKERRQ(ierr);
+  ierr = DMCreateMatrix(pdm, MATAIJ, &pA);CHKERRQ(ierr);
+  ierr = MatComputeBandwidth(A, 0.0, &bw);CHKERRQ(ierr);
+  ierr = MatComputeBandwidth(pA, 0.0, &pbw);CHKERRQ(ierr);
+  ierr = MatDestroy(&A);CHKERRQ(ierr);
+  ierr = MatDestroy(&pA);CHKERRQ(ierr);
+  ierr = DMDestroy(&pdm);CHKERRQ(ierr);
+  if (pbw > bw) {
+    ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Ordering method %s increased bandwidth from %d to %d\n", order, bw, pbw);CHKERRQ(ierr);
+  } else {
+    ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Ordering method %s reduced bandwidth from %d to %d\n", order, bw, pbw);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char **argv)
+{
+  DM             dm;
+  PetscSection   s;
+  AppCtx         user;
+  PetscErrorCode ierr;
+
+  ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
+  ierr = ProcessOptions(&user);CHKERRQ(ierr);
+  ierr = DMPlexCreateDoublet(PETSC_COMM_WORLD, user.dim, user.cellSimplex, user.interpolate, user.refinementUniform, user.refinementLimit, &dm);CHKERRQ(ierr);
+  ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
+  ierr = DMPlexCreateSection(dm, user.dim, user.numFields, user.numComponents, user.numDof, 0, NULL, NULL, &s);CHKERRQ(ierr);
+  ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);
+  ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
+  ierr = TestReordering(dm, &user);CHKERRQ(ierr);
+  ierr = DMDestroy(&dm);CHKERRQ(ierr);
+  ierr = CleanupContext(&user);CHKERRQ(ierr);
+  ierr = PetscFinalize();
+  return 0;
+}

src/dm/impls/plex/examples/tests/output/ex10_0.out

+Mat Object: 1 MPI processes
+  type: seqaij
+row 0: (0, 0)  (1, 0)  (2, 0) 
+row 1: (0, 0)  (1, 0)  (2, 0)  (3, 0) 
+row 2: (0, 0)  (1, 0)  (2, 0)  (3, 0) 
+row 3: (1, 0)  (2, 0)  (3, 0) 
+Mat Object: 1 MPI processes
+  type: seqaij
+row 0: (0, 0)  (1, 0)  (2, 0) 
+row 1: (0, 0)  (1, 0)  (2, 0)  (3, 0) 
+row 2: (0, 0)  (1, 0)  (2, 0)  (3, 0) 
+row 3: (1, 0)  (2, 0)  (3, 0) 
+Ordering method rcm reduced bandwidth from 5 to 5

src/dm/impls/plex/examples/tests/output/ex10_1.out

+Mat Object: 1 MPI processes
+  type: seqaij
+row 0: (0, 0)  (1, 0)  (2, 0)  (3, 0) 
+row 1: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0) 
+row 2: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0) 
+row 3: (0, 0)  (1, 0)  (2, 0)  (3, 0) 
+row 4: (1, 0)  (2, 0)  (4, 0)  (5, 0) 
+row 5: (1, 0)  (2, 0)  (4, 0)  (5, 0) 
+Mat Object: 1 MPI processes
+  type: seqaij
+row 0: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0) 
+row 1: (0, 0)  (1, 0)  (2, 0)  (3, 0) 
+row 2: (0, 0)  (1, 0)  (2, 0)  (3, 0) 
+row 3: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0) 
+row 4: (0, 0)  (3, 0)  (4, 0)  (5, 0) 
+row 5: (0, 0)  (3, 0)  (4, 0)  (5, 0) 
+Ordering method rcm increased bandwidth from 9 to 11

src/dm/impls/plex/examples/tests/output/ex10_2.out

+Mat Object: 1 MPI processes
+  type: seqaij
+row 0: (0, 0)  (1, 0)  (2, 0)  (3, 0) 
+row 1: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0) 
+row 2: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0) 
+row 3: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0) 
+row 4: (1, 0)  (2, 0)  (3, 0)  (4, 0) 
+Mat Object: 1 MPI processes
+  type: seqaij
+row 0: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0) 
+row 1: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0) 
+row 2: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0) 
+row 3: (0, 0)  (1, 0)  (2, 0)  (3, 0) 
+row 4: (0, 0)  (1, 0)  (2, 0)  (4, 0) 
+Ordering method rcm increased bandwidth from 7 to 9

src/dm/impls/plex/examples/tests/output/ex10_3.out

+Mat Object: 1 MPI processes
+  type: seqaij
+row 0: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0) 
+row 1: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0) 
+row 2: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 3: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 4: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0) 
+row 5: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 6: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 7: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0) 
+row 8: (2, 0)  (3, 0)  (5, 0)  (6, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 9: (2, 0)  (3, 0)  (5, 0)  (6, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 10: (2, 0)  (3, 0)  (5, 0)  (6, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 11: (2, 0)  (3, 0)  (5, 0)  (6, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+Mat Object: 1 MPI processes
+  type: seqaij
+row 0: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 1: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 2: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0) 
+row 3: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0) 
+row 4: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 5: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0) 
+row 6: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0) 
+row 7: (0, 0)  (1, 0)  (2, 0)  (3, 0)  (4, 0)  (5, 0)  (6, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 8: (0, 0)  (1, 0)  (4, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 9: (0, 0)  (1, 0)  (4, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 10: (0, 0)  (1, 0)  (4, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+row 11: (0, 0)  (1, 0)  (4, 0)  (7, 0)  (8, 0)  (9, 0)  (10, 0)  (11, 0) 
+Ordering method rcm increased bandwidth from 19 to 23

src/dm/impls/plex/examples/tests/output/ex10_4.out

+Ordering method rcm reduced bandwidth from 155 to 43

src/dm/impls/plex/examples/tests/output/ex10_5.out

+Ordering method rcm reduced bandwidth from 27 to 13

src/dm/impls/plex/examples/tests/output/ex10_6.out

+Ordering method rcm reduced bandwidth from 71 to 45

src/dm/impls/plex/makefile

 CPPFLAGS =
 CFLAGS   =
 FFLAGS   =
-SOURCEC  = plexcreate.c plex.c plexinterpolate.c plexpreallocate.c plexgeometry.c plexlabel.c plexsubmesh.c plexexodusii.c plexcgns.c plexvtk.c plexpoint.c plexvtu.c plexfem.c
+SOURCEC  = plexcreate.c plex.c plexinterpolate.c plexpreallocate.c plexreorder.c plexgeometry.c plexlabel.c plexsubmesh.c plexexodusii.c plexcgns.c plexvtk.c plexpoint.c plexvtu.c plexfem.c
 SOURCEF  =
 SOURCEH  =
 DIRS     = examples

src/dm/impls/plex/plexcreate.c

 }
 
 #undef __FUNCT__
+#define __FUNCT__ "DMPlexCreateDoublet"
+/*@
+  DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement.
+
+  Collective on MPI_Comm
+
+  Input Parameters:
++ comm - The communicator for the DM object
+. dim - The spatial dimension
+. simplex - Flag for simplicial cells, otherwise they are tensor product cells
+. interpolate - Flag to create intermediate mesh pieces (edges, faces)
+. refinementUniform - Flag for uniform parallel refinement
+- refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell
+
+  Output Parameter:
+. dm  - The DM object
+
+  Level: beginner
+
+.keywords: DM, create
+.seealso: DMSetType(), DMCreate()
+@*/
+PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm)
+{
+  DM             dm;
+  PetscInt       p;
+  PetscMPIInt    rank;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
+  ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
+  ierr = DMPlexSetDimension(dm, dim);CHKERRQ(ierr);
+  ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
+  switch (dim) {
+  case 2:
+    if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);}
+    else         {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);}
+    break;
+  case 3:
+    if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);}
+    else         {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);}
+    break;
+  default:
+    SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
+  }
+  if (rank) {
+    PetscInt numPoints[2] = {0, 0};
+    ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
+  } else {
+    switch (dim) {
+    case 2:
+      if (simplex) {
+        PetscInt    numPoints[2]        = {4, 2};
+        PetscInt    coneSize[6]         = {3, 3, 0, 0, 0, 0};
+        PetscInt    cones[6]            = {2, 3, 4,  5, 4, 3};
+        PetscInt    coneOrientations[6] = {0, 0, 0,  0, 0, 0};
+        PetscScalar vertexCoords[8]     = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
+        PetscInt    markerPoints[8]     = {2, 1, 3, 1, 4, 1, 5, 1};
+
+        ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
+        for (p = 0; p < 4; ++p) {ierr = DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
+      } else {
+        PetscInt    numPoints[2]        = {6, 2};
+        PetscInt    coneSize[8]         = {4, 4, 0, 0, 0, 0, 0, 0};
+        PetscInt    cones[8]            = {2, 3, 4, 5,  3, 6, 7, 4};
+        PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
+        PetscScalar vertexCoords[12]    = {-1.0, -0.5,  0.0, -0.5,  0.0, 0.5,  -1.0, 0.5,  1.0, -0.5,  1.0, 0.5};
+
+        ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
+      }
+      break;
+    case 3:
+      if (simplex) {
+        PetscInt    numPoints[2]        = {5, 2};
+        PetscInt    coneSize[7]         = {4, 4, 0, 0, 0, 0, 0};
+        PetscInt    cones[8]            = {4, 3, 5, 2,  5, 3, 4, 6};
+        PetscInt    coneOrientations[8] = {0, 0, 0, 0,  0, 0, 0, 0};
+        PetscScalar vertexCoords[15]    = {-1.0, 0.0, 0.0,  0.0, -1.0, 0.0,  0.0, 0.0, 1.0,  0.0, 1.0, 0.0,  1.0, 0.0, 0.0};
+        PetscInt    markerPoints[10]    = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};
+
+        ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
+        for (p = 0; p < 5; ++p) {ierr = DMPlexSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
+      } else {
+        PetscInt    numPoints[2]         = {12, 2};
+        PetscInt    coneSize[14]         = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+        PetscInt    cones[16]            = {2, 3, 4, 5, 6, 7, 8, 9,  5, 4, 10, 11, 7, 12, 13, 8};
+        PetscInt    coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0,  0, 0,  0,  0, 0,  0,  0, 0};
+        PetscScalar vertexCoords[36]     = {-1.0, -0.5, -0.5,  -1.0,  0.5, -0.5,  0.0,  0.5, -0.5,   0.0, -0.5, -0.5,
+                                            -1.0, -0.5,  0.5,   0.0, -0.5,  0.5,  0.0,  0.5,  0.5,  -1.0,  0.5,  0.5,
+                                             1.0,  0.5, -0.5,   1.0, -0.5, -0.5,  1.0, -0.5,  0.5,   1.0,  0.5,  0.5};
+
+        ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
+      }
+      break;
+    default:
+      SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %d", dim);
+    }
+  }
+  *newdm = dm;
+  if (refinementLimit > 0.0) {
+    DM rdm;
+    const char *name;
+
+    ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr);
+    ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr);
+    ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr);
+    ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr);
+    ierr = PetscObjectSetName((PetscObject)    rdm,  name);CHKERRQ(ierr);
+    ierr = DMDestroy(newdm);CHKERRQ(ierr);
+    *newdm = rdm;
+  }
+  if (interpolate) {
+    DM idm;
+    const char *name;
+
+    ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr);
+    ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr);
+    ierr = PetscObjectSetName((PetscObject)    idm,  name);CHKERRQ(ierr);
+    ierr = DMPlexCopyCoordinates(*newdm, idm);CHKERRQ(ierr);
+    ierr = DMPlexCopyLabels(*newdm, idm);CHKERRQ(ierr);
+    ierr = DMDestroy(newdm);CHKERRQ(ierr);
+    *newdm = idm;
+  }
+  {
+    DM refinedMesh     = NULL;
+    DM distributedMesh = NULL;
+
+    /* Distribute mesh over processes */
+    ierr = DMPlexDistribute(*newdm, NULL, 0, NULL, &distributedMesh);CHKERRQ(ierr);
+    if (distributedMesh) {
+      ierr = DMDestroy(newdm);CHKERRQ(ierr);
+      *newdm = distributedMesh;
+    }
+    if (refinementUniform) {
+      ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr);
+      ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr);
+      if (refinedMesh) {
+        ierr = DMDestroy(newdm);CHKERRQ(ierr);
+        *newdm = refinedMesh;
+      }
+    }
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "DMPlexCreateSquareBoundary"
 /*
  Simple square boundary:

src/dm/impls/plex/plexinterpolate.c

   Note: This is typically used when adding pieces other than vertices to a mesh
 
 .keywords: mesh
-.seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMPlexGetCoordinateSection()
+.seealso: DMCopyLabels(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMPlexGetCoordinateSection()
 @*/
 PetscErrorCode DMPlexCopyCoordinates(DM dmA, DM dmB)
 {
   ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexCopyLabels"
+/*@
+  DMPlexCopyLabels - Copy labels from one mesh to another with a superset of the points
+
+  Collective on DM
+
+  Input Parameter:
+. dmA - The DMPlex object with initial labels
+
+  Output Parameter:
+. dmB - The DMPlex object with copied labels
+
+  Level: intermediate
+
+  Note: This is typically used when interpolating or otherwise adding to a mesh
+
+.keywords: mesh
+.seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMPlexGetCoordinateSection()
+@*/
+PetscErrorCode DMPlexCopyLabels(DM dmA, DM dmB)
+{
+  PetscInt       numLabels, l;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  if (dmA == dmB) PetscFunctionReturn(0);
+  ierr = DMPlexGetNumLabels(dmA, &numLabels);CHKERRQ(ierr);
+  for (l = 0; l < numLabels; ++l) {
+    DMLabel     label, labelNew;
+    const char *name;
+    PetscBool   flg;
+
+    ierr = DMPlexGetLabelName(dmA, l, &name);CHKERRQ(ierr);
+    ierr = PetscStrcmp(name, "depth", &flg);CHKERRQ(ierr);
+    if (flg) continue;
+    ierr = DMPlexGetLabel(dmA, name, &label);CHKERRQ(ierr);
+    ierr = DMLabelDuplicate(label, &labelNew);CHKERRQ(ierr);
+    ierr = DMPlexAddLabel(dmB, labelNew);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}

src/dm/impls/plex/plexlabel.c

   PetscFunctionReturn(0);
 }
 
+#undef __FUNCT__
+#define __FUNCT__ "DMLabelPermute"
+PetscErrorCode DMLabelPermute(DMLabel label, IS permutation, DMLabel *labelNew)
+{
+  const PetscInt *perm;
+  PetscInt        numValues, numPoints, v, q;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  ierr = DMLabelDuplicate(label, labelNew);CHKERRQ(ierr);
+  ierr = DMLabelGetNumValues(*labelNew, &numValues);CHKERRQ(ierr);
+  ierr = ISGetLocalSize(permutation, &numPoints);CHKERRQ(ierr);
+  ierr = ISGetIndices(permutation, &perm);CHKERRQ(ierr);
+  for (v = 0; v < numValues; ++v) {
+    const PetscInt offset = (*labelNew)->stratumOffsets[v];
+    const PetscInt size   = (*labelNew)->stratumSizes[v];
+
+    for (q = offset; q < offset+size; ++q) {
+      const PetscInt point = (*labelNew)->points[q];
+
+      if ((point < 0) || (point >= numPoints)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label point %d is not in [0, %d) for the remapping", point, numPoints);
+      (*labelNew)->points[q] = perm[point];
+    }
+    ierr = PetscSortInt(size, &(*labelNew)->points[offset]);CHKERRQ(ierr);
+  }
+  ierr = ISRestoreIndices(permutation, &perm);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
 
 
 #undef __FUNCT__

src/dm/impls/plex/plexreorder.c

+#include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
+#include <petsc-private/matorderimpl.h> /*I      "petscmat.h"      I*/
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexCreateOrderingClosure_Static"
+PetscErrorCode DMPlexCreateOrderingClosure_Static(DM dm, PetscInt numPoints, const PetscInt pperm[], PetscInt **clperm, PetscInt **invclperm)
+{
+  PetscInt      *perm, *iperm;
+  PetscInt       depth, d, pStart, pEnd, fStart, fMax, fEnd, p;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
+  ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
+  ierr = PetscMalloc2(pEnd-pStart,PetscInt,&perm,pEnd-pStart,PetscInt,&iperm);CHKERRQ(ierr);
+  for (p = pStart; p < pEnd; ++p) iperm[p] = -1;
+  for (d = depth; d > 0; --d) {
+    ierr = DMPlexGetDepthStratum(dm, d,   &pStart, &pEnd);CHKERRQ(ierr);
+    ierr = DMPlexGetDepthStratum(dm, d-1, &fStart, &fEnd);CHKERRQ(ierr);
+    fMax = fStart;
+    for (p = pStart; p < pEnd; ++p) {
+      const PetscInt *cone;
+      PetscInt        point, coneSize, c;
+
+      if (d == depth) {
+        perm[p]         = pperm[p];
+        iperm[pperm[p]] = p;
+      }
+      point = perm[p];
+      ierr = DMPlexGetConeSize(dm, point, &coneSize);CHKERRQ(ierr);
+      ierr = DMPlexGetCone(dm, point, &cone);CHKERRQ(ierr);
+      for (c = 0; c < coneSize; ++c) {
+        const PetscInt oldc = cone[c];
+        const PetscInt newc = iperm[oldc];
+
+        if (newc < 0) {
+          perm[fMax]  = oldc;
+          iperm[oldc] = fMax++;
+        }
+      }
+    }
+    if (fMax != fEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of depth %d faces %d does not match permuted nubmer %d", d, fEnd-fStart, fMax-fStart);
+  }
+  *clperm    = perm;
+  *invclperm = iperm;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexGetOrdering"
+/*@
+  DMPlexGetOrdering - Calculate a reordering of the mesh
+
+  Collective on DM
+
+  Input Parameter:
++ dm - The DMPlex object
+- otype - type of reordering, one of the following:
+$     MATORDERINGNATURAL - Natural
+$     MATORDERINGND - Nested Dissection
+$     MATORDERING1WD - One-way Dissection
+$     MATORDERINGRCM - Reverse Cuthill-McKee
+$     MATORDERINGQMD - Quotient Minimum Degree
+
+
+  Output Parameter:
+. perm - The point permutation as an IS
+
+  Level: intermediate
+
+.keywords: mesh
+.seealso: MatGetOrdering()
+@*/
+PetscErrorCode DMPlexGetOrdering(DM dm, MatOrderingType otype, IS *perm)
+{
+  PetscInt       numCells = 0;
+  PetscInt      *start = NULL, *adjacency = NULL, *cperm, *clperm, *invclperm, *mask, *xls, pStart, pEnd, c, i;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  PetscValidPointer(perm, 3);
+  ierr = DMPlexCreateNeighborCSR(dm, 0, &numCells, &start, &adjacency);CHKERRQ(ierr);
+  ierr = PetscMalloc3(numCells,PetscInt,&cperm,numCells,PetscInt,&mask,numCells*2,PetscInt,&xls);CHKERRQ(ierr);
+  /* Shift for Fortran numbering */
+  for (i = 0; i < start[numCells]; ++i) ++adjacency[i];
+  for (i = 0; i <= numCells; ++i)       ++start[i];
+  ierr = SPARSEPACKgenrcm(&numCells, start, adjacency, cperm, mask, xls);CHKERRQ(ierr);
+  ierr = PetscFree(start);CHKERRQ(ierr);
+  ierr = PetscFree(adjacency);CHKERRQ(ierr);
+  /* Shift for Fortran numbering */
+  for (c = 0; c < numCells; ++c) --cperm[c];
+  /* Construct closure */
+  ierr = DMPlexCreateOrderingClosure_Static(dm, numCells, cperm, &clperm, &invclperm);
+  ierr = PetscFree3(cperm,mask,xls);CHKERRQ(ierr);
+  ierr = PetscFree(clperm);CHKERRQ(ierr);
+  /* Invert permutation */
+  ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
+  ierr = ISCreateGeneral(PetscObjectComm((PetscObject) dm), pEnd-pStart, invclperm, PETSC_OWN_POINTER, perm);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMPlexPermute"
+/*@
+  DMPlexPermute - Reorder the mesh according to the input permutation
+
+  Collective on DM
+
+  Input Parameter:
++ dm - The DMPlex object
+- perm - The point permutation
+
+  Output Parameter:
+. pdm - The permuted DM
+
+  Level: intermediate
+
+.keywords: mesh
+.seealso: MatPermute()
+@*/
+PetscErrorCode DMPlexPermute(DM dm, IS perm, DM *pdm)
+{
+  DM_Plex       *plex = (DM_Plex *) dm->data, *plexNew;
+  PetscSection   section, sectionNew;
+  PetscInt       dim;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
+  PetscValidHeaderSpecific(perm, IS_CLASSID, 2);
+  PetscValidPointer(pdm, 3);
+  ierr = DMCreate(PetscObjectComm((PetscObject) dm), pdm);CHKERRQ(ierr);
+  ierr = DMSetType(*pdm, DMPLEX);CHKERRQ(ierr);
+  ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
+  ierr = DMPlexSetDimension(*pdm, dim);CHKERRQ(ierr);
+  ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
+  if (section) {
+    ierr = PetscSectionPermute(section, perm, &sectionNew);CHKERRQ(ierr);
+    ierr = DMSetDefaultSection(*pdm, sectionNew);CHKERRQ(ierr);
+    ierr = PetscSectionDestroy(&sectionNew);CHKERRQ(ierr);
+  }
+  plexNew = (DM_Plex *) (*pdm)->data;
+  /* Ignore ltogmap, ltogmapb */
+  /* Ignore sf, defaultSF */
+  /* Ignore globalVertexNumbers, globalCellNumbers */
+  /* Remap coordinates */
+  {
+    DM           cdm, cdmNew;
+    PetscSection csection, csectionNew;
+    Vec          coordinates, coordinatesNew;
+    PetscScalar *coords, *coordsNew;
+    PetscInt     pStart, pEnd, p;
+
+    ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
+    ierr = DMGetDefaultSection(cdm, &csection);CHKERRQ(ierr);
+    ierr = PetscSectionPermute(csection, perm, &csectionNew);CHKERRQ(ierr);
+    ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
+    ierr = VecDuplicate(coordinates, &coordinatesNew);CHKERRQ(ierr);
+    ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
+    ierr = VecGetArray(coordinatesNew, &coordsNew);CHKERRQ(ierr);
+    ierr = PetscSectionGetChart(csectionNew, &pStart, &pEnd);CHKERRQ(ierr);
+    for (p = pStart; p < pEnd; ++p) {
+      PetscInt dof, off, offNew, d;
+
+      ierr = PetscSectionGetDof(csectionNew, p, &dof);CHKERRQ(ierr);
+      ierr = PetscSectionGetOffset(csection, p, &off);CHKERRQ(ierr);
+      ierr = PetscSectionGetOffset(csectionNew, p, &offNew);CHKERRQ(ierr);
+      for (d = 0; d < dof; ++d) coordsNew[offNew+d] = coords[off+d];
+    }
+    ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
+    ierr = VecRestoreArray(coordinatesNew, &coordsNew);CHKERRQ(ierr);
+    ierr = DMGetCoordinateDM(*pdm, &cdmNew);CHKERRQ(ierr);
+    ierr = DMSetDefaultSection(cdmNew, csectionNew);CHKERRQ(ierr);
+    ierr = DMSetCoordinatesLocal(*pdm, coordinatesNew);CHKERRQ(ierr);
+    ierr = PetscSectionDestroy(&csectionNew);CHKERRQ(ierr);
+    ierr = VecDestroy(&coordinatesNew);CHKERRQ(ierr);
+  }
+  /* Reorder labels */
+  {
+    DMLabel label = plex->labels, labelNew = NULL;
+
+    while (label) {
+      if (!plexNew->labels) {
+        ierr = DMLabelPermute(label, perm, &plexNew->labels);CHKERRQ(ierr);
+        labelNew = plexNew->labels;
+      } else {
+        ierr = DMLabelPermute(label, perm, &labelNew->next);CHKERRQ(ierr);
+        labelNew = labelNew->next;
+      }
+      label = label->next;
+    }
+    if (plex->subpointMap) {ierr = DMLabelPermute(plex->subpointMap, perm, &plexNew->subpointMap);CHKERRQ(ierr);}
+  }
+  /* Reorder topology */
+  {
+    const PetscInt *pperm;
+    PetscInt        maxConeSize, maxSupportSize, n, pStart, pEnd, p;
+
+    ierr = DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);CHKERRQ(ierr);
+    plexNew->maxConeSize    = maxConeSize;
+    plexNew->maxSupportSize = maxSupportSize;
+    ierr = PetscSectionDestroy(&plexNew->coneSection);CHKERRQ(ierr);
+    ierr = PetscSectionPermute(plex->coneSection, perm, &plexNew->coneSection);CHKERRQ(ierr);
+    ierr = PetscSectionGetStorageSize(plexNew->coneSection, &n);CHKERRQ(ierr);
+    ierr = PetscMalloc(n * sizeof(PetscInt), &plexNew->cones);CHKERRQ(ierr);
+    ierr = PetscMalloc(n * sizeof(PetscInt), &plexNew->coneOrientations);CHKERRQ(ierr);
+    ierr = ISGetIndices(perm, &pperm);CHKERRQ(ierr);
+    ierr = PetscSectionGetChart(plex->coneSection, &pStart, &pEnd);CHKERRQ(ierr);
+    for (p = pStart; p < pEnd; ++p) {
+      PetscInt dof, off, offNew, d;
+
+      ierr = PetscSectionGetDof(plexNew->coneSection, pperm[p], &dof);CHKERRQ(ierr);
+      ierr = PetscSectionGetOffset(plex->coneSection, p, &off);CHKERRQ(ierr);
+      ierr = PetscSectionGetOffset(plexNew->coneSection, pperm[p], &offNew);CHKERRQ(ierr);
+      for (d = 0; d < dof; ++d) {
+        plexNew->cones[offNew+d]            = pperm[plex->cones[off+d]];
+        plexNew->coneOrientations[offNew+d] = plex->coneOrientations[off+d];
+      }
+    }
+    ierr = PetscSectionDestroy(&plexNew->supportSection);CHKERRQ(ierr);
+    ierr = PetscSectionPermute(plex->supportSection, perm, &plexNew->supportSection);CHKERRQ(ierr);
+    ierr = PetscSectionGetStorageSize(plexNew->supportSection, &n);CHKERRQ(ierr);
+    ierr = PetscMalloc(n * sizeof(PetscInt), &plexNew->supports);CHKERRQ(ierr);
+    ierr = PetscSectionGetChart(plex->supportSection, &pStart, &pEnd);CHKERRQ(ierr);
+    for (p = pStart; p < pEnd; ++p) {
+      PetscInt dof, off, offNew, d;
+
+      ierr = PetscSectionGetDof(plexNew->supportSection, pperm[p], &dof);CHKERRQ(ierr);
+      ierr = PetscSectionGetOffset(plex->supportSection, p, &off);CHKERRQ(ierr);
+      ierr = PetscSectionGetOffset(plexNew->supportSection, pperm[p], &offNew);CHKERRQ(ierr);
+      for (d = 0; d < dof; ++d) {
+        plexNew->supports[offNew+d] = pperm[plex->supports[off+d]];
+      }
+    }
+    ierr = ISRestoreIndices(perm, &pperm);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}

src/mat/order/amd/amd.c

 
 #include <petscmat.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 #define UF_long long long
 #define UF_long_max LONG_LONG_MAX
 #define UF_long_id "%lld"

src/mat/order/degree.c

 
 /* degree.f -- translated by f2c (version 19931217).*/
 
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*****************************************************************/
 /*********     DEGREE ..... DEGREE IN MASKED COMPONENT   *********/

src/mat/order/fn1wd.c

 
 /* fn1wd.f -- translated by f2c (version 19931217).*/
 
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*****************************************************************/
 /********     FN1WD ..... FIND ONE-WAY DISSECTORS        *********/

src/mat/order/fndsep.c

 /* fndsep.f -- translated by f2c (version 19931217).
 */
 
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*****************************************************************/
 /*************     FNDSEP ..... FIND SEPARATOR       *************/

src/mat/order/fnroot.c

 /* fnroot.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*****************************************************************/
 /********     FNROOT ..... FIND PSEUDO-PERIPHERAL NODE    ********/

src/mat/order/gen1wd.c

 /* gen1wd.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*****************************************************************/
 /***********     GEN1WD ..... GENERAL ONE-WAY DISSECTION  ********/

src/mat/order/gennd.c

 /* gennd.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 #undef __FUNCT__
 #define __FUNCT__ "SPARSEPACKrevrse"

src/mat/order/genqmd.c

 /* genqmd.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /******************************************************************/
 /***********    GENQMD ..... QUOT MIN DEGREE ORDERING    **********/

src/mat/order/genrcm.c

 /* genrcm.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*****************************************************************/
 /*****************************************************************/

src/mat/order/makefile

 SOURCEC  = sp1wd.c spnd.c spqmd.c sprcm.c sorder.c sregis.c\
            degree.c  fnroot.c genqmd.c qmdqt.c rcm.c fn1wd.c gen1wd.c \
            genrcm.c qmdrch.c rootls.c fndsep.c gennd.c qmdmrg.c qmdupd.c
-SOURCEH  = order.h
+SOURCEH  = ../../../include/petsc-private/matorderimpl.h
 LIBBASE  = libpetscmat
 DIRS     = amd
 LOCDIR   = src/mat/order/

src/mat/order/order.h

-#include <petscmat.h>
-#include <petsc-private/petscimpl.h>
-/*
-   Defines the interface to the SparsePack routines, translated into C.
-*/
-PETSC_EXTERN PetscErrorCode SPARSEPACKgen1wd(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKgennd(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKgenrcm(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKgenqmd(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-
-PETSC_EXTERN PetscErrorCode SPARSEPACKqmdrch(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKqmdmrg(const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKqmdqt(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKqmdupd(const PetscInt*,const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKfnroot(PetscInt*,const PetscInt*,const PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKrootls(const PetscInt*,const PetscInt*,const PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKfn1wd(PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKrevrse(const PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKrootls(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKfndsep(PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKdegree(const PetscInt*,const PetscInt*,const PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscInt*);
-PETSC_EXTERN PetscErrorCode SPARSEPACKrcm(const PetscInt*,const PetscInt*,const PetscInt*,PetscInt*,PetscInt*,PetscInt*,PetscInt*);
-
-
-

src/mat/order/qmdmrg.c

 /* qmdmrg.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /******************************************************************/
 /***********     QMDMRG ..... QUOT MIN DEG MERGE       ************/

src/mat/order/qmdqt.c

 /* qmdqt.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /***************************************************************/
 /********     QMDQT  ..... QUOT MIN DEG QUOT TRANSFORM  ********/

src/mat/order/qmdrch.c

 /* qmdrch.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*****************************************************************/
 /**********     QMDRCH ..... QUOT MIN DEG REACH SET    ***********/

src/mat/order/qmdupd.c

 /* qmdupd.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /******************************************************************/
 /***********     QMDUPD ..... QUOT MIN DEG UPDATE      ************/

src/mat/order/rcm.c

 /* rcm.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*****************************************************************/
 /*********     RCM ..... REVERSE CUTHILL-MCKEE ORDERING   *******/

src/mat/order/rootls.c

 /* rootls.f -- translated by f2c (version 19931217).*/
 
 #include <petscsys.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*****************************************************************/
 /*********     ROOTLS ..... ROOTED LEVEL STRUCTURE      **********/

src/mat/order/sorder.c

 
   PetscFunctionBegin;
   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
-  PetscValidPointer(rperm,2);
-  PetscValidPointer(cperm,3);
+  PetscValidPointer(rperm,3);
+  PetscValidPointer(cperm,4);
   if (!mat->assembled) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
   if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
 

src/mat/order/sp1wd.c

 
 #include <petscmat.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*
     MatGetOrdering_1WD - Find the 1-way dissection ordering of a given matrix.

src/mat/order/spnd.c

 
 #include <petscmat.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*
     MatGetOrdering_ND - Find the nested dissection ordering of a given matrix.

src/mat/order/spqmd.c

 
 #include <petscmat.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*
     MatGetOrdering_QMD - Find the Quotient Minimum Degree ordering of a given matrix.

src/mat/order/sprcm.c

 
 #include <petscmat.h>
-#include <../src/mat/order/order.h>
+#include <petsc-private/matorderimpl.h>
 
 /*
     MatGetOrdering_RCM - Find the Reverse Cuthill-McKee ordering of a given matrix.

src/mat/utils/bandwidth.c

+#include <petsc-private/matimpl.h>       /*I  "petscmat.h"  I*/
+
+#undef __FUNCT__
+#define __FUNCT__ "MatComputeBandwidth"
+/*@
+  MatComputeBandwidth - Calculate the full bandwidth of the matrix, meaning the width 2k+1 where k diagonals on either side are sufficient to contain all the matrix nonzeros.
+
+  Collective on Mat
+
+  Input Parameters:
++ A - The Mat
+- fraction - An optional percentage of the Frobenius norm of the matrix that the bandwidth should enclose
+
+  Output Parameter:
+. bw - The matrix bandwidth
+
+  Level: beginner
+
+.seealso: DMPlexCreate(), DMPlexSetConeSize(), DMPlexSetChart()
+@*/
+PetscErrorCode MatComputeBandwidth(Mat A, PetscReal fraction, PetscInt *bw)
+{
+  PetscInt       lbw[2] = {0, 0}, gbw[2];
+  PetscInt       rStart, rEnd, r;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
+  PetscValidLogicalCollectiveReal(A,fraction,2);
+  PetscValidPointer(bw, 3);
+  if ((fraction > 0.0) && (fraction < 1.0)) SETERRQ(PetscObjectComm((PetscObject) A), PETSC_ERR_SUP, "We do not yet support a fractional bandwidth");
+  ierr = MatGetOwnershipRange(A, &rStart, &rEnd);CHKERRQ(ierr);
+  for (r = rStart; r < rEnd; ++r) {
+    const PetscInt *cols;
+    PetscInt        ncols;
+
+    ierr = MatGetRow(A, r, &ncols, &cols, NULL);CHKERRQ(ierr);
+    if (ncols) {
+      lbw[0] = PetscMax(lbw[0], r - cols[0]);
+      lbw[1] = PetscMax(lbw[1], cols[ncols-1] - r);
+    }
+    ierr = MatRestoreRow(A, r, &ncols, &cols, NULL);CHKERRQ(ierr);
+  }
+  ierr = MPI_Allreduce(lbw, gbw, 2, MPIU_INT, MPIU_MAX, PetscObjectComm((PetscObject) A));CHKERRQ(ierr);
+  *bw = 2*PetscMax(gbw[0], gbw[1]) + 1;
+  PetscFunctionReturn(0);
+}

src/mat/utils/makefile

 FFLAGS   =
 SOURCEC  = convert.c matstash.c axpy.c zerodiag.c \
            getcolv.c gcreate.c freespace.c compressedrow.c multequal.c \
-           matstashspace.c pheap.c
+           matstashspace.c pheap.c bandwidth.c
 SOURCEF  =
 SOURCEH  = freespace.h petscheap.h
 LIBBASE  = libpetscmat

src/vec/is/utils/vsectionis.c

   PetscFunctionReturn(0);
 }
 
+#undef __FUNCT__
+#define __FUNCT__ "PetscSectionPermute"
+/*@
+  PetscSectionPermute - Reorder the section according to the input point permutation
+
+  Collective on PetscSection
+
+  Input Parameter:
++ section - The PetscSection object
+- perm - The point permutation
+
+  Output Parameter:
+. sectionNew - The permuted PetscSection
+
+  Level: intermediate
+
+.keywords: mesh
+.seealso: MatPermute()
+@*/
+PetscErrorCode PetscSectionPermute(PetscSection section, IS permutation, PetscSection *sectionNew)
+{
+  PetscSection    s = section, sNew;
+  const PetscInt *perm;
+  PetscInt        numFields, f, numPoints, pStart, pEnd, p;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 1);
+  PetscValidHeaderSpecific(permutation, IS_CLASSID, 2);
+  PetscValidPointer(sectionNew, 3);
+  ierr = PetscSectionCreate(PetscObjectComm((PetscObject) s), &sNew);CHKERRQ(ierr);
+  ierr = PetscSectionGetNumFields(s, &numFields);CHKERRQ(ierr);
+  if (numFields) {ierr = PetscSectionSetNumFields(sNew, numFields);CHKERRQ(ierr);}
+  for (f = 0; f < numFields; ++f) {
+    const char *name;
+    PetscInt    numComp;
+
+    ierr = PetscSectionGetFieldName(s, f, &name);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldName(sNew, f, name);CHKERRQ(ierr);
+    ierr = PetscSectionGetFieldComponents(s, f, &numComp);CHKERRQ(ierr);
+    ierr = PetscSectionSetFieldComponents(sNew, f, numComp);CHKERRQ(ierr);
+  }
+  ierr = ISGetLocalSize(permutation, &numPoints);CHKERRQ(ierr);
+  ierr = ISGetIndices(permutation, &perm);CHKERRQ(ierr);
+  ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
+  ierr = PetscSectionSetChart(sNew, pStart, pEnd);CHKERRQ(ierr);
+  if (numPoints < pEnd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Permutation size %d is less than largest Section point %d", numPoints, pEnd);
+  for (p = pStart; p < pEnd; ++p) {
+    PetscInt dof, cdof;
+
+    ierr = PetscSectionGetDof(s, p, &dof);CHKERRQ(ierr);
+    ierr = PetscSectionSetDof(sNew, perm[p], dof);CHKERRQ(ierr);
+    ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr);
+    if (cdof) {ierr = PetscSectionSetConstraintDof(sNew, perm[p], cdof);CHKERRQ(ierr);}
+    for (f = 0; f < numFields; ++f) {
+      ierr = PetscSectionGetFieldDof(s, p, f, &dof);CHKERRQ(ierr);
+      ierr = PetscSectionSetFieldDof(sNew, perm[p], f, dof);CHKERRQ(ierr);
+      ierr = PetscSectionGetFieldConstraintDof(s, p, f, &cdof);CHKERRQ(ierr);
+      if (cdof) {ierr = PetscSectionSetFieldConstraintDof(sNew, perm[p], f, cdof);CHKERRQ(ierr);}
+    }
+  }
+  ierr = PetscSectionSetUp(sNew);CHKERRQ(ierr);
+  for (p = pStart; p < pEnd; ++p) {
+    const PetscInt *cind;
+    PetscInt        cdof;
+
+    ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr);
+    if (cdof) {
+      ierr = PetscSectionGetConstraintIndices(s, p, &cind);CHKERRQ(ierr);
+      ierr = PetscSectionSetConstraintIndices(sNew, perm[p], cind);CHKERRQ(ierr);
+    }
+    for (f = 0; f < numFields; ++f) {
+      ierr = PetscSectionGetFieldConstraintDof(s, p, f, &cdof);CHKERRQ(ierr);
+      if (cdof) {
+        ierr = PetscSectionGetFieldConstraintIndices(s, p, f, &cind);CHKERRQ(ierr);
+        ierr = PetscSectionSetFieldConstraintIndices(sNew, perm[p], f, cind);CHKERRQ(ierr);
+      }
+    }
+  }
+  ierr = ISRestoreIndices(permutation, &perm);CHKERRQ(ierr);
+  *sectionNew = sNew;
+  PetscFunctionReturn(0);
+}
+
 /*
   I need extract and merge routines for section based on fields. This should be trivial except for updating the
   constraint indices.
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.