/*BHEADER**********************************************************************
 * (c) 1997   The Regents of the University of California
 *
 * See the file COPYRIGHT_and_DISCLAIMER for a complete copyright
 * notice, contact person, and disclaimer.
 *
 * $Revision$
 *********************************************************************EHEADER*/
/******************************************************************************
 * Cyclic reduction algorithm (coded as if it were a 1D MG method)
 *
 *****************************************************************************/

#include "headers.h"

#define DEBUG 0

/*--------------------------------------------------------------------------
 * Macros
 *--------------------------------------------------------------------------*/

#define hypre_CycRedSetCIndex(base_index, base_stride, level, cdir, cindex) \
{\
   if (level > 0)\
      hypre_SetIndex(cindex, 0, 0, 0);\
   else\
      hypre_CopyIndex(base_index,  cindex);\
   hypre_IndexD(cindex, cdir) += 0;\
}

#define hypre_CycRedSetFIndex(base_index, base_stride, level, cdir, findex) \
{\
   if (level > 0)\
      hypre_SetIndex(findex, 0, 0, 0);\
   else\
      hypre_CopyIndex(base_index,  findex);\
   hypre_IndexD(findex, cdir) += 1;\
}

#define hypre_CycRedSetStride(base_index, base_stride, level, cdir, stride) \
{\
   if (level > 0)\
      hypre_SetIndex(stride, 1, 1, 1);\
   else\
      hypre_CopyIndex(base_stride, stride);\
   hypre_IndexD(stride, cdir) *= 2;\
}

/*--------------------------------------------------------------------------
 * hypre_CyclicReductionData data structure
 *--------------------------------------------------------------------------*/

typedef struct
{
   MPI_Comm              comm;
                      
   int                   num_levels;
                      
   int                   cdir;         /* coarsening direction */
   hypre_Index           base_index;
   hypre_Index           base_stride;

   hypre_StructGrid    **grid_l;
                    
   hypre_BoxArray       *base_points;
   hypre_BoxArray      **fine_points_l;

   double               *data;
   hypre_StructMatrix  **A_l;
   hypre_StructVector  **x_l;

   hypre_ComputePkg    **down_compute_pkg_l;
   hypre_ComputePkg    **up_compute_pkg_l;

   int                   time_index;
   int                   solve_flops;

} hypre_CyclicReductionData;

/*--------------------------------------------------------------------------
 * hypre_CyclicReductionCreate
 *--------------------------------------------------------------------------*/

void *
hypre_CyclicReductionCreate( MPI_Comm  comm )
{
   hypre_CyclicReductionData *cyc_red_data;

   cyc_red_data = hypre_CTAlloc(hypre_CyclicReductionData, 1);
   
   (cyc_red_data -> comm) = comm;
   (cyc_red_data -> cdir) = 0;
   (cyc_red_data -> time_index)  = hypre_InitializeTiming("CyclicReduction");

   /* set defaults */
   hypre_SetIndex((cyc_red_data -> base_index), 0, 0, 0);
   hypre_SetIndex((cyc_red_data -> base_stride), 1, 1, 1);

   return (void *) cyc_red_data;
}

/*--------------------------------------------------------------------------
 * hypre_CycRedCreateCoarseOp
 *--------------------------------------------------------------------------*/

hypre_StructMatrix *
hypre_CycRedCreateCoarseOp( hypre_StructMatrix *A,
                            hypre_StructGrid   *coarse_grid,
                            int                 cdir        )
{
   hypre_StructMatrix    *Ac;

   hypre_Index           *Ac_stencil_shape;
   hypre_StructStencil   *Ac_stencil;
   int                    Ac_stencil_size;
   int                    Ac_stencil_dim;
   int                    Ac_num_ghost[] = {0, 0, 0, 0, 0, 0};
                       
   int                    i;
   int                    stencil_rank;
 
   Ac_stencil_dim = 1;

   /*-----------------------------------------------
    * Define Ac_stencil
    *-----------------------------------------------*/

   stencil_rank = 0;

   /*-----------------------------------------------
    * non-symmetric case:
    *
    * 3 point fine grid stencil produces 3 point Ac
    *-----------------------------------------------*/

   if (!hypre_StructMatrixSymmetric(A))
   {
      Ac_stencil_size = 3;
      Ac_stencil_shape = hypre_CTAlloc(hypre_Index, Ac_stencil_size);
      for (i = -1; i < 2; i++)
      {
         /* Storage for 3 elements (c,w,e) */
         hypre_SetIndex(Ac_stencil_shape[stencil_rank],i,0,0);
         stencil_rank++;
      }
   }

   /*-----------------------------------------------
    * symmetric case:
    *
    * 3 point fine grid stencil produces 3 point Ac
    *
    * Only store the lower triangular part + diagonal = 2 entries,
    * lower triangular means the lower triangular part on the matrix
    * in the standard lexicalgraphic ordering.
    *-----------------------------------------------*/

   else
   {
      Ac_stencil_size = 2;
      Ac_stencil_shape = hypre_CTAlloc(hypre_Index, Ac_stencil_size);
      for (i = -1; i < 1; i++)
      {

         /* Storage for 2 elements in (c,w) */
         hypre_SetIndex(Ac_stencil_shape[stencil_rank],i,0,0);
         stencil_rank++;
      }
   }

   Ac_stencil = hypre_StructStencilCreate(Ac_stencil_dim, Ac_stencil_size,
                                       Ac_stencil_shape);

   Ac = hypre_StructMatrixCreate(hypre_StructMatrixComm(A),
                              coarse_grid, Ac_stencil);

   hypre_StructStencilDestroy(Ac_stencil);

   /*-----------------------------------------------
    * Coarse operator in symmetric iff fine operator is
    *-----------------------------------------------*/

   hypre_StructMatrixSymmetric(Ac) = hypre_StructMatrixSymmetric(A);

   /*-----------------------------------------------
    * Set number of ghost points
    *-----------------------------------------------*/

   Ac_num_ghost[2*cdir]     = 1;
   if (!hypre_StructMatrixSymmetric(A))
   {
     Ac_num_ghost[2*cdir + 1] = 1;
   }
   hypre_StructMatrixSetNumGhost(Ac, Ac_num_ghost);

   hypre_StructMatrixInitializeShell(Ac);
 
   return Ac;
}

