Add -dfs_assemble_reduced, uncouples components in matrix assembly.

This is particularly useful to improve AMG performance.  ML and
BoomerAMG don't have block matrix formats so inter-component coupling
significantly adds to the cost.  PETSc has block smoothers and a more
efficient storage scheme so these extra entries are normally very cheap
to leave in, even if the coupling is weak (leaving them in trades
metadata for useful data).  It would be really nice if the AMG packages
could just provide interpolation matrices, respecting the nodal
structure, and let PETSc handle the rest.

include/private/fsimpl.h

@@ -58,6 +58,7 @@ struct _p_dFS {
   dMeshTag     bdyConstraintTag; /**< User-defined context for enforcing boundary constraints */
   dTruth       spacebuilt;
   dTruth       assemblefull;    /**< Use full order constraints for assembly */
+  dTruth       assemblereduced; /**< Assemble only diagonal part of blocks, only matters for bs>1 and MATAIJ */
   dInt         ruleStrength;
   dInt         bs;              /**< Block size (number of dofs per node) */
   dMeshEH     *ents;            /**< All entities in active set */
@@ -82,6 +83,7 @@ struct _p_dFS {
   Vec          dcache;          /**< All Dirichlet values, this is only a cache so that we can project a vector into the inhomogeneous space */
   VecScatter   dscat;           /**< Scatter from global closure to \a dcache. */
   ISLocalToGlobalMapping bmapping; /**< Block mapping, Dirichlet blocks have negative global index */
+  ISLocalToGlobalMapping mapping;  /**< Scalar mapping, mapping[i] = bmapping[i/bs]*bs+i%bs; */
   dInt         maxQ;
   dFSRotation  rot;             /**< Rotation for local vector */
   s_dFSWorkspace workspace[dFS_MAX_WORKSPACES];

src/fs/impls/cont/cont.c

@@ -311,6 +311,13 @@ static dErr dFSBuildSpace_Cont(dFS fs)
     err = VecDestroy(g);dCHK(err);
     err = ISLocalToGlobalMappingCreateNC(((dObject)fs)->comm,nc+ngh,globals,&fs->bmapping);dCHK(err);
     /* Don't free \a globals because we used the no-copy variant, so the IS takes ownership. */
+    {
+      dInt *sglobals;        /* Scalar globals */
+      err = dMallocA((nc+ngh)*bs,&sglobals);dCHK(err);
+      for (i=0; i<(nc+ngh)*bs; i++) sglobals[i] = globals[i/bs]*bs + i%bs;
+      err = ISLocalToGlobalMappingCreateNC(((dObject)fs)->comm,(nc+ngh)*bs,sglobals,&fs->mapping);dCHK(err);
+      /* Don't free \a sglobals because we used the no-copy variant, so the IS takes ownership. */
+    }
   }
 
   /* Create a cache for Dirichlet part of closure vector, this could be local but it doesn't cost anything to make it global */

src/fs/interface/fs.c

@@ -186,6 +186,7 @@ dErr dFSDestroy(dFS fs)
   err = MatDestroy(fs->E);dCHK(err);
   err = MatDestroy(fs->Ep);dCHK(err);
   err = ISLocalToGlobalMappingDestroy(fs->bmapping);dCHK(err);
+  err = ISLocalToGlobalMappingDestroy(fs->mapping);dCHK(err);
   err = dFree3(fs->rule,fs->efs,fs->off);dCHK(err);
   err = dMeshRestoreVertexCoords(fs->mesh,fs->nelem,NULL,&fs->vtxoff,&fs->vtx);dCHK(err);
   err = PetscHeaderDestroy(fs);dCHK(err);
@@ -560,11 +561,17 @@ dErr dFSGetMatrix(dFS fs,const MatType mtype,Mat *inJ)
   err = MatSetType(J,mtype);dCHK(err);
   err = MatSeqBAIJSetPreallocation(J,bs,27,NULL);dCHK(err);         /* \bug incorrect for unstructured meshes */
   err = MatMPIBAIJSetPreallocation(J,bs,27,NULL,25,NULL);dCHK(err); /* \todo this wastes a lot of space in parallel */
-  err = MatSeqAIJSetPreallocation(J,bs*27,NULL);dCHK(err);
-  err = MatMPIAIJSetPreallocation(J,bs*27,NULL,bs*25,NULL);dCHK(err);
+  if (fs->assemblereduced) {
+    err = MatSeqAIJSetPreallocation(J,27,NULL);dCHK(err);
+    err = MatMPIAIJSetPreallocation(J,27,NULL,25,NULL);dCHK(err);
+  } else {
+    err = MatSeqAIJSetPreallocation(J,bs*27,NULL);dCHK(err);
+    err = MatMPIAIJSetPreallocation(J,bs*27,NULL,bs*25,NULL);dCHK(err);
+  }
   err = MatHasOperation(J,MATOP_SET_BLOCK_SIZE,&hassetbs);dCHK(err);
   if (hassetbs) {err = MatSetBlockSize(J,bs);dCHK(err);}
   err = MatSetLocalToGlobalMappingBlock(J,fs->bmapping);dCHK(err);
+  err = MatSetLocalToGlobalMapping(J,fs->mapping);dCHK(err);
 
   /* We want the resulting matrices to be usable with matrix-free operations based on this FS */
   err = PetscObjectCompose((dObject)J,"DohpFS",(dObject)fs);dCHK(err);
@@ -647,7 +654,22 @@ dErr dFSMatSetValuesBlockedExpanded(dFS fs,Mat A,dInt m,const dInt idxm[],dInt n
   err = MatRestoreArray_SeqAIJ(E,&ca);dCHK(err);
   err = MatRestoreRowIJ(E,0,dFALSE,dFALSE,&cn,&ci,&cj,&done);dCHK(err);
   if (!done) dERROR(1,"Failed to return indices");
-  err = MatSetValuesBlockedLocal(A,lm,lidxm,ln,lidxn,lv,imode);dCHK(err);
+  if (fs->assemblereduced) {
+    dInt brow[lm],bcol[ln];
+    dScalar bval[lm*ln];
+    for (dInt k=0; k<bs; k++) {
+      for (i=0; i<lm; i++) {
+        for (j=0; j<ln; j++) {
+          bval[i*ln+j] = lv[(i*bs+k)*ln*bs+(j*bs+k)];
+        }
+      }
+      for (i=0; i<lm; i++) brow[i] = lidxm[i]*bs+k;
+      for (j=0; j<ln; j++) bcol[j] = lidxn[j]*bs+k;
+      err = MatSetValuesLocal(A,lm,brow,ln,bcol,bval,imode);dCHK(err);
+    }
+  } else {
+    err = MatSetValuesBlockedLocal(A,lm,lidxm,ln,lidxn,lv,imode);dCHK(err);
+  }
 
   if (lidxm != lidxms) {err = dFree(lidxm);dCHK(err);}
   if (lidxn != lidxns) {err = dFree(lidxn);dCHK(err);}

src/fs/interface/fsreg.c

@@ -91,6 +91,7 @@ dErr dFSSetFromOptions(dFS fs)
   }
   err = PetscOptionsInt("-dfs_rule_strength","Choose rules that are stronger than necessary","dFSRuleStrength",fs->ruleStrength,&fs->ruleStrength,NULL);dCHK(err);
   err = PetscOptionsList("-dfs_ordering_type","Function Space ordering, usually to reduce bandwidth","dFSBuildSpace",MatOrderingList,fs->orderingtype,fs->orderingtype,256,NULL);dCHK(err);
+  err = PetscOptionsTruth("-dfs_assemble_reduced","Assemble only diagonal part of blocks","",fs->assemblereduced,&fs->assemblereduced,NULL);dCHK(err);
   if (fs->ops->setfromoptions) {
     err = (*fs->ops->setfromoptions)(fs);dCHK(err);
   }