Recmat performance

src/struct_ls/pfmg_setup.c

@@ -151,11 +151,13 @@ hypre_PFMGSetup( void               *pfmg_vdata,
    A_l[0] = hypre_StructMatrixRef(A);
    b_l[0] = hypre_StructVectorRef(b);
    x_l[0] = hypre_StructVectorRef(x);
+   hypre_StructVectorResizeBehavior(x_l[0]) = 1;
 
    tx_l[0] = hypre_StructVectorCreate(comm, grid_l[0]);
    hypre_StructVectorSetNumGhost(tx_l[0], x_num_ghost);
    hypre_StructVectorInitialize(tx_l[0]);
    hypre_StructVectorAssemble(tx_l[0]);
+   hypre_StructVectorResizeBehavior(tx_l[0]) = 1;
 
    //   /* RDF AP Debug */
    //   hypre_StructAssumedPartitionPrint("zAP", hypre_BoxManAssumedPartition(
@@ -230,11 +232,13 @@ hypre_PFMGSetup( void               *pfmg_vdata,
       hypre_StructVectorSetNumGhost(x_l[l + 1], x_num_ghost);
       hypre_StructVectorInitialize(x_l[l + 1]);
       hypre_StructVectorAssemble(x_l[l + 1]);
+      hypre_StructVectorResizeBehavior(x_l[l + 1]) = 1;
 
       tx_l[l + 1] = hypre_StructVectorCreate(comm, grid_l[l + 1]);
       hypre_StructVectorSetNumGhost(tx_l[l + 1], x_num_ghost);
       hypre_StructVectorInitialize(tx_l[l + 1]);
       hypre_StructVectorAssemble(tx_l[l + 1]);
+      hypre_StructVectorResizeBehavior(tx_l[l + 1]) = 1;
 
       HYPRE_ANNOTATE_MGLEVEL_END(l);
    }

src/struct_mv/HYPRE_struct_mv.h

@@ -514,6 +514,17 @@ HYPRE_Int HYPRE_StructVectorCreate(MPI_Comm            comm,
  **/
 HYPRE_Int HYPRE_StructVectorDestroy(HYPRE_StructVector vector);
 
+/**
+ * Set the resize behavior for the struct vector.
+ *    0 = resize/restore on every matvec
+ *    1 = resize only when necessary, retain
+ *        initial data space in save_data_space
+ *    2 = resize only when necessary, retain
+ *        most recent data space in save_data_space
+ **/
+HYPRE_Int HYPRE_StructVectorSetResizeBehavior(HYPRE_StructVector vector,
+                                              HYPRE_Int          resize_behavior);
+
 /* RDF: Need a good user interface for setting the grid. */
 
 /**

src/struct_mv/HYPRE_struct_vector.c

@@ -37,6 +37,18 @@ HYPRE_StructVectorDestroy( HYPRE_StructVector struct_vector )
    return ( hypre_StructVectorDestroy(struct_vector) );
 }
 
+/*--------------------------------------------------------------------------
+ * HYPRE_StructVectorSetResizeBehavior
+ *--------------------------------------------------------------------------*/
+
+HYPRE_Int
+HYPRE_StructVectorSetResizeBehavior(HYPRE_StructVector vector,
+                                    HYPRE_Int          resize_behavior)
+{
+   hypre_StructVectorResizeBehavior(vector) = resize_behavior;
+   return hypre_error_flag;
+}
+
 /*--------------------------------------------------------------------------
  *--------------------------------------------------------------------------*/
 

src/struct_mv/_hypre_struct_mv.h

@@ -1668,6 +1668,11 @@ typedef struct hypre_StructVector_struct
    HYPRE_Int             ref_count;
 
    /* Information needed to Restore() after Reindex() and Resize() */
+   HYPRE_Int             resize_behavior;             /* 0 = resize/restore on every matvec
+                                                         1 = resize when necessary, keep most
+                                                             recent data space in save_data_space
+                                                         2 = resize when necessary, keep initial
+                                                             data space in save_data_space */
    hypre_StructGrid     *save_grid;
    hypre_Index           save_stride;
    HYPRE_Complex        *save_data;                   /* Only needed to support InitializeData() */
@@ -1695,6 +1700,7 @@ typedef struct hypre_StructVector_struct
 #define hypre_StructVectorBGhostNotClear(vector) ((vector) -> bghost_not_clear)
 #define hypre_StructVectorGlobalSize(vector)     ((vector) -> global_size)
 #define hypre_StructVectorRefCount(vector)       ((vector) -> ref_count)
+#define hypre_StructVectorResizeBehavior(vector) ((vector) -> resize_behavior)
 #define hypre_StructVectorSaveGrid(vector)      ((vector) -> save_grid)
 #define hypre_StructVectorSaveStride(vector)    ((vector) -> save_stride)
 #define hypre_StructVectorSaveData(vector)      ((vector) -> save_data)
@@ -1857,6 +1863,7 @@ HYPRE_Int hypre_SubtractBoxArrays ( hypre_BoxArray *box_array1, hypre_BoxArray *
                                     hypre_BoxArray *tmp_box_array );
 HYPRE_Int hypre_UnionBoxes ( hypre_BoxArray *boxes );
 HYPRE_Int hypre_MinUnionBoxes ( hypre_BoxArray *boxes );
+HYPRE_Int hypre_BoxArrayContains( hypre_BoxArray *box_array1, hypre_BoxArray *box_array2 );
 
 /* box_boundary.c */
 HYPRE_Int hypre_BoxBoundaryIntersect ( hypre_Box *box, hypre_StructGrid *grid, HYPRE_Int d,

src/struct_mv/box_algebra.c

@@ -662,3 +662,56 @@ hypre_MinUnionBoxes( hypre_BoxArray *boxes )
 
    return hypre_error_flag;
 }
+
+/*--------------------------------------------------------------------------
+ * hypre_BoxArrayContains
+ *
+ * Check whether box_array1 contains (is a superset of) box_array2
+ * WM: NOTE - this is really checking whether there is a box in box_array1
+ *            that can cover each box in box_array2... so we don't check
+ *            whether a union of multiple boxes from box_array1 will cover
+ *            boxes in box_array2.
+ *--------------------------------------------------------------------------*/
+
+HYPRE_Int
+hypre_BoxArrayContains( hypre_BoxArray *box_array1,
+                        hypre_BoxArray *box_array2 )
+{
+   /* WM: TODO - test this function! */
+   HYPRE_Int i, j, d, covered = 0;
+   HYPRE_Int ndim = hypre_BoxArrayNDim(box_array2);
+   hypre_Box *box1;
+   hypre_Box *box2;
+
+   hypre_ForBoxI(i, box_array2)
+   {
+      box2 = hypre_BoxArrayBox(box_array2, i);
+      covered = 0;
+      hypre_ForBoxI(j, box_array1)
+      {
+         box1 = hypre_BoxArrayBox(box_array1, j);
+         for (d = 0; d < ndim; d++)
+         {
+            if ( hypre_BoxIMinD(box1, d) <= hypre_BoxIMinD(box2, d)
+                 && hypre_BoxIMaxD(box1, d) >= hypre_BoxIMaxD(box2, d) )
+            {
+               covered++;
+            }
+         }
+         if (covered == ndim)
+         {
+            break;
+         }
+         else
+         {
+            covered = 0;
+         }
+      }
+      if (!covered)
+      {
+         return 0;
+      }
+   }
+
+   return 1;
+}

src/struct_mv/protos.h

@@ -49,6 +49,7 @@ HYPRE_Int hypre_SubtractBoxArrays ( hypre_BoxArray *box_array1, hypre_BoxArray *
                                     hypre_BoxArray *tmp_box_array );
 HYPRE_Int hypre_UnionBoxes ( hypre_BoxArray *boxes );
 HYPRE_Int hypre_MinUnionBoxes ( hypre_BoxArray *boxes );
+HYPRE_Int hypre_BoxArrayContains( hypre_BoxArray *box_array1, hypre_BoxArray *box_array2 );
 
 /* box_boundary.c */
 HYPRE_Int hypre_BoxBoundaryIntersect ( hypre_Box *box, hypre_StructGrid *grid, HYPRE_Int d,

src/struct_mv/struct_matvec.c

@@ -14,7 +14,7 @@
 #include "_hypre_struct_mv.h"
 #include "_hypre_struct_mv.hpp"
 
-#define UNROLL_MAXDEPTH 4
+#define UNROLL_MAXDEPTH 9
 
 /*--------------------------------------------------------------------------
  * Matvec data structure
@@ -236,6 +236,7 @@ hypre_StructMatvecCompute( void               *matvec_vdata,
    HYPRE_Int                ndim        = hypre_StructMatrixNDim(A);
 
    hypre_BoxArray          *data_space;
+   hypre_BoxArray          *old_data_space;
    hypre_BoxArrayArray     *compute_box_aa;
    hypre_BoxArray          *compute_box_a;
    hypre_Box               *compute_box;
@@ -256,6 +257,7 @@ hypre_StructMatvecCompute( void               *matvec_vdata,
    hypre_StructStencil     *stencil;
    hypre_Index             *stencil_shape;
    HYPRE_Int                stencil_size;
+   HYPRE_Int                old_data_size;
 
    hypre_StructGrid        *grid;
    HYPRE_Int                ran_nboxes;
@@ -349,7 +351,30 @@ hypre_StructMatvecCompute( void               *matvec_vdata,
    /* This resizes the data for x using the data_space computed during setup */
    data_space = hypre_BoxArrayClone(matvec_data -> data_space);
    hypre_StructVectorReindex(x, grid, dom_stride);
-   hypre_StructVectorResize(x, data_space);
+   if (hypre_StructVectorResizeBehavior(x))
+   {
+      if (!hypre_BoxArrayContains(hypre_StructVectorDataSpace(x), data_space))
+      {
+         if (hypre_StructVectorResizeBehavior(x) > 1)
+         {
+            old_data_space = hypre_StructVectorSaveDataSpace(x);
+            old_data_size = hypre_StructVectorSaveDataSize(x);
+         }
+
+         hypre_StructVectorForget(x);
+         hypre_StructVectorResize(x, data_space);
+
+         if (hypre_StructVectorResizeBehavior(x) > 1)
+         {
+            hypre_StructVectorSaveDataSpace(x) = old_data_space;
+            hypre_StructVectorSaveDataSize(x)  = old_data_size;
+         }
+      }
+   }
+   else
+   {
+      hypre_StructVectorResize(x, data_space);
+   }
 
    stencil       = hypre_StructMatrixStencil(A);
    stencil_shape = hypre_StructStencilShape(stencil);
@@ -600,8 +625,19 @@ hypre_StructMatvecCompute( void               *matvec_vdata,
    }
    else
    {
-      /* Restore the original grid and data layout for x */
-      hypre_StructVectorRestore(x);
+      if (hypre_StructVectorResizeBehavior(x))
+      {
+         /* Just restore the grid (call restore with SaveDataSpace set to NULL) */
+         old_data_space = hypre_StructVectorSaveDataSpace(x);
+         hypre_StructVectorSaveDataSpace(x) = NULL;
+         hypre_StructVectorRestore(x);
+         hypre_StructVectorSaveDataSpace(x) = old_data_space;
+      }
+      else
+      {
+         /* Restore the original grid and data layout for x */
+         hypre_StructVectorRestore(x);
+      }
    }
    hypre_BoxDestroy(compute_box);
 

src/struct_mv/struct_vector.c

@@ -178,6 +178,7 @@ hypre_StructVectorCreate( MPI_Comm          comm,
    hypre_StructVectorDataAlloced(vector)    = 0;
    hypre_StructVectorBGhostNotClear(vector) = 0;
    hypre_StructVectorRefCount(vector)       = 1;
+   hypre_StructVectorResizeBehavior(vector) = 0;
 
    /* set defaults */
    for (i = 0; i < 2 * ndim; i++)
@@ -405,13 +406,15 @@ hypre_StructVectorResize( hypre_StructVector *vector,
       }
 
       /* Free up some things */
+      /* WM: question - won't the if below always be true (since DataAlloced is check above in order to get here)? */
       if (hypre_StructVectorDataAlloced(vector) == 1)
       {
          hypre_TFree(old_data, memory_location);
       }
       hypre_TFree(old_data_indices, HYPRE_MEMORY_HOST);
 
       /* Save old data */
+      /* WM: quesetion - related to the above, isn't old_data always just a null pointer? */
       hypre_StructVectorSaveData(vector)      = old_data;
       hypre_StructVectorSaveDataSpace(vector) = old_data_space;
       hypre_StructVectorSaveDataSize(vector)  = old_data_size;

src/struct_mv/struct_vector.h

@@ -52,6 +52,11 @@ typedef struct hypre_StructVector_struct
    HYPRE_Int             ref_count;
 
    /* Information needed to Restore() after Reindex() and Resize() */
+   HYPRE_Int             resize_behavior;             /* 0 = resize/restore on every matvec
+                                                         1 = resize when necessary, keep most
+                                                             recent data space in save_data_space
+                                                         2 = resize when necessary, keep initial
+                                                             data space in save_data_space */
    hypre_StructGrid     *save_grid;
    hypre_Index           save_stride;
    HYPRE_Complex        *save_data;                   /* Only needed to support InitializeData() */
@@ -79,6 +84,7 @@ typedef struct hypre_StructVector_struct
 #define hypre_StructVectorBGhostNotClear(vector) ((vector) -> bghost_not_clear)
 #define hypre_StructVectorGlobalSize(vector)     ((vector) -> global_size)
 #define hypre_StructVectorRefCount(vector)       ((vector) -> ref_count)
+#define hypre_StructVectorResizeBehavior(vector) ((vector) -> resize_behavior)
 #define hypre_StructVectorSaveGrid(vector)      ((vector) -> save_grid)
 #define hypre_StructVectorSaveStride(vector)    ((vector) -> save_stride)
 #define hypre_StructVectorSaveData(vector)      ((vector) -> save_data)