/*--------------------------------------------------------------------------
 * hypre_CycRedSetupCoarseOp
 *--------------------------------------------------------------------------*/

int
hypre_CycRedSetupCoarseOp( hypre_StructMatrix *A,
                           hypre_StructMatrix *Ac,
                           hypre_Index         cindex,
                           hypre_Index         cstride )

{
   hypre_Index             index;

   hypre_StructGrid       *fgrid;
   int                    *fgrid_ids;
   hypre_StructGrid       *cgrid;
   hypre_BoxArray         *cgrid_boxes;
   int                    *cgrid_ids;
   hypre_Box              *cgrid_box;
   hypre_IndexRef          cstart;
   hypre_Index             stridec;
   hypre_Index             fstart;
   hypre_IndexRef          stridef;
   hypre_Index             loop_size;

   int                     fi, ci;
   int                     loopi, loopj, loopk;

   hypre_Box              *A_dbox;
   hypre_Box              *Ac_dbox;

   double                 *a_cc, *a_cw, *a_ce;
   double                 *ac_cc, *ac_cw, *ac_ce;
                    
   int                     iA, iAm1, iAp1;
   int                     iAc;
                         
   int                     xOffsetA; 
                         
   int                     ierr = 0;

   stridef = cstride;
   hypre_SetIndex(stridec, 1, 1, 1);

   fgrid = hypre_StructMatrixGrid(A);
   fgrid_ids = hypre_StructGridIDs(fgrid);

   cgrid = hypre_StructMatrixGrid(Ac);
   cgrid_boxes = hypre_StructGridBoxes(cgrid);
   cgrid_ids = hypre_StructGridIDs(cgrid);

   fi = 0;
   hypre_ForBoxI(ci, cgrid_boxes)
      {
         while (fgrid_ids[fi] != cgrid_ids[ci])
         {
            fi++;
         }

         cgrid_box = hypre_BoxArrayBox(cgrid_boxes, ci);

         cstart = hypre_BoxIMin(cgrid_box);
         hypre_StructMapCoarseToFine(cstart, cindex, cstride, fstart);

         A_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(A), fi);
         Ac_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(Ac), ci);

         /*-----------------------------------------------
          * Extract pointers for 3-point fine grid operator:
          * 
          * a_cc is pointer for center coefficient
          * a_cw is pointer for west coefficient
          * a_ce is pointer for east coefficient
          *-----------------------------------------------*/

         hypre_SetIndex(index,0,0,0);
         a_cc = hypre_StructMatrixExtractPointerByIndex(A, fi, index);

         hypre_SetIndex(index,-1,0,0);
         a_cw = hypre_StructMatrixExtractPointerByIndex(A, fi, index);

         hypre_SetIndex(index,1,0,0);
         a_ce = hypre_StructMatrixExtractPointerByIndex(A, fi, index);

         /*-----------------------------------------------
          * Extract pointers for coarse grid operator - always 3-point:
          *
          * If A is symmetric so is Ac.  We build only the
          * lower triangular part (plus diagonal).
          * 
          * ac_cc is pointer for center coefficient (etc.)
          *-----------------------------------------------*/

         hypre_SetIndex(index,0,0,0);
         ac_cc = hypre_StructMatrixExtractPointerByIndex(Ac, ci, index);

         hypre_SetIndex(index,-1,0,0);
         ac_cw = hypre_StructMatrixExtractPointerByIndex(Ac, ci, index);

         if(!hypre_StructMatrixSymmetric(A))
         {
            hypre_SetIndex(index,1,0,0);
            ac_ce = hypre_StructMatrixExtractPointerByIndex(Ac, ci, index);
         }

         /*-----------------------------------------------
          * Define offsets for fine grid stencil and interpolation
          *
          * In the BoxLoop below I assume iA and iP refer
          * to data associated with the point which we are
          * building the stencil for.  The below offsets
          * are used in refering to data associated with
          * other points. 
          *-----------------------------------------------*/

         hypre_SetIndex(index,1,0,0);
         xOffsetA = hypre_BoxOffsetDistance(A_dbox,index); 

         /*-----------------------------------------------
          * non-symmetric case
          *-----------------------------------------------*/

         if(!hypre_StructMatrixSymmetric(A))
         {
            hypre_BoxGetSize(cgrid_box, loop_size);

            hypre_BoxLoop2Begin(loop_size,
                                A_dbox, fstart, stridef, iA,
                                Ac_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iA,iAc,iAm1,iAp1
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop2For(loopi, loopj, loopk, iA, iAc)
               {
                  iAm1 = iA - xOffsetA;
                  iAp1 = iA + xOffsetA;

                  ac_cw[iAc] = - a_cw[iA] *a_cw[iAm1] / a_cc[iAm1];

                  ac_cc[iAc] = a_cc[iA]
                     - a_cw[iA] * a_ce[iAm1] / a_cc[iAm1]   
                     - a_ce[iA] * a_cw[iAp1] / a_cc[iAp1];   

                  ac_ce[iAc] = - a_ce[iA] *a_ce[iAp1] / a_cc[iAp1];

               }
            hypre_BoxLoop2End(iA, iAc);
         }

         /*-----------------------------------------------
          * symmetric case
          *-----------------------------------------------*/

         else
         {
            hypre_BoxGetSize(cgrid_box, loop_size);

            hypre_BoxLoop2Begin(loop_size,
                                A_dbox, fstart, stridef, iA,
                                Ac_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iA,iAc,iAm1,iAp1
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop2For(loopi, loopj, loopk, iA, iAc)
               {
                  iAm1 = iA - xOffsetA;
                  iAp1 = iA + xOffsetA;

                  ac_cw[iAc] = - a_cw[iA] *a_cw[iAm1] / a_cc[iAm1];

                  ac_cc[iAc] = a_cc[iA]
                     - a_cw[iA] * a_ce[iAm1] / a_cc[iAm1]   
                     - a_ce[iA] * a_cw[iAp1] / a_cc[iAp1];   
               }
            hypre_BoxLoop2End(iA, iAc);
         }

      } /* end ForBoxI */

   hypre_StructMatrixAssemble(Ac);

   /*-----------------------------------------------------------------------
    * Collapse stencil in periodic direction on coarsest grid.
    *-----------------------------------------------------------------------*/

   if (hypre_IndexX(hypre_StructGridPeriodic(cgrid)) == 1)
   {
      hypre_ForBoxI(ci, cgrid_boxes)
      {
         cgrid_box = hypre_BoxArrayBox(cgrid_boxes, ci);

         cstart = hypre_BoxIMin(cgrid_box);

         Ac_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(Ac), ci);

         /*-----------------------------------------------
          * Extract pointers for coarse grid operator - always 3-point:
          *
          * If A is symmetric so is Ac.  We build only the
          * lower triangular part (plus diagonal).
          *
          * ac_cc is pointer for center coefficient (etc.)
          *-----------------------------------------------*/

         hypre_SetIndex(index,0,0,0);
         ac_cc = hypre_StructMatrixExtractPointerByIndex(Ac, ci, index);

         hypre_SetIndex(index,-1,0,0);
         ac_cw = hypre_StructMatrixExtractPointerByIndex(Ac, ci, index);

         if(!hypre_StructMatrixSymmetric(A))
         {
            hypre_SetIndex(index,1,0,0);
            ac_ce = hypre_StructMatrixExtractPointerByIndex(Ac, ci, index);
         }


         /*-----------------------------------------------
          * non-symmetric case
          *-----------------------------------------------*/

         if(!hypre_StructMatrixSymmetric(A))
         {
            hypre_BoxGetSize(cgrid_box, loop_size);

            hypre_BoxLoop1Begin(loop_size,
                                Ac_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iAc
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop1For(loopi, loopj, loopk, iAc)
               {
                  ac_cc[iAc] += (ac_cw[iAc] + ac_ce[iAc]);
                  ac_cw[iAc]  =  0.0;
                  ac_ce[iAc]  =  0.0;
               }
            hypre_BoxLoop1End(iAc);
         }

         /*-----------------------------------------------
          * symmetric case
          *-----------------------------------------------*/

         else
         {
            hypre_BoxGetSize(cgrid_box, loop_size);

            hypre_BoxLoop1Begin(loop_size,
                                Ac_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iAc
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop1For(loopi, loopj, loopk, iAc)
               {
                  ac_cc[iAc] += (2.0  *  ac_cw[iAc]);
                  ac_cw[iAc]  =  0.0;
               }
            hypre_BoxLoop1End(iAc);
         }

      } /* end ForBoxI */

   }

   hypre_StructMatrixAssemble(Ac);

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_CyclicReductionSetup
 *--------------------------------------------------------------------------*/

int
hypre_CyclicReductionSetup( void               *cyc_red_vdata,
                            hypre_StructMatrix *A,
                            hypre_StructVector *b,
                            hypre_StructVector *x             )
{
   hypre_CyclicReductionData *cyc_red_data = cyc_red_vdata;

   MPI_Comm                comm        = (cyc_red_data -> comm);
   int                     cdir        = (cyc_red_data -> cdir);
   hypre_IndexRef          base_index  = (cyc_red_data -> base_index);
   hypre_IndexRef          base_stride = (cyc_red_data -> base_stride);

   int                     num_levels;
   hypre_StructGrid      **grid_l;
   hypre_BoxArray         *base_points;
   hypre_BoxArray        **fine_points_l;
   double                 *data;
   int                     data_size = 0;
   hypre_StructMatrix    **A_l;
   hypre_StructVector    **x_l;
   hypre_ComputePkg      **down_compute_pkg_l;
   hypre_ComputePkg      **up_compute_pkg_l;

   hypre_Index             cindex;
   hypre_Index             findex;
   hypre_Index             stride;

   hypre_BoxArrayArray    *send_boxes;
   hypre_BoxArrayArray    *recv_boxes;
   int                   **send_processes;
   int                   **recv_processes;
   hypre_BoxArrayArray    *indt_boxes;
   hypre_BoxArrayArray    *dept_boxes;
                       
   hypre_StructGrid       *grid;

   hypre_Box              *cbox;
                    
   int                     l;
   int                     flop_divisor;
                         
   int                     x_num_ghost[] = {0, 0, 0, 0, 0, 0};
                         
   int                     ierr = 0;

   /*-----------------------------------------------------
    * Set up coarse grids
    *-----------------------------------------------------*/

   grid = hypre_StructMatrixGrid(A);

   /* Compute a preliminary num_levels value based on the grid */
   cbox = hypre_BoxDuplicate(hypre_StructGridBoundingBox(grid));
   num_levels = hypre_Log2(hypre_BoxSizeD(cbox, cdir)) + 2;

   grid_l    = hypre_TAlloc(hypre_StructGrid *, num_levels);
   hypre_StructGridRef(grid, &grid_l[0]);
   for (l = 0; ; l++)
   {
      /* set cindex and stride */
      hypre_CycRedSetCIndex(base_index, base_stride, l, cdir, cindex);
      hypre_CycRedSetStride(base_index, base_stride, l, cdir, stride);

      /* check to see if we should coarsen */
      if ( hypre_BoxIMinD(cbox, cdir) == hypre_BoxIMaxD(cbox, cdir) )
      {
         /* stop coarsening */
         break;
      }

      /* coarsen cbox */
      hypre_ProjectBox(cbox, cindex, stride);
      hypre_StructMapFineToCoarse(hypre_BoxIMin(cbox), cindex, stride,
                                  hypre_BoxIMin(cbox));
      hypre_StructMapFineToCoarse(hypre_BoxIMax(cbox), cindex, stride,
                                  hypre_BoxIMax(cbox));

      /* coarsen the grid */
      hypre_StructCoarsen(grid_l[l], cindex, stride, 1, &grid_l[l+1]);
   }
   num_levels = l + 1;

   /* free up some things */
   hypre_BoxDestroy(cbox);

   (cyc_red_data -> num_levels)      = num_levels;
   (cyc_red_data -> grid_l)          = grid_l;

   /*-----------------------------------------------------
    * Set up base points
    *-----------------------------------------------------*/

   base_points = hypre_BoxArrayDuplicate(hypre_StructGridBoxes(grid_l[0]));
   hypre_ProjectBoxArray(base_points, base_index, base_stride);

   (cyc_red_data -> base_points) = base_points;

   /*-----------------------------------------------------
    * Set up fine points
    *-----------------------------------------------------*/

   fine_points_l   = hypre_TAlloc(hypre_BoxArray *,  num_levels);

   for (l = 0; l < (num_levels - 1); l++)
   {
      hypre_CycRedSetCIndex(base_index, base_stride, l, cdir, cindex);
      hypre_CycRedSetFIndex(base_index, base_stride, l, cdir, findex);
      hypre_CycRedSetStride(base_index, base_stride, l, cdir, stride);

      fine_points_l[l] =
         hypre_BoxArrayDuplicate(hypre_StructGridBoxes(grid_l[l]));
      hypre_ProjectBoxArray(fine_points_l[l], findex, stride);
   }
  
   fine_points_l[l] =
      hypre_BoxArrayDuplicate(hypre_StructGridBoxes(grid_l[l]));
   if (num_levels == 1)
   {
      hypre_ProjectBoxArray(fine_points_l[l], base_index, base_stride);
   }

   (cyc_red_data -> fine_points_l)   = fine_points_l;

   /*-----------------------------------------------------
    * Set up matrix and vector structures
    *-----------------------------------------------------*/

   A_l  = hypre_TAlloc(hypre_StructMatrix *, num_levels);
   x_l  = hypre_TAlloc(hypre_StructVector *, num_levels);

   A_l[0] = hypre_StructMatrixRef(A);
   x_l[0] = hypre_StructVectorRef(x);

   x_num_ghost[2*cdir]     = 1;
   x_num_ghost[2*cdir + 1] = 1;

   for (l = 0; l < (num_levels - 1); l++)
   {
      A_l[l+1] = hypre_CycRedCreateCoarseOp(A_l[l], grid_l[l+1], cdir);
      data_size += hypre_StructMatrixDataSize(A_l[l+1]);

      x_l[l+1] = hypre_StructVectorCreate(comm, grid_l[l+1]);
      hypre_StructVectorSetNumGhost(x_l[l+1], x_num_ghost);
      hypre_StructVectorInitializeShell(x_l[l+1]);
      data_size += hypre_StructVectorDataSize(x_l[l+1]);
   }

   data = hypre_SharedCTAlloc(double, data_size);

   (cyc_red_data -> data) = data;

   for (l = 0; l < (num_levels - 1); l++)
   {
      hypre_StructMatrixInitializeData(A_l[l+1], data);
      data += hypre_StructMatrixDataSize(A_l[l+1]);
      hypre_StructVectorInitializeData(x_l[l+1], data);
      hypre_StructVectorAssemble(x_l[l+1]);
      data += hypre_StructVectorDataSize(x_l[l+1]);
   }

   (cyc_red_data -> A_l)  = A_l;
   (cyc_red_data -> x_l)  = x_l;

   /*-----------------------------------------------------
    * Set up coarse grid operators
    *-----------------------------------------------------*/

   for (l = 0; l < (num_levels - 1); l++)
   {
      hypre_CycRedSetCIndex(base_index, base_stride, l, cdir, cindex);
      hypre_CycRedSetStride(base_index, base_stride, l, cdir, stride);

      hypre_CycRedSetupCoarseOp(A_l[l], A_l[l+1], cindex, stride);
   }

   /*----------------------------------------------------------
    * Set up compute packages
    *----------------------------------------------------------*/

   down_compute_pkg_l = hypre_TAlloc(hypre_ComputePkg *, (num_levels - 1));
   up_compute_pkg_l   = hypre_TAlloc(hypre_ComputePkg *, (num_levels - 1));

   for (l = 0; l < (num_levels - 1); l++)
   {
      hypre_CycRedSetCIndex(base_index, base_stride, l, cdir, cindex);
      hypre_CycRedSetFIndex(base_index, base_stride, l, cdir, findex);
      hypre_CycRedSetStride(base_index, base_stride, l, cdir, stride);

      hypre_CreateComputeInfo(grid_l[l], hypre_StructMatrixStencil(A_l[l]),
                              &send_boxes, &recv_boxes,
                              &send_processes, &recv_processes,
                              &indt_boxes, &dept_boxes);
 
      /* down-cycle */
      hypre_ProjectBoxArrayArray(send_boxes, findex, stride);
      hypre_ProjectBoxArrayArray(recv_boxes, findex, stride);
      hypre_ProjectBoxArrayArray(indt_boxes, cindex, stride);
      hypre_ProjectBoxArrayArray(dept_boxes, cindex, stride);
      hypre_ComputePkgCreate(send_boxes, recv_boxes,
                             stride, stride,
                             send_processes, recv_processes,
                             indt_boxes, dept_boxes,
                             stride, grid_l[l],
                             hypre_StructVectorDataSpace(x_l[l]), 1,
                             &down_compute_pkg_l[l]);

      hypre_CreateComputeInfo(grid_l[l], hypre_StructMatrixStencil(A_l[l]),
                              &send_boxes, &recv_boxes,
                              &send_processes, &recv_processes,
                              &indt_boxes, &dept_boxes);

      /* up-cycle */
      hypre_ProjectBoxArrayArray(send_boxes, cindex, stride);
      hypre_ProjectBoxArrayArray(recv_boxes, cindex, stride);
      hypre_ProjectBoxArrayArray(indt_boxes, findex, stride);
      hypre_ProjectBoxArrayArray(dept_boxes, findex, stride);
      hypre_ComputePkgCreate(send_boxes, recv_boxes,
                             stride, stride,
                             send_processes, recv_processes,
                             indt_boxes, dept_boxes,
                             stride, grid_l[l],
                             hypre_StructVectorDataSpace(x_l[l]), 1,
                             &up_compute_pkg_l[l]);
   }

   (cyc_red_data -> down_compute_pkg_l) = down_compute_pkg_l;
   (cyc_red_data -> up_compute_pkg_l)   = up_compute_pkg_l;

   /*-----------------------------------------------------
    * Compute solve flops
    *-----------------------------------------------------*/

   flop_divisor = (hypre_IndexX(base_stride) *
                   hypre_IndexY(base_stride) *
                   hypre_IndexZ(base_stride)  );
   (cyc_red_data -> solve_flops) =
      hypre_StructVectorGlobalSize(x_l[0])/2/flop_divisor;
   (cyc_red_data -> solve_flops) +=
      5*hypre_StructVectorGlobalSize(x_l[0])/2/flop_divisor;
   for (l = 1; l < (num_levels - 1); l++)
   {
      (cyc_red_data -> solve_flops) +=
         10*hypre_StructVectorGlobalSize(x_l[l])/2;
   }

   if (num_levels > 1)
   {
      (cyc_red_data -> solve_flops) +=  
          hypre_StructVectorGlobalSize(x_l[l])/2;
   }
   

   /*-----------------------------------------------------
    * Finalize some things
    *-----------------------------------------------------*/

#if DEBUG
   {
      char  filename[255];

      /* debugging stuff */
      for (l = 0; l < num_levels; l++)
      {
         sprintf(filename, "yout_A.%02d", l);
         hypre_StructMatrixPrint(filename, A_l[l], 0);
      }
   }
#endif

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_CyclicReduction
 *
 * The solution vectors on each level are also used to store the
 * right-hand-side data.  We can do this because of the red-black
 * nature of the algorithm and the fact that the method is exact,
 * allowing one to assume initial guesses of zero on all grid levels.
 *--------------------------------------------------------------------------*/

int
hypre_CyclicReduction( void               *cyc_red_vdata,
                       hypre_StructMatrix *A,
                       hypre_StructVector *b,
                       hypre_StructVector *x             )
{
   hypre_CyclicReductionData *cyc_red_data = cyc_red_vdata;

   int                   num_levels      = (cyc_red_data -> num_levels);
   int                   cdir            = (cyc_red_data -> cdir);
   hypre_IndexRef        base_index      = (cyc_red_data -> base_index);
   hypre_IndexRef        base_stride     = (cyc_red_data -> base_stride);
   hypre_BoxArray       *base_points     = (cyc_red_data -> base_points);
   hypre_BoxArray      **fine_points_l   = (cyc_red_data -> fine_points_l);
   hypre_StructMatrix  **A_l             = (cyc_red_data -> A_l);
   hypre_StructVector  **x_l             = (cyc_red_data -> x_l);
   hypre_ComputePkg    **down_compute_pkg_l =
      (cyc_red_data -> down_compute_pkg_l);
   hypre_ComputePkg    **up_compute_pkg_l   =
      (cyc_red_data -> up_compute_pkg_l);
                    
   hypre_StructGrid     *fgrid;
   int                  *fgrid_ids;
   hypre_StructGrid     *cgrid;
   hypre_BoxArray       *cgrid_boxes;
   int                  *cgrid_ids;

   hypre_CommHandle     *comm_handle;
                     
   hypre_BoxArrayArray  *compute_box_aa;
   hypre_BoxArray       *compute_box_a;
   hypre_Box            *compute_box;
                     
   hypre_Box            *A_dbox;
   hypre_Box            *x_dbox;
   hypre_Box            *b_dbox;
   hypre_Box            *xc_dbox;
                     
   double               *Ap, *Awp, *Aep;
   double               *xp, *xwp, *xep;
   double               *bp;
   double               *xcp;
                       
   int                   Ai;
   int                   xi;
   int                   bi;
   int                   xci;
                     
   hypre_Index           cindex;
   hypre_Index           stride;
                       
   hypre_Index           index;
   hypre_Index           loop_size;
   hypre_Index           start;
   hypre_Index           startc;
   hypre_Index           stridec;
                     
   int                   compute_i, fi, ci, j, l;
   int                   loopi, loopj, loopk;
                      
   int                   ierr = 0;

   hypre_BeginTiming(cyc_red_data -> time_index);


   /*--------------------------------------------------
    * Initialize some things
    *--------------------------------------------------*/

   hypre_SetIndex(stridec, 1, 1, 1);

   hypre_StructMatrixDestroy(A_l[0]);
   hypre_StructVectorDestroy(x_l[0]);
   A_l[0] = hypre_StructMatrixRef(A);
   x_l[0] = hypre_StructVectorRef(x);

   /*--------------------------------------------------
    * Copy b into x
    *--------------------------------------------------*/

   compute_box_a = base_points;
   hypre_ForBoxI(fi, compute_box_a)
      {
         compute_box = hypre_BoxArrayBox(compute_box_a, fi);

         x_dbox = hypre_BoxArrayBox(hypre_StructVectorDataSpace(x), fi);
         b_dbox = hypre_BoxArrayBox(hypre_StructVectorDataSpace(b), fi);

         xp = hypre_StructVectorBoxData(x, fi);
         bp = hypre_StructVectorBoxData(b, fi);

         hypre_CopyIndex(hypre_BoxIMin(compute_box), start);
         hypre_BoxGetStrideSize(compute_box, base_stride, loop_size);

         hypre_BoxLoop2Begin(loop_size,
                                x_dbox, start, base_stride, xi,
                                b_dbox, start, base_stride, bi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,xi,bi
#include "hypre_box_smp_forloop.h"
         hypre_BoxLoop2For(loopi, loopj, loopk, xi, bi)
            {
               xp[xi] = bp[bi];
            }
         hypre_BoxLoop2End(xi, bi);
      }

   /*--------------------------------------------------
    * Down cycle:
    *
    * 1) Do an F-relaxation sweep with zero initial guess
    * 2) Compute and inject residual at C-points
    *    - computations are at C-points
    *    - communications are at F-points
    *
    * Notes:
    * - Before these two steps are executed, the
    * fine-grid solution vector contains the right-hand-side.
    * - After these two steps are executed, the fine-grid
    * solution vector contains the right-hand side at
    * C-points and the current solution approximation at
    * F-points.  The coarse-grid solution vector contains
    * the restricted (injected) fine-grid residual.
    * - The coarsest grid solve is built into this loop
    * because it involves the same code as step 1.
    *--------------------------------------------------*/
 
   /* The break out of this loop is just before step 2 below */
   for (l = 0; ; l++)
   {
      /* set cindex and stride */
      hypre_CycRedSetCIndex(base_index, base_stride, l, cdir, cindex);
      hypre_CycRedSetStride(base_index, base_stride, l, cdir, stride);

      /* Step 1 */
      compute_box_a = fine_points_l[l];
      hypre_ForBoxI(fi, compute_box_a)
         {
            compute_box = hypre_BoxArrayBox(compute_box_a, fi);

            A_dbox =
               hypre_BoxArrayBox(hypre_StructMatrixDataSpace(A_l[l]), fi);
            x_dbox =
               hypre_BoxArrayBox(hypre_StructVectorDataSpace(x_l[l]), fi);

            hypre_SetIndex(index, 0, 0, 0);
            Ap = hypre_StructMatrixExtractPointerByIndex(A_l[l], fi, index);
            xp = hypre_StructVectorBoxData(x_l[l], fi);

            hypre_CopyIndex(hypre_BoxIMin(compute_box), start);
            hypre_BoxGetStrideSize(compute_box, stride, loop_size);

            hypre_BoxLoop2Begin(loop_size,
                                A_dbox, start, stride, Ai,
                                x_dbox, start, stride, xi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,Ai,xi
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop2For(loopi, loopj, loopk, Ai, xi)
               {
                  xp[xi] /= Ap[Ai]; 
               }
            hypre_BoxLoop2End(Ai, xi);
         }

      if (l == (num_levels - 1))
         break;

      /* Step 2 */
      fgrid = hypre_StructVectorGrid(x_l[l]);
      fgrid_ids = hypre_StructGridIDs(fgrid);
      cgrid = hypre_StructVectorGrid(x_l[l+1]);
      cgrid_boxes = hypre_StructGridBoxes(cgrid);
      cgrid_ids = hypre_StructGridIDs(cgrid);
      for (compute_i = 0; compute_i < 2; compute_i++)
      {
         switch(compute_i)
         {
            case 0:
            {
               xp = hypre_StructVectorData(x_l[l]);
               hypre_InitializeIndtComputations(down_compute_pkg_l[l], xp,
                                                &comm_handle);
               compute_box_aa =
                  hypre_ComputePkgIndtBoxes(down_compute_pkg_l[l]);
            }
            break;

            case 1:
            {
               hypre_FinalizeIndtComputations(comm_handle);
               compute_box_aa =
                  hypre_ComputePkgDeptBoxes(down_compute_pkg_l[l]);
            }
            break;
         }

         fi = 0;
         hypre_ForBoxArrayI(ci, cgrid_boxes)
            {
               while (fgrid_ids[fi] != cgrid_ids[ci])
               {
                  fi++;
               }

               compute_box_a =
                  hypre_BoxArrayArrayBoxArray(compute_box_aa, fi);

               A_dbox =
                  hypre_BoxArrayBox(hypre_StructMatrixDataSpace(A_l[l]), fi);
               x_dbox =
                  hypre_BoxArrayBox(hypre_StructVectorDataSpace(x_l[l]), fi);
               xc_dbox =
                  hypre_BoxArrayBox(hypre_StructVectorDataSpace(x_l[l+1]), ci);

               xp  = hypre_StructVectorBoxData(x_l[l], fi);
               xcp = hypre_StructVectorBoxData(x_l[l+1], ci);

               hypre_SetIndex(index, -1, 0, 0);
               Awp =
                  hypre_StructMatrixExtractPointerByIndex(A_l[l], fi, index);
               xwp = hypre_StructVectorBoxData(x_l[l], fi) +
                  hypre_BoxOffsetDistance(x_dbox, index);

               hypre_SetIndex(index,  1, 0, 0);
               Aep =
                  hypre_StructMatrixExtractPointerByIndex(A_l[l], fi, index);
               xep = hypre_StructVectorBoxData(x_l[l], fi) +
                  hypre_BoxOffsetDistance(x_dbox, index);

               hypre_ForBoxI(j, compute_box_a)
                  {
                     compute_box = hypre_BoxArrayBox(compute_box_a, j);

                     hypre_CopyIndex(hypre_BoxIMin(compute_box), start);
                     hypre_StructMapFineToCoarse(start, cindex, stride,
                                                 startc);

                     hypre_BoxGetStrideSize(compute_box, stride, loop_size);

                     hypre_BoxLoop3Begin(loop_size,
                                         A_dbox, start, stride, Ai,
                                         x_dbox, start, stride, xi,
                                         xc_dbox, startc, stridec, xci);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,Ai,xi,xci
#include "hypre_box_smp_forloop.h"
                     hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, xci)
                        {
                           xcp[xci] = xp[xi] -
                              Awp[Ai]*xwp[xi] -
                              Aep[Ai]*xep[xi];
                        }
                     hypre_BoxLoop3End(Ai, xi, xci);
                  }
            }
      }
   }

   /*--------------------------------------------------
    * Up cycle:
    *
    * 1) Inject coarse error into fine-grid solution
    *    vector (this is the solution at the C-points)
    * 2) Do an F-relaxation sweep on Ax = 0 and update
    *    solution at F-points
    *    - computations are at F-points
    *    - communications are at C-points
    *--------------------------------------------------*/

   for (l = (num_levels - 2); l >= 0; l--)
   {
      /* set cindex and stride */
      hypre_CycRedSetCIndex(base_index, base_stride, l, cdir, cindex);
      hypre_CycRedSetStride(base_index, base_stride, l, cdir, stride);

      /* Step 1 */
      fgrid = hypre_StructVectorGrid(x_l[l]);
      fgrid_ids = hypre_StructGridIDs(fgrid);
      cgrid = hypre_StructVectorGrid(x_l[l+1]);
      cgrid_boxes = hypre_StructGridBoxes(cgrid);
      cgrid_ids = hypre_StructGridIDs(cgrid);
      fi = 0;
      hypre_ForBoxI(ci, cgrid_boxes)
         {
            while (fgrid_ids[fi] != cgrid_ids[ci])
            {
               fi++;
            }

            compute_box = hypre_BoxArrayBox(cgrid_boxes, ci);

            hypre_CopyIndex(hypre_BoxIMin(compute_box), startc);
            hypre_StructMapCoarseToFine(startc, cindex, stride, start);

            x_dbox =
               hypre_BoxArrayBox(hypre_StructVectorDataSpace(x_l[l]), fi);
            xc_dbox =
               hypre_BoxArrayBox(hypre_StructVectorDataSpace(x_l[l+1]), ci);

            xp  = hypre_StructVectorBoxData(x_l[l], fi);
            xcp = hypre_StructVectorBoxData(x_l[l+1], ci);

            hypre_BoxGetSize(compute_box, loop_size);

            hypre_BoxLoop2Begin(loop_size,
                                x_dbox, start, stride, xi,
                                xc_dbox, startc, stridec, xci);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,xi,xci
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop2For(loopi, loopj, loopk, xi, xci)
               {
                  xp[xi] = xcp[xci];
               }
            hypre_BoxLoop2End(xi, xci);
         }

      /* Step 2 */
      for (compute_i = 0; compute_i < 2; compute_i++)
      {
         switch(compute_i)
         {
            case 0:
            {
               xp = hypre_StructVectorData(x_l[l]);
               hypre_InitializeIndtComputations(up_compute_pkg_l[l], xp,
                                                &comm_handle);
               compute_box_aa =
                  hypre_ComputePkgIndtBoxes(up_compute_pkg_l[l]);
            }
            break;

            case 1:
            {
               hypre_FinalizeIndtComputations(comm_handle);
               compute_box_aa =
                  hypre_ComputePkgDeptBoxes(up_compute_pkg_l[l]);
            }
            break;
         }

         hypre_ForBoxArrayI(fi, compute_box_aa)
            {
               compute_box_a =
                  hypre_BoxArrayArrayBoxArray(compute_box_aa, fi);

               A_dbox =
                  hypre_BoxArrayBox(hypre_StructMatrixDataSpace(A_l[l]), fi);
               x_dbox =
                  hypre_BoxArrayBox(hypre_StructVectorDataSpace(x_l[l]), fi);

               hypre_SetIndex(index, 0, 0, 0);
               Ap = hypre_StructMatrixExtractPointerByIndex(A_l[l], fi, index);
               xp = hypre_StructVectorBoxData(x_l[l], fi);

               hypre_SetIndex(index, -1, 0, 0);
               Awp =
                  hypre_StructMatrixExtractPointerByIndex(A_l[l], fi, index);
               xwp = hypre_StructVectorBoxData(x_l[l], fi) +
                  hypre_BoxOffsetDistance(x_dbox, index);

               hypre_SetIndex(index,  1, 0, 0);
               Aep =
                  hypre_StructMatrixExtractPointerByIndex(A_l[l], fi, index);
               xep = hypre_StructVectorBoxData(x_l[l], fi) +
                  hypre_BoxOffsetDistance(x_dbox, index);

               hypre_ForBoxI(j, compute_box_a)
                  {
                     compute_box = hypre_BoxArrayBox(compute_box_a, j);

                     hypre_CopyIndex(hypre_BoxIMin(compute_box), start);
                     hypre_BoxGetStrideSize(compute_box, stride, loop_size);

                     hypre_BoxLoop2Begin(loop_size,
                                         A_dbox, start, stride, Ai,
                                         x_dbox, start, stride, xi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,Ai,xi
#include "hypre_box_smp_forloop.h"
                     hypre_BoxLoop2For(loopi, loopj, loopk, Ai, xi)
                        {
                           xp[xi] -= (Awp[Ai]*xwp[xi] +
                                      Aep[Ai]*xep[xi]  ) / Ap[Ai];
                        }
                     hypre_BoxLoop2End(Ai, xi);
                  }
            }
      }
   }

   /*-----------------------------------------------------
    * Finalize some things
    *-----------------------------------------------------*/

   hypre_IncFLOPCount(cyc_red_data -> solve_flops);
   hypre_EndTiming(cyc_red_data -> time_index);

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_CyclicReductionSetBase
 *--------------------------------------------------------------------------*/
 
int
hypre_CyclicReductionSetBase( void        *cyc_red_vdata,
                              hypre_Index  base_index,
                              hypre_Index  base_stride )
{
   hypre_CyclicReductionData *cyc_red_data = cyc_red_vdata;
   int                      d;
   int                      ierr = 0;
 
   for (d = 0; d < 3; d++)
   {
      hypre_IndexD((cyc_red_data -> base_index),  d) =
         hypre_IndexD(base_index,  d);
      hypre_IndexD((cyc_red_data -> base_stride), d) =
         hypre_IndexD(base_stride, d);
   }
 
   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_CyclicReductionDestroy
 *--------------------------------------------------------------------------*/

int
hypre_CyclicReductionDestroy( void *cyc_red_vdata )
{
   hypre_CyclicReductionData *cyc_red_data = cyc_red_vdata;

   int l;
   int ierr = 0;

   if (cyc_red_data)
   {
      hypre_BoxArrayDestroy(cyc_red_data -> base_points);
      hypre_StructGridDestroy(cyc_red_data -> grid_l[0]);
      hypre_StructMatrixDestroy(cyc_red_data -> A_l[0]);
      hypre_StructVectorDestroy(cyc_red_data -> x_l[0]);
      for (l = 0; l < ((cyc_red_data -> num_levels) - 1); l++)
      {
         hypre_StructGridDestroy(cyc_red_data -> grid_l[l+1]);
         hypre_BoxArrayDestroy(cyc_red_data -> fine_points_l[l]);
         hypre_StructMatrixDestroy(cyc_red_data -> A_l[l+1]);
         hypre_StructVectorDestroy(cyc_red_data -> x_l[l+1]);
         hypre_ComputePkgDestroy(cyc_red_data -> down_compute_pkg_l[l]);
         hypre_ComputePkgDestroy(cyc_red_data -> up_compute_pkg_l[l]);
      }
      hypre_BoxArrayDestroy(cyc_red_data -> fine_points_l[l]);
      hypre_SharedTFree(cyc_red_data -> data); 
      hypre_TFree(cyc_red_data -> grid_l);
      hypre_TFree(cyc_red_data -> fine_points_l);
      hypre_TFree(cyc_red_data -> A_l);
      hypre_TFree(cyc_red_data -> x_l);
      hypre_TFree(cyc_red_data -> down_compute_pkg_l);
      hypre_TFree(cyc_red_data -> up_compute_pkg_l);

      hypre_FinalizeTiming(cyc_red_data -> time_index);
      hypre_TFree(cyc_red_data);
   }

   return ierr;
}