/*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*/
/******************************************************************************
 *
 *
 *****************************************************************************/

#include "headers.h"
#include "smg.h"

/*--------------------------------------------------------------------------
 * hypre_SMG3CreateRAPOp 
 *    Sets up new coarse grid operator stucture.
 *--------------------------------------------------------------------------*/
 
hypre_StructMatrix *
hypre_SMG3CreateRAPOp( hypre_StructMatrix *R,
                       hypre_StructMatrix *A,
                       hypre_StructMatrix *PT,
                       hypre_StructGrid   *coarse_grid )
{
   hypre_StructMatrix    *RAP;

   hypre_Index           *RAP_stencil_shape;
   hypre_StructStencil   *RAP_stencil;
   int                    RAP_stencil_size;
   int                    RAP_stencil_dim;
   int                    RAP_num_ghost[] = {1, 1, 1, 1, 1, 1};

   hypre_StructStencil   *A_stencil;
   int                    A_stencil_size;

   int                    k, j, i;
   int                    stencil_rank;

   RAP_stencil_dim = 3;
 
   A_stencil = hypre_StructMatrixStencil(A);
   A_stencil_size = hypre_StructStencilSize(A_stencil);
 
   /*-----------------------------------------------------------------------
    * Define RAP_stencil
    *-----------------------------------------------------------------------*/

   stencil_rank = 0;

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

   if (!hypre_StructMatrixSymmetric(A))
   {

      /*--------------------------------------------------------------------
       * 7 or 15 point fine grid stencil produces 15 point RAP
       *--------------------------------------------------------------------*/
      if( A_stencil_size <= 15)
      {
         RAP_stencil_size = 15;
         RAP_stencil_shape = hypre_CTAlloc(hypre_Index, RAP_stencil_size);
         for (k = -1; k < 2; k++)
         {
            for (j = -1; j < 2; j++)
            {
               for (i = -1; i < 2; i++)
               {

                  /*--------------------------------------------------------
                   * Storage for c,w,e,n,s elements in each plane
                   *--------------------------------------------------------*/
                  if( i*j == 0 )
                  {
                     hypre_SetIndex(RAP_stencil_shape[stencil_rank],i,j,k);
                     stencil_rank++;
                  }
               }
            }
         }
      }

      /*--------------------------------------------------------------------
       * 19 or 27 point fine grid stencil produces 27 point RAP
       *--------------------------------------------------------------------*/
      else
      {
         RAP_stencil_size = 27;
         RAP_stencil_shape = hypre_CTAlloc(hypre_Index, RAP_stencil_size);
         for (k = -1; k < 2; k++)
         {
            for (j = -1; j < 2; j++)
            {
               for (i = -1; i < 2; i++)
               {

                  /*--------------------------------------------------------
                   * Storage for 9 elements (c,w,e,n,s,sw,se,nw,se) in
                   * each plane
                   *--------------------------------------------------------*/
                  hypre_SetIndex(RAP_stencil_shape[stencil_rank],i,j,k);
                  stencil_rank++;
               }
            }
         }
      }
   }

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

   else
   {

      /*--------------------------------------------------------------------
       * 7 or 15 point fine grid stencil produces 15 point RAP
       * Only store the lower triangular part + diagonal = 8 entries,
       * lower triangular means the lower triangular part on the matrix
       * in the standard lexicalgraphic ordering.
       *--------------------------------------------------------------------*/
      if( A_stencil_size <= 15)
      {
         RAP_stencil_size = 8;
         RAP_stencil_shape = hypre_CTAlloc(hypre_Index, RAP_stencil_size);
         for (k = -1; k < 1; k++)
         {
            for (j = -1; j < 2; j++)
            {
               for (i = -1; i < 2; i++)
               {

                  /*--------------------------------------------------------
                   * Store  5 elements in lower plane (c,w,e,s,n)  
                   * and 3 elements in same plane (c,w,s)
                   *--------------------------------------------------------*/
                  if( i*j == 0 && i+j+k <= 0)
                  {
                     hypre_SetIndex(RAP_stencil_shape[stencil_rank],i,j,k);
                     stencil_rank++;
                  }
               }
            }
         }
      }

      /*--------------------------------------------------------------------
       * 19 or 27 point fine grid stencil produces 27 point RAP
       * Only store the lower triangular part + diagonal = 14 entries,
       * lower triangular means the lower triangular part on the matrix
       * in the standard lexicalgraphic ordering.
       *--------------------------------------------------------------------*/
      else
      {
         RAP_stencil_size = 14;
         RAP_stencil_shape = hypre_CTAlloc(hypre_Index, RAP_stencil_size);
         for (k = -1; k < 1; k++)
         {
            for (j = -1; j < 2; j++)
            {
               for (i = -1; i < 2; i++)
               {

                  /*--------------------------------------------------------
                   * Store  9 elements in lower plane (c,w,e,s,n,sw,se,nw,ne)  
                   * and 5 elements in same plane (c,w,s,sw,se)
                   *--------------------------------------------------------*/
                  if( k < 0 || (i+j+k <=0 && j < 1) )
                  {
                     hypre_SetIndex(RAP_stencil_shape[stencil_rank],i,j,k);
                     stencil_rank++;
                  }
               }
            }
         }
      }
   }

   RAP_stencil = hypre_StructStencilCreate(RAP_stencil_dim, RAP_stencil_size,
                                           RAP_stencil_shape);
   RAP = hypre_StructMatrixCreate(hypre_StructMatrixComm(A),
                                  coarse_grid, RAP_stencil);

   hypre_StructStencilDestroy(RAP_stencil);

   /*-----------------------------------------------------------------------
    * Coarse operator in symmetric iff fine operator is
    *-----------------------------------------------------------------------*/
   hypre_StructMatrixSymmetric(RAP) = hypre_StructMatrixSymmetric(A);

   /*-----------------------------------------------------------------------
    * Set number of ghost points
    *-----------------------------------------------------------------------*/
   if (hypre_StructMatrixSymmetric(A))
   {
      RAP_num_ghost[1] = 0;
      RAP_num_ghost[3] = 0;
      RAP_num_ghost[5] = 0;
   }
   hypre_StructMatrixSetNumGhost(RAP, RAP_num_ghost);

   return RAP;
}

/*--------------------------------------------------------------------------
 * Routines to build RAP. These routines are fairly general
 *  1) No assumptions about symmetry of A
 *  2) No assumption that R = transpose(P)
 *  3) 7,15,19 or 27-point fine grid A 
 *
 * I am, however, assuming that the c-to-c interpolation is the identity.
 *
 * I've written a two routines - hypre_SMG3BuildRAPSym to build the lower
 * triangular part of RAP (including the diagonal) and
 * hypre_SMG3BuildRAPNoSym to build the upper triangular part of RAP
 * (excluding the diagonal). So using symmetric storage, only the first
 * routine would be called. With full storage both would need to be called.
 *
 *--------------------------------------------------------------------------*/

int
hypre_SMG3BuildRAPSym( hypre_StructMatrix *A,
                       hypre_StructMatrix *PT,
                       hypre_StructMatrix *R,
                       hypre_StructMatrix *RAP,
                       hypre_Index         cindex,
                       hypre_Index         cstride )

{

   hypre_Index           index;

   hypre_StructStencil  *fine_stencil;
   int                   fine_stencil_size;

   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            *PT_dbox;
   hypre_Box            *R_dbox;
   hypre_Box            *RAP_dbox;

   double               *pa, *pb;
   double               *ra, *rb;

   double               *a_cc, *a_cw, *a_ce, *a_cs, *a_cn;
   double               *a_ac, *a_aw, *a_as;
   double               *a_bc, *a_bw, *a_be, *a_bs, *a_bn;
   double               *a_csw, *a_cse, *a_cnw, *a_cne;
   double               *a_asw, *a_ase;
   double               *a_bsw, *a_bse, *a_bnw, *a_bne;

   double               *rap_cc, *rap_cw, *rap_cs;
   double               *rap_bc, *rap_bw, *rap_be, *rap_bs, *rap_bn;
   double               *rap_csw, *rap_cse;
   double               *rap_bsw, *rap_bse, *rap_bnw, *rap_bne;

   int                   iA, iAm1, iAp1;
   int                   iAc;
   int                   iP, iP1;
   int                   iR;
                        
   int                   zOffsetA; 
   int                   xOffsetP; 
   int                   yOffsetP; 
   int                   zOffsetP; 
                        
   int                   ierr = 0;

   fine_stencil = hypre_StructMatrixStencil(A);
   fine_stencil_size = hypre_StructStencilSize(fine_stencil);

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

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

   cgrid = hypre_StructMatrixGrid(RAP);
   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);
         PT_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(PT), fi);
         R_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(R), fi);
         RAP_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(RAP), ci);

         /*-----------------------------------------------------------------
          * Extract pointers for interpolation operator:
          * pa is pointer for weight for f-point above c-point 
          * pb is pointer for weight for f-point below c-point 
          *-----------------------------------------------------------------*/

         hypre_SetIndex(index,0,0,1);
         pa = hypre_StructMatrixExtractPointerByIndex(PT, fi, index);

         hypre_SetIndex(index,0,0,-1);
         pb = hypre_StructMatrixExtractPointerByIndex(PT, fi, index);
 
         /*-----------------------------------------------------------------
          * Extract pointers for restriction operator:
          * ra is pointer for weight for f-point above c-point 
          * rb is pointer for weight for f-point below c-point 
          *-----------------------------------------------------------------*/

         hypre_SetIndex(index,0,0,1);
         ra = hypre_StructMatrixExtractPointerByIndex(R, fi, index);

         hypre_SetIndex(index,0,0,-1);
         rb = hypre_StructMatrixExtractPointerByIndex(R, fi, index);
 
         /*-----------------------------------------------------------------
          * Extract pointers for 7-point fine grid operator:
          * 
          * a_cc is pointer for center coefficient
          * a_cw is pointer for west coefficient in same plane
          * a_ce is pointer for east coefficient in same plane
          * a_cs is pointer for south coefficient in same plane
          * a_cn is pointer for north coefficient in same plane
          * a_ac is pointer for center coefficient in plane above
          * a_bc is pointer for center coefficient in plane below
          *-----------------------------------------------------------------*/

         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);

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

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

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

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

         /*-----------------------------------------------------------------
          * Extract additional pointers for 15-point fine grid operator:
          *
          * a_aw is pointer for west coefficient in plane above
          * a_ae is pointer for east coefficient in plane above
          * a_as is pointer for south coefficient in plane above
          * a_an is pointer for north coefficient in plane above
          * a_bw is pointer for west coefficient in plane below
          * a_be is pointer for east coefficient in plane below
          * a_bs is pointer for south coefficient in plane below
          * a_bn is pointer for north coefficient in plane below
          *-----------------------------------------------------------------*/

         if(fine_stencil_size > 7)
         {
            hypre_SetIndex(index,-1,0,1);
            a_aw = hypre_StructMatrixExtractPointerByIndex(A, fi, index);

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

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

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

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

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

         }
  
         /*-----------------------------------------------------------------
          * Extract additional pointers for 19-point fine grid operator:
          *
          * a_csw is pointer for southwest coefficient in same plane
          * a_cse is pointer for southeast coefficient in same plane
          * a_cnw is pointer for northwest coefficient in same plane
          * a_cne is pointer for northeast coefficient in same plane
          *-----------------------------------------------------------------*/

         if(fine_stencil_size > 15)
         {
            hypre_SetIndex(index,-1,-1,0);
            a_csw = hypre_StructMatrixExtractPointerByIndex(A, fi, index);

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

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

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

         }

         /*-----------------------------------------------------------------
          * Extract additional pointers for 27-point fine grid operator:
          *
          * a_asw is pointer for southwest coefficient in plane above
          * a_ase is pointer for southeast coefficient in plane above
          * a_anw is pointer for northwest coefficient in plane above
          * a_ane is pointer for northeast coefficient in plane above
          * a_bsw is pointer for southwest coefficient in plane below
          * a_bse is pointer for southeast coefficient in plane below
          * a_bnw is pointer for northwest coefficient in plane below
          * a_bne is pointer for northeast coefficient in plane below
          *-----------------------------------------------------------------*/

         if(fine_stencil_size > 19)
         {
            hypre_SetIndex(index,-1,-1,1);
            a_asw = hypre_StructMatrixExtractPointerByIndex(A, fi, index);

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

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

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

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

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

         }

         /*-----------------------------------------------------------------
          * Extract pointers for 15-point coarse grid operator:
          *
          * We build only the lower triangular part (plus diagonal).
          * 
          * rap_cc is pointer for center coefficient (etc.)
          *-----------------------------------------------------------------*/

         hypre_SetIndex(index,0,0,0);
         rap_cc = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,-1,0,0);
         rap_cw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,0,-1,0);
         rap_cs = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,0,0,-1);
         rap_bc = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,-1,0,-1);
         rap_bw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,1,0,-1);
         rap_be = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,0,-1,-1);
         rap_bs = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,0,1,-1);
         rap_bn = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         /*-----------------------------------------------------------------
          * Extract additional pointers for 27-point coarse grid operator:
          *
          * A 27-point coarse grid operator is produced when the fine grid 
          * stencil is 19 or 27 point.
          *
          * We build only the lower triangular part.
          *
          * rap_csw is pointer for southwest coefficient in same plane (etc.)
          *-----------------------------------------------------------------*/

         if(fine_stencil_size > 15)
         {
            hypre_SetIndex(index,-1,-1,0);
            rap_csw =
               hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,1,-1,0);
            rap_cse =
               hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,-1,-1,-1);
            rap_bsw =
               hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,1,-1,-1);
            rap_bse =
               hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,-1,1,-1);
            rap_bnw =
               hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,1,1,-1);
            rap_bne =
               hypre_StructMatrixExtractPointerByIndex(RAP, 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,0,0,1);
         zOffsetA = hypre_BoxOffsetDistance(A_dbox,index); 
         zOffsetP = hypre_BoxOffsetDistance(PT_dbox,index); 
         hypre_SetIndex(index,0,1,0);
         yOffsetP = hypre_BoxOffsetDistance(PT_dbox,index); 
         hypre_SetIndex(index,1,0,0);
         xOffsetP = hypre_BoxOffsetDistance(PT_dbox,index); 

         /*--------------------------------------------------------------------
          * Switch statement to direct control to apropriate BoxLoop depending
          * on stencil size. Default is full 27-point.
          *-----------------------------------------------------------------*/

         switch (fine_stencil_size)
         {

            /*--------------------------------------------------------------
             * Loop for symmetric 7-point fine grid operator; produces a
             * symmetric 15-point coarse grid operator. We calculate only the
             * lower triangular stencil entries: (below-south, below-west,
             * below-center, below-east, below-north, center-south,
             * center-west, and center-center).
             *--------------------------------------------------------------*/

            case 7:

            hypre_BoxGetSize(cgrid_box, loop_size);
            hypre_BoxLoop4Begin(loop_size,
                                PT_dbox,  cstart, stridec, iP,
                                R_dbox,   cstart, stridec, iR,
                                A_dbox,   fstart, stridef, iA,
                                RAP_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iP,iR,iA,iAc,iAm1,iAp1,iP1
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop4For(loopi, loopj, loopk, iP, iR, iA, iAc)
               {
                  iAm1 = iA - zOffsetA;
                  iAp1 = iA + zOffsetA;

                  iP1 = iP - zOffsetP - yOffsetP;
                  rap_bs[iAc] = rb[iR] * a_cs[iAm1] * pa[iP1];

                  iP1 = iP - zOffsetP - xOffsetP;
                  rap_bw[iAc] = rb[iR] * a_cw[iAm1] * pa[iP1];
 
                  iP1 = iP - zOffsetP; 
                  rap_bc[iAc] =          a_bc[iA]   * pa[iP1]
                     +          rb[iR] * a_cc[iAm1] * pa[iP1]
                     +          rb[iR] * a_bc[iAm1];
 
                  iP1 = iP - zOffsetP + xOffsetP;
                  rap_be[iAc] = rb[iR] * a_ce[iAm1] * pa[iP1];
 
                  iP1 = iP - zOffsetP + yOffsetP;
                  rap_bn[iAc] = rb[iR] * a_cn[iAm1] * pa[iP1];
 
                  iP1 = iP - yOffsetP;
                  rap_cs[iAc] =          a_cs[iA]
                     +          rb[iR] * a_cs[iAm1] * pb[iP1]
                     +          ra[iR] * a_cs[iAp1] * pa[iP1];
 
                  iP1 = iP - xOffsetP;
                  rap_cw[iAc] =          a_cw[iA]
                     +          rb[iR] * a_cw[iAm1] * pb[iP1]
                     +          ra[iR] * a_cw[iAp1] * pa[iP1];
 
                  rap_cc[iAc] =          a_cc[iA]
                     +          rb[iR] * a_cc[iAm1] * pb[iP]
                     +          ra[iR] * a_cc[iAp1] * pa[iP]
                     +          rb[iR] * a_ac[iAm1]
                     +          ra[iR] * a_bc[iAp1]
                     +                   a_bc[iA]   * pb[iP]
                     +                   a_ac[iA]   * pa[iP];

               }
            hypre_BoxLoop4End(iP, iR, iA, iAc);

            break;

            /*--------------------------------------------------------------
             * Loop for symmetric 15-point fine grid operator; produces a
             * symmetric 15-point coarse grid operator. We calculate only the
             * lower triangular stencil entries: (below-south, below-west,
             * below-center, below-east, below-north, center-south,
             * center-west, and center-center).
             *--------------------------------------------------------------*/

            case 15:

            hypre_BoxGetSize(cgrid_box, loop_size);
            hypre_BoxLoop4Begin(loop_size,
                                PT_dbox,  cstart, stridec, iP,
                                R_dbox,   cstart, stridec, iR,
                                A_dbox,   fstart, stridef, iA,
                                RAP_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iP,iR,iA,iAc,iAm1,iAp1,iP1
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop4For(loopi, loopj, loopk, iP, iR, iA, iAc)
               {
                  iAm1 = iA - zOffsetA;
                  iAp1 = iA + zOffsetA;

                  iP1 = iP - zOffsetP - yOffsetP;
                  rap_bs[iAc] = rb[iR] * a_cs[iAm1] * pa[iP1]
                     +          rb[iR] * a_bs[iAm1]
                     +                   a_bs[iA]   * pa[iP1];

                  iP1 = iP - zOffsetP - xOffsetP;
                  rap_bw[iAc] = rb[iR] * a_cw[iAm1] * pa[iP1]
                     +          rb[iR] * a_bw[iAm1]
                     +                   a_bw[iA]   * pa[iP1];
 
                  iP1 = iP - zOffsetP; 
                  rap_bc[iAc] =          a_bc[iA]   * pa[iP1]
                     +          rb[iR] * a_cc[iAm1] * pa[iP1]
                     +          rb[iR] * a_bc[iAm1];
 
                  iP1 = iP - zOffsetP + xOffsetP;
                  rap_be[iAc] = rb[iR] * a_ce[iAm1] * pa[iP1]
                     +          rb[iR] * a_be[iAm1]
                     +                   a_be[iA]   * pa[iP1];
 
                  iP1 = iP - zOffsetP + yOffsetP;
                  rap_bn[iAc] = rb[iR] * a_cn[iAm1] * pa[iP1]
                     +          rb[iR] * a_bn[iAm1]
                     +                   a_bn[iA]   * pa[iP1];
 
                  iP1 = iP - yOffsetP;
                  rap_cs[iAc] =          a_cs[iA]
                     +          rb[iR] * a_cs[iAm1] * pb[iP1]
                     +          ra[iR] * a_cs[iAp1] * pa[iP1]
                     +                   a_bs[iA]   * pb[iP1]
                     +                   a_as[iA]   * pa[iP1]
                     +          rb[iR] * a_as[iAm1]
                     +          ra[iR] * a_bs[iAp1];
 
                  iP1 = iP - xOffsetP;
                  rap_cw[iAc] =          a_cw[iA]
                     +          rb[iR] * a_cw[iAm1] * pb[iP1]
                     +          ra[iR] * a_cw[iAp1] * pa[iP1]
                     +                   a_bw[iA]   * pb[iP1]
                     +                   a_aw[iA]   * pa[iP1]
                     +          rb[iR] * a_aw[iAm1]
                     +          ra[iR] * a_bw[iAp1];
 
                  rap_cc[iAc] =          a_cc[iA]
                     +          rb[iR] * a_cc[iAm1] * pb[iP]
                     +          ra[iR] * a_cc[iAp1] * pa[iP]
                     +          rb[iR] * a_ac[iAm1]
                     +          ra[iR] * a_bc[iAp1]
                     +                   a_bc[iA]   * pb[iP]
                     +                   a_ac[iA]   * pa[iP];

               }
            hypre_BoxLoop4End(iP, iR, iA, iAc);

            break;

            /*--------------------------------------------------------------
             * Loop for symmetric 19-point fine grid operator; produces a
             * symmetric 27-point coarse grid operator. We calculate only the
             * lower triangular stencil entries: (below-southwest, below-south,
             * below-southeast, below-west, below-center, below-east,
             * below-northwest, below-north, below-northeast, center-southwest,
             * center-south, center-southeast, center-west, and center-center).
             *--------------------------------------------------------------*/

            case 19:

            hypre_BoxGetSize(cgrid_box, loop_size);
            hypre_BoxLoop4Begin(loop_size,
                                PT_dbox,  cstart, stridec, iP,
                                R_dbox,   cstart, stridec, iR,
                                A_dbox,   fstart, stridef, iA,
                                RAP_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iP,iR,iA,iAc,iAm1,iAp1,iP1
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop4For(loopi, loopj, loopk, iP, iR, iA, iAc)
               {
                  iAm1 = iA - zOffsetA;
                  iAp1 = iA + zOffsetA;

                  iP1 = iP - zOffsetP - yOffsetP - xOffsetP;
                  rap_bsw[iAc] = rb[iR] * a_csw[iAm1] * pa[iP1];

                  iP1 = iP - zOffsetP - yOffsetP;
                  rap_bs[iAc] = rb[iR] * a_cs[iAm1] * pa[iP1]
                     +          rb[iR] * a_bs[iAm1]
                     +                   a_bs[iA]   * pa[iP1];

                  iP1 = iP - zOffsetP - yOffsetP + xOffsetP;
                  rap_bse[iAc] = rb[iR] * a_cse[iAm1] * pa[iP1];

                  iP1 = iP - zOffsetP - xOffsetP;
                  rap_bw[iAc] = rb[iR] * a_cw[iAm1] * pa[iP1]
                     +          rb[iR] * a_bw[iAm1]
                     +                   a_bw[iA]   * pa[iP1];
 
                  iP1 = iP - zOffsetP; 
                  rap_bc[iAc] =          a_bc[iA] * pa[iP1]
                     +          rb[iR] * a_cc[iAm1] * pa[iP1]
                     +          rb[iR] * a_bc[iAm1];
 
                  iP1 = iP - zOffsetP + xOffsetP;
                  rap_be[iAc] = rb[iR] * a_ce[iAm1] * pa[iP1]
                     +          rb[iR] * a_be[iAm1]
                     +                   a_be[iA]   * pa[iP1];
 
                  iP1 = iP - zOffsetP + yOffsetP - xOffsetP;
                  rap_bnw[iAc] = rb[iR] * a_cnw[iAm1] * pa[iP1];

                  iP1 = iP - zOffsetP + yOffsetP;
                  rap_bn[iAc] = rb[iR] * a_cn[iAm1] * pa[iP1]
                     +          rb[iR] * a_bn[iAm1]
                     +                   a_bn[iA]   * pa[iP1];
 
                  iP1 = iP - zOffsetP + yOffsetP + xOffsetP;
                  rap_bne[iAc] = rb[iR] * a_cne[iAm1] * pa[iP1];

                  iP1 = iP - yOffsetP - xOffsetP;
                  rap_csw[iAc] =         a_csw[iA]
                     +          rb[iR] * a_csw[iAm1] * pb[iP1]
                     +          ra[iR] * a_csw[iAp1] * pa[iP1];

                  iP1 = iP - yOffsetP;
                  rap_cs[iAc] =          a_cs[iA]
                     +          rb[iR] * a_cs[iAm1] * pb[iP1]
                     +          ra[iR] * a_cs[iAp1] * pa[iP1]
                     +                   a_bs[iA]   * pb[iP1]
                     +                   a_as[iA]   * pa[iP1]
                     +          rb[iR] * a_as[iAm1]
                     +          ra[iR] * a_bs[iAp1];
 
                  iP1 = iP - yOffsetP + xOffsetP;
                  rap_cse[iAc] =          a_cse[iA]
                     +          rb[iR] * a_cse[iAm1] * pb[iP1]
                     +          ra[iR] * a_cse[iAp1] * pa[iP1];

                  iP1 = iP - xOffsetP;
                  rap_cw[iAc] =          a_cw[iA]
                     +          rb[iR] * a_cw[iAm1] * pb[iP1]
                     +          ra[iR] * a_cw[iAp1] * pa[iP1]
                     +                   a_bw[iA]   * pb[iP1]
                     +                   a_aw[iA]   * pa[iP1]
                     +          rb[iR] * a_aw[iAm1]
                     +          ra[iR] * a_bw[iAp1];
 
                  rap_cc[iAc] =          a_cc[iA]
                     +          rb[iR] * a_cc[iAm1] * pb[iP]
                     +          ra[iR] * a_cc[iAp1] * pa[iP]
                     +          rb[iR] * a_ac[iAm1]
                     +          ra[iR] * a_bc[iAp1]
                     +                   a_bc[iA]   * pb[iP]
                     +                   a_ac[iA]   * pa[iP];

               }
            hypre_BoxLoop4End(iP, iR, iA, iAc);

            break;

            /*--------------------------------------------------------------
             * Loop for symmetric 27-point fine grid operator; produces a
             * symmetric 27-point coarse grid operator. We calculate only the
             * lower triangular stencil entries: (below-southwest, below-south,
             * below-southeast, below-west, below-center, below-east,
             * below-northwest, below-north, below-northeast, center-southwest,
             * center-south, center-southeast, center-west, and center-center).
             *--------------------------------------------------------------*/

            default:

            hypre_BoxGetSize(cgrid_box, loop_size);
            hypre_BoxLoop4Begin(loop_size,
                                PT_dbox,  cstart, stridec, iP,
                                R_dbox,   cstart, stridec, iR,
                                A_dbox,   fstart, stridef, iA,
                                RAP_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iP,iR,iA,iAc,iAm1,iAp1,iP1
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop4For(loopi, loopj, loopk, iP, iR, iA, iAc)
               {
                  iAm1 = iA - zOffsetA;
                  iAp1 = iA + zOffsetA;

                  iP1 = iP - zOffsetP - yOffsetP - xOffsetP;
                  rap_bsw[iAc] = rb[iR] * a_csw[iAm1] * pa[iP1]
                     +           rb[iR] * a_bsw[iAm1]
                     +                    a_bsw[iA]   * pa[iP1];

                  iP1 = iP - zOffsetP - yOffsetP;
                  rap_bs[iAc] = rb[iR] * a_cs[iAm1] * pa[iP1]
                     +          rb[iR] * a_bs[iAm1]
                     +                   a_bs[iA]   * pa[iP1];

                  iP1 = iP - zOffsetP - yOffsetP + xOffsetP;
                  rap_bse[iAc] = rb[iR] * a_cse[iAm1] * pa[iP1]
                     +           rb[iR] * a_bse[iAm1]
                     +                    a_bse[iA]   * pa[iP1];

                  iP1 = iP - zOffsetP - xOffsetP;
                  rap_bw[iAc] = rb[iR] * a_cw[iAm1] * pa[iP1]
                     +          rb[iR] * a_bw[iAm1]
                     +                   a_bw[iA]   * pa[iP1];
 
                  iP1 = iP - zOffsetP; 
                  rap_bc[iAc] =          a_bc[iA] * pa[iP1]
                     +          rb[iR] * a_cc[iAm1] * pa[iP1]
                     +          rb[iR] * a_bc[iAm1];
 
                  iP1 = iP - zOffsetP + xOffsetP;
                  rap_be[iAc] = rb[iR] * a_ce[iAm1] * pa[iP1]
                     +          rb[iR] * a_be[iAm1]
                     +                   a_be[iA]   * pa[iP1];
 
                  iP1 = iP - zOffsetP + yOffsetP - xOffsetP;
                  rap_bnw[iAc] = rb[iR] * a_cnw[iAm1] * pa[iP1]
                     +           rb[iR] * a_bnw[iAm1]
                     +                    a_bnw[iA]   * pa[iP1];

                  iP1 = iP - zOffsetP + yOffsetP;
                  rap_bn[iAc] = rb[iR] * a_cn[iAm1] * pa[iP1]
                     +          rb[iR] * a_bn[iAm1]
                     +                   a_bn[iA]   * pa[iP1];
 
                  iP1 = iP - zOffsetP + yOffsetP + xOffsetP;
                  rap_bne[iAc] = rb[iR] * a_cne[iAm1] * pa[iP1]
                     +           rb[iR] * a_bne[iAm1]
                     +                    a_bne[iA]   * pa[iP1];

                  iP1 = iP - yOffsetP - xOffsetP;
                  rap_csw[iAc] =          a_csw[iA]
                     +          rb[iR] * a_csw[iAm1] * pb[iP1]
                     +          ra[iR] * a_csw[iAp1] * pa[iP1]
                     +                   a_bsw[iA]   * pb[iP1]
                     +                   a_asw[iA]   * pa[iP1]
                     +          rb[iR] * a_asw[iAm1]
                     +          ra[iR] * a_bsw[iAp1];

                  iP1 = iP - yOffsetP;
                  rap_cs[iAc] =          a_cs[iA]
                     +          rb[iR] * a_cs[iAm1] * pb[iP1]
                     +          ra[iR] * a_cs[iAp1] * pa[iP1]
                     +                   a_bs[iA]   * pb[iP1]
                     +                   a_as[iA]   * pa[iP1]
                     +          rb[iR] * a_as[iAm1]
                     +          ra[iR] * a_bs[iAp1];
 
                  iP1 = iP - yOffsetP + xOffsetP;
                  rap_cse[iAc] =          a_cse[iA]
                     +          rb[iR] * a_cse[iAm1] * pb[iP1]
                     +          ra[iR] * a_cse[iAp1] * pa[iP1]
                     +                   a_bse[iA]   * pb[iP1]
                     +                   a_ase[iA]   * pa[iP1]
                     +          rb[iR] * a_ase[iAm1]
                     +          ra[iR] * a_bse[iAp1];

                  iP1 = iP - xOffsetP;
                  rap_cw[iAc] =          a_cw[iA]
                     +          rb[iR] * a_cw[iAm1] * pb[iP1]
                     +          ra[iR] * a_cw[iAp1] * pa[iP1]
                     +                   a_bw[iA]   * pb[iP1]
                     +                   a_aw[iA]   * pa[iP1]
                     +          rb[iR] * a_aw[iAm1]
                     +          ra[iR] * a_bw[iAp1];
 
                  rap_cc[iAc] =          a_cc[iA]
                     +          rb[iR] * a_cc[iAm1] * pb[iP]
                     +          ra[iR] * a_cc[iAp1] * pa[iP]
                     +          rb[iR] * a_ac[iAm1]
                     +          ra[iR] * a_bc[iAp1]
                     +                   a_bc[iA]   * pb[iP]
                     +                   a_ac[iA]   * pa[iP];

               }
            hypre_BoxLoop4End(iP, iR, iA, iAc);

            break;

         } /* end switch statement */

      } /* end ForBoxI */

   return ierr;
}

/*--------------------------------------------------------------------------
 *--------------------------------------------------------------------------*/

int
hypre_SMG3BuildRAPNoSym( hypre_StructMatrix *A,
                         hypre_StructMatrix *PT,
                         hypre_StructMatrix *R,
                         hypre_StructMatrix *RAP,
                         hypre_Index         cindex,
                         hypre_Index         cstride )

{

   hypre_Index           index;

   hypre_StructStencil  *fine_stencil;
   int                   fine_stencil_size;

   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            *PT_dbox;
   hypre_Box            *R_dbox;
   hypre_Box            *RAP_dbox;

   double               *pa, *pb;
   double               *ra, *rb;

   double               *a_cc, *a_cw, *a_ce, *a_cs, *a_cn;
   double               *a_ac, *a_aw, *a_ae, *a_as, *a_an;
   double               *a_be, *a_bn;
   double               *a_csw, *a_cse, *a_cnw, *a_cne;
   double               *a_asw, *a_ase, *a_anw, *a_ane;
   double               *a_bnw, *a_bne;

   double               *rap_ce, *rap_cn;
   double               *rap_ac, *rap_aw, *rap_ae, *rap_as, *rap_an;
   double               *rap_cnw, *rap_cne;
   double               *rap_asw, *rap_ase, *rap_anw, *rap_ane;

   int                  iA, iAm1, iAp1;
   int                  iAc;
   int                  iP, iP1;
   int                  iR;

   int                  zOffsetA;
   int                  xOffsetP;
   int                  yOffsetP;
   int                  zOffsetP;

   int                  ierr = 0;

   fine_stencil = hypre_StructMatrixStencil(A);
   fine_stencil_size = hypre_StructStencilSize(fine_stencil);

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

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

   cgrid = hypre_StructMatrixGrid(RAP);
   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);
         PT_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(PT), fi);
         R_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(R), fi);
         RAP_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(RAP), ci);

         /*-----------------------------------------------------------------
          * Extract pointers for interpolation operator:
          * pa is pointer for weight for f-point above c-point 
          * pb is pointer for weight for f-point below c-point 
          *-----------------------------------------------------------------*/

         hypre_SetIndex(index,0,0,1);
         pa = hypre_StructMatrixExtractPointerByIndex(PT, fi, index);

         hypre_SetIndex(index,0,0,-1);
         pb = hypre_StructMatrixExtractPointerByIndex(PT, fi, index);

 
         /*-----------------------------------------------------------------
          * Extract pointers for restriction operator:
          * ra is pointer for weight for f-point above c-point 
          * rb is pointer for weight for f-point below c-point 
          *-----------------------------------------------------------------*/

         hypre_SetIndex(index,0,0,1);
         ra = hypre_StructMatrixExtractPointerByIndex(R, fi, index);

         hypre_SetIndex(index,0,0,-1);
         rb = hypre_StructMatrixExtractPointerByIndex(R, fi, index);

 
         /*-----------------------------------------------------------------
          * Extract pointers for 7-point fine grid operator:
          * 
          * a_cc is pointer for center coefficient
          * a_cw is pointer for west coefficient in same plane
          * a_ce is pointer for east coefficient in same plane
          * a_cs is pointer for south coefficient in same plane
          * a_cn is pointer for north coefficient in same plane
          * a_ac is pointer for center coefficient in plane above
          * a_bc is pointer for center coefficient in plane below
          *-----------------------------------------------------------------*/

         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);

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

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

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

         /*-----------------------------------------------------------------
          * Extract additional pointers for 15-point fine grid operator:
          *
          * a_aw is pointer for west coefficient in plane above
          * a_ae is pointer for east coefficient in plane above
          * a_as is pointer for south coefficient in plane above
          * a_an is pointer for north coefficient in plane above
          * a_bw is pointer for west coefficient in plane below
          * a_be is pointer for east coefficient in plane below
          * a_bs is pointer for south coefficient in plane below
          * a_bn is pointer for north coefficient in plane below
          *-----------------------------------------------------------------*/

         if(fine_stencil_size > 7)
         {
            hypre_SetIndex(index,-1,0,1);
            a_aw = hypre_StructMatrixExtractPointerByIndex(A, fi, index);

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

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

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

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

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

         }
  
         /*-----------------------------------------------------------------
          * Extract additional pointers for 19-point fine grid operator:
          *
          * a_csw is pointer for southwest coefficient in same plane
          * a_cse is pointer for southeast coefficient in same plane
          * a_cnw is pointer for northwest coefficient in same plane
          * a_cne is pointer for northeast coefficient in same plane
          *-----------------------------------------------------------------*/

         if(fine_stencil_size > 15)
         {
            hypre_SetIndex(index,-1,-1,0);
            a_csw = hypre_StructMatrixExtractPointerByIndex(A, fi, index);

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

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

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

         }

         /*-----------------------------------------------------------------
          * Extract additional pointers for 27-point fine grid operator:
          *
          * a_asw is pointer for southwest coefficient in plane above
          * a_ase is pointer for southeast coefficient in plane above
          * a_anw is pointer for northwest coefficient in plane above
          * a_ane is pointer for northeast coefficient in plane above
          * a_bsw is pointer for southwest coefficient in plane below
          * a_bse is pointer for southeast coefficient in plane below
          * a_bnw is pointer for northwest coefficient in plane below
          * a_bne is pointer for northeast coefficient in plane below
          *-----------------------------------------------------------------*/

         if(fine_stencil_size > 19)
         {
            hypre_SetIndex(index,-1,-1,1);
            a_asw = hypre_StructMatrixExtractPointerByIndex(A, fi, index);

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

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

            hypre_SetIndex(index,1,1,1);
            a_ane = hypre_StructMatrixExtractPointerByIndex(A, fi, index);

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

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

         }

         /*-----------------------------------------------------------------
          * Extract pointers for 15-point coarse grid operator:
          *
          * We build only the upper triangular part (excluding diagonal).
          * 
          * rap_ce is pointer for east coefficient in same plane (etc.)
          *-----------------------------------------------------------------*/

         hypre_SetIndex(index,1,0,0);
         rap_ce = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,0,1,0);
         rap_cn = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,0,0,1);
         rap_ac = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,-1,0,1);
         rap_aw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,1,0,1);
         rap_ae = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,0,-1,1);
         rap_as = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         hypre_SetIndex(index,0,1,1);
         rap_an = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

         /*-----------------------------------------------------------------
          * Extract additional pointers for 27-point coarse grid operator:
          *
          * A 27-point coarse grid operator is produced when the fine grid 
          * stencil is 19 or 27 point.
          *
          * We build only the upper triangular part.
          *
          * rap_cnw is pointer for northwest coefficient in same plane (etc.)
          *-----------------------------------------------------------------*/

         if(fine_stencil_size > 15)
         {
            hypre_SetIndex(index,-1,1,0);
            rap_cnw =
               hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,1,1,0);
            rap_cne =
               hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,-1,-1,1);
            rap_asw =
               hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,1,-1,1);
            rap_ase =
               hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,-1,1,1);
            rap_anw =
               hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,1,1,1);
            rap_ane =
               hypre_StructMatrixExtractPointerByIndex(RAP, 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,0,0,1);
         zOffsetA = hypre_BoxOffsetDistance(A_dbox,index); 
         zOffsetP = hypre_BoxOffsetDistance(PT_dbox,index); 
         hypre_SetIndex(index,0,1,0);
         yOffsetP = hypre_BoxOffsetDistance(PT_dbox,index); 
         hypre_SetIndex(index,1,0,0);
         xOffsetP = hypre_BoxOffsetDistance(PT_dbox,index); 

         /*-----------------------------------------------------------------
          * Switch statement to direct control to apropriate BoxLoop depending
          * on stencil size. Default is full 27-point.
          *-----------------------------------------------------------------*/

         switch (fine_stencil_size)
         {

            /*--------------------------------------------------------------
             * Loop for 7-point fine grid operator; produces upper triangular
             * part of 15-point coarse grid operator. stencil entries:
             * (above-north, above-east, above-center, above-west,
             * above-south, center-north, and center-east).
             *--------------------------------------------------------------*/

            case 7:

            hypre_BoxGetSize(cgrid_box, loop_size);
            hypre_BoxLoop4Begin(loop_size,
                                PT_dbox,  cstart, stridec, iP,
                                R_dbox,   cstart, stridec, iR,
                                A_dbox,   fstart, stridef, iA,
                                RAP_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iP,iR,iA,iAc,iAm1,iAp1,iP1
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop4For(loopi, loopj, loopk, iP, iR, iA, iAc)
               {
                  iAm1 = iA - zOffsetA;
                  iAp1 = iA + zOffsetA;

                  iP1 = iP + zOffsetP + yOffsetP;
                  rap_an[iAc] = ra[iR] * a_cn[iAp1] * pb[iP1];

                  iP1 = iP + zOffsetP + xOffsetP;
                  rap_ae[iAc] = ra[iR] * a_ce[iAp1] * pb[iP1];
 
                  iP1 = iP + zOffsetP; 
                  rap_ac[iAc] =          a_ac[iA]   * pb[iP1]
                     +          ra[iR] * a_cc[iAp1] * pb[iP1]
                     +          ra[iR] * a_ac[iAp1];
 
                  iP1 = iP + zOffsetP - xOffsetP;
                  rap_aw[iAc] = ra[iR] * a_cw[iAp1] * pb[iP1];
 
                  iP1 = iP + zOffsetP - yOffsetP;
                  rap_as[iAc] = ra[iR] * a_cs[iAp1] * pb[iP1];
 
                  iP1 = iP + yOffsetP;
                  rap_cn[iAc] =          a_cn[iA]
                     +          rb[iR] * a_cn[iAm1] * pb[iP1]
                     +          ra[iR] * a_cn[iAp1] * pa[iP1];
 
                  iP1 = iP + xOffsetP;
                  rap_ce[iAc] =          a_ce[iA]
                     +          rb[iR] * a_ce[iAm1] * pb[iP1]
                     +          ra[iR] * a_ce[iAp1] * pa[iP1];
 
               }
            hypre_BoxLoop4End(iP, iR, iA, iAc);

            break;

            /*--------------------------------------------------------------
             * Loop for 15-point fine grid operator; produces upper triangular
             * part of 15-point coarse grid operator. stencil entries:
             * (above-north, above-east, above-center, above-west,
             * above-south, center-north, and center-east).
             *--------------------------------------------------------------*/

            case 15:

            hypre_BoxGetSize(cgrid_box, loop_size);
            hypre_BoxLoop4Begin(loop_size,
                                PT_dbox,  cstart, stridec, iP,
                                R_dbox,   cstart, stridec, iR,
                                A_dbox,   fstart, stridef, iA,
                                RAP_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iP,iR,iA,iAc,iAm1,iAp1,iP1
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop4For(loopi, loopj, loopk, iP, iR, iA, iAc)
               {
                  iAm1 = iA - zOffsetA;
                  iAp1 = iA + zOffsetA;

                  iP1 = iP + zOffsetP + yOffsetP;
                  rap_an[iAc] = ra[iR] * a_cn[iAp1] * pb[iP1]
                     +          ra[iR] * a_an[iAp1]
                     +                   a_an[iA]   * pb[iP1];

                  iP1 = iP + zOffsetP + xOffsetP;
                  rap_ae[iAc] = ra[iR] * a_ce[iAp1] * pb[iP1]
                     +          ra[iR] * a_ae[iAp1]
                     +                   a_ae[iA]   * pb[iP1];
 
                  iP1 = iP + zOffsetP; 
                  rap_ac[iAc] =          a_ac[iA]   * pb[iP1]
                     +          ra[iR] * a_cc[iAp1] * pb[iP1]
                     +          ra[iR] * a_ac[iAp1];
 
                  iP1 = iP + zOffsetP - xOffsetP;
                  rap_aw[iAc] = ra[iR] * a_cw[iAp1] * pb[iP1]
                     +          ra[iR] * a_aw[iAp1]
                     +                   a_aw[iA]   * pb[iP1];
 
                  iP1 = iP + zOffsetP - yOffsetP;
                  rap_as[iAc] = ra[iR] * a_cs[iAp1] * pb[iP1]
                     +          ra[iR] * a_as[iAp1]
                     +                   a_as[iA]   * pb[iP1];
 
                  iP1 = iP + yOffsetP;
                  rap_cn[iAc] =          a_cn[iA]
                     +          rb[iR] * a_cn[iAm1] * pb[iP1]
                     +          ra[iR] * a_cn[iAp1] * pa[iP1]
                     +                   a_bn[iA]   * pb[iP1]
                     +                   a_an[iA]   * pa[iP1]
                     +          rb[iR] * a_an[iAm1]
                     +          ra[iR] * a_bn[iAp1];
 
                  iP1 = iP + xOffsetP;
                  rap_ce[iAc] =          a_ce[iA]
                     +          rb[iR] * a_ce[iAm1] * pb[iP1]
                     +          ra[iR] * a_ce[iAp1] * pa[iP1]
                     +                   a_be[iA]   * pb[iP1]
                     +                   a_ae[iA]   * pa[iP1]
                     +          rb[iR] * a_ae[iAm1]
                     +          ra[iR] * a_be[iAp1];
 
               }
            hypre_BoxLoop4End(iP, iR, iA, iAc);

            break;


            /*--------------------------------------------------------------
             * Loop for 19-point fine grid operator; produces upper triangular
             * part of 27-point coarse grid operator. stencil entries:
             * (above-northeast, above-north, above-northwest, above-east,
             * above-center, above-west, above-southeast, above-south,
             * above-southwest, center-northeast, center-north,
             * center-northwest, and center-east).
             *--------------------------------------------------------------*/

            case 19:

            hypre_BoxGetSize(cgrid_box, loop_size);
            hypre_BoxLoop4Begin(loop_size,
                                PT_dbox,  cstart, stridec, iP,
                                R_dbox,   cstart, stridec, iR,
                                A_dbox,   fstart, stridef, iA,
                                RAP_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iP,iR,iA,iAc,iAm1,iAp1,iP1
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop4For(loopi, loopj, loopk, iP, iR, iA, iAc)
               {
                  iAm1 = iA - zOffsetA;
                  iAp1 = iA + zOffsetA;

                  iP1 = iP + zOffsetP + yOffsetP + xOffsetP;
                  rap_ane[iAc] = ra[iR] * a_cne[iAp1] * pb[iP1];

                  iP1 = iP + zOffsetP + yOffsetP;
                  rap_an[iAc] = ra[iR] * a_cn[iAp1] * pb[iP1]
                     +          ra[iR] * a_an[iAp1]
                     +                   a_an[iA]   * pb[iP1];

                  iP1 = iP + zOffsetP + yOffsetP - xOffsetP;
                  rap_anw[iAc] = ra[iR] * a_cnw[iAp1] * pb[iP1];

                  iP1 = iP + zOffsetP + xOffsetP;
                  rap_ae[iAc] = ra[iR] * a_ce[iAp1] * pb[iP1]
                     +          ra[iR] * a_ae[iAp1]
                     +                   a_ae[iA]   * pb[iP1];

                  iP1 = iP + zOffsetP; 
                  rap_ac[iAc] =          a_ac[iA]   * pb[iP1]
                     +          ra[iR] * a_cc[iAp1] * pb[iP1]
                     +          ra[iR] * a_ac[iAp1];
 
                  iP1 = iP + zOffsetP - xOffsetP;
                  rap_aw[iAc] = ra[iR] * a_cw[iAp1] * pb[iP1]
                     +          ra[iR] * a_aw[iAp1]
                     +                   a_aw[iA]   * pb[iP1];
 
                  iP1 = iP + zOffsetP - yOffsetP + xOffsetP;
                  rap_ase[iAc] = ra[iR] * a_cse[iAp1] * pb[iP1];

                  iP1 = iP + zOffsetP - yOffsetP;
                  rap_as[iAc] = ra[iR] * a_cs[iAp1] * pb[iP1]
                     +          ra[iR] * a_as[iAp1]
                     +                   a_as[iA]   * pb[iP1];
 
                  iP1 = iP + zOffsetP - yOffsetP - xOffsetP;
                  rap_asw[iAc] = ra[iR] * a_csw[iAp1] * pb[iP1];

                  iP1 = iP + yOffsetP + xOffsetP;
                  rap_cne[iAc] =         a_cne[iA]
                     +          rb[iR] * a_cne[iAm1] * pb[iP1]
                     +          ra[iR] * a_cne[iAp1] * pa[iP1];

                  iP1 = iP + yOffsetP;
                  rap_cn[iAc] =          a_cn[iA]
                     +          rb[iR] * a_cn[iAm1] * pb[iP1]
                     +          ra[iR] * a_cn[iAp1] * pa[iP1]
                     +                   a_bn[iA]   * pb[iP1]
                     +                   a_an[iA]   * pa[iP1]
                     +          rb[iR] * a_an[iAm1]
                     +          ra[iR] * a_bn[iAp1];
 
                  iP1 = iP + yOffsetP - xOffsetP;
                  rap_cnw[iAc] =         a_cnw[iA]
                     +          rb[iR] * a_cnw[iAm1] * pb[iP1]
                     +          ra[iR] * a_cnw[iAp1] * pa[iP1];

                  iP1 = iP + xOffsetP;
                  rap_ce[iAc] =          a_ce[iA]
                     +          rb[iR] * a_ce[iAm1] * pb[iP1]
                     +          ra[iR] * a_ce[iAp1] * pa[iP1]
                     +                   a_be[iA]   * pb[iP1]
                     +                   a_ae[iA]   * pa[iP1]
                     +          rb[iR] * a_ae[iAm1]
                     +          ra[iR] * a_be[iAp1];
 
               }
            hypre_BoxLoop4End(iP, iR, iA, iAc);

            break;

            /*--------------------------------------------------------------
             * Loop for 27-point fine grid operator; produces upper triangular
             * part of 27-point coarse grid operator. stencil entries:
             * (above-northeast, above-north, above-northwest, above-east,
             * above-center, above-west, above-southeast, above-south,
             * above-southwest, center-northeast, center-north,
             * center-northwest, and center-east).
             *--------------------------------------------------------------*/

            default:

            hypre_BoxGetSize(cgrid_box, loop_size);
            hypre_BoxLoop4Begin(loop_size,
                                PT_dbox,  cstart, stridec, iP,
                                R_dbox,   cstart, stridec, iR,
                                A_dbox,   fstart, stridef, iA,
                                RAP_dbox, cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iP,iR,iA,iAc,iAm1,iAp1,iP1
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop4For(loopi, loopj, loopk, iP, iR, iA, iAc)
               {
                  iAm1 = iA - zOffsetA;
                  iAp1 = iA + zOffsetA;

                  iP1 = iP + zOffsetP + yOffsetP + xOffsetP;
                  rap_ane[iAc] = ra[iR] * a_cne[iAp1] * pb[iP1]
                     +           ra[iR] * a_ane[iAp1]
                     +                    a_ane[iA]   * pb[iP1];

                  iP1 = iP + zOffsetP + yOffsetP;
                  rap_an[iAc] = ra[iR] * a_cn[iAp1] * pb[iP1]
                     +          ra[iR] * a_an[iAp1]
                     +                   a_an[iA]   * pb[iP1];

                  iP1 = iP + zOffsetP + yOffsetP - xOffsetP;
                  rap_anw[iAc] = ra[iR] * a_cnw[iAp1] * pb[iP1]
                     +           ra[iR] * a_anw[iAp1]
                     +                    a_anw[iA]   * pb[iP1];

                  iP1 = iP + zOffsetP + xOffsetP;
                  rap_ae[iAc] = ra[iR] * a_ce[iAp1] * pb[iP1]
                     +          ra[iR] * a_ae[iAp1]
                     +                   a_ae[iA]   * pb[iP1];

                  iP1 = iP + zOffsetP; 
                  rap_ac[iAc] =          a_ac[iA]   * pb[iP1]
                     +          ra[iR] * a_cc[iAp1] * pb[iP1]
                     +          ra[iR] * a_ac[iAp1];
 
                  iP1 = iP + zOffsetP - xOffsetP;
                  rap_aw[iAc] = ra[iR] * a_cw[iAp1] * pb[iP1]
                     +          ra[iR] * a_aw[iAp1]
                     +                   a_aw[iA]   * pb[iP1];
 
                  iP1 = iP + zOffsetP - yOffsetP + xOffsetP;
                  rap_ase[iAc] = ra[iR] * a_cse[iAp1] * pb[iP1]
                     +           ra[iR] * a_ase[iAp1]
                     +                    a_ase[iA]   * pb[iP1];

                  iP1 = iP + zOffsetP - yOffsetP;
                  rap_as[iAc] = ra[iR] * a_cs[iAp1] * pb[iP1]
                     +          ra[iR] * a_as[iAp1]
                     +                   a_as[iA]   * pb[iP1];
 
                  iP1 = iP + zOffsetP - yOffsetP - xOffsetP;
                  rap_asw[iAc] = ra[iR] * a_csw[iAp1] * pb[iP1]
                     +           ra[iR] * a_asw[iAp1]
                     +                    a_asw[iA]   * pb[iP1];


                  iP1 = iP + yOffsetP + xOffsetP;
                  rap_cne[iAc] =         a_cne[iA]
                     +          rb[iR] * a_cne[iAm1] * pb[iP1]
                     +          ra[iR] * a_cne[iAp1] * pa[iP1]
                     +                   a_bne[iA]   * pb[iP1]
                     +                   a_ane[iA]   * pa[iP1]
                     +          rb[iR] * a_ane[iAm1]
                     +          ra[iR] * a_bne[iAp1];

                  iP1 = iP + yOffsetP;
                  rap_cn[iAc] =          a_cn[iA]
                     +          rb[iR] * a_cn[iAm1] * pb[iP1]
                     +          ra[iR] * a_cn[iAp1] * pa[iP1]
                     +                   a_bn[iA]   * pb[iP1]
                     +                   a_an[iA]   * pa[iP1]
                     +          rb[iR] * a_an[iAm1]
                     +          ra[iR] * a_bn[iAp1];
 
                  iP1 = iP + yOffsetP - xOffsetP;
                  rap_cnw[iAc] =         a_cnw[iA]
                     +          rb[iR] * a_cnw[iAm1] * pb[iP1]
                     +          ra[iR] * a_cnw[iAp1] * pa[iP1]
                     +                   a_bnw[iA]   * pb[iP1]
                     +                   a_anw[iA]   * pa[iP1]
                     +          rb[iR] * a_anw[iAm1]
                     +          ra[iR] * a_bnw[iAp1];

                  iP1 = iP + xOffsetP;
                  rap_ce[iAc] =          a_ce[iA]
                     +          rb[iR] * a_ce[iAm1] * pb[iP1]
                     +          ra[iR] * a_ce[iAp1] * pa[iP1]
                     +                   a_be[iA]   * pb[iP1]
                     +                   a_ae[iA]   * pa[iP1]
                     +          rb[iR] * a_ae[iAm1]
                     +          ra[iR] * a_be[iAp1];
 
               }
            hypre_BoxLoop4End(iP, iR, iA, iAc);

            break;

         } /* end switch statement */

      } /* end ForBoxI */

   return ierr;
}


/*--------------------------------------------------------------------------
 * hypre_SMG3RAPPeriodicSym 
 *    Collapses stencil in periodic direction on coarsest grid.
 *--------------------------------------------------------------------------*/
 
int
hypre_SMG3RAPPeriodicSym( hypre_StructMatrix *RAP,
                          hypre_Index         cindex,
                          hypre_Index         cstride )

{

   hypre_Index             index;

   hypre_StructGrid       *cgrid;
   hypre_BoxArray         *cgrid_boxes;
   hypre_Box              *cgrid_box;
   hypre_IndexRef          cstart;
   hypre_Index             stridec;
   hypre_Index             loop_size;

   int                  ci;
   int                  loopi, loopj, loopk;

   hypre_Box           *RAP_dbox;

   double               *rap_bc, *rap_bw, *rap_be, *rap_bs, *rap_bn;
   double               *rap_cc, *rap_cw,  *rap_cs;
   double               *rap_bsw, *rap_bse, *rap_bnw, *rap_bne;
   double               *rap_csw, *rap_cse;

   int                  iAc;
   int                  iAcmx;
   int                  iAcmy;
   int                  iAcmxmy;
   int                  iAcpxmy;

   int                  xOffset;
   int                  yOffset;

   double               zero = 0.0;

   hypre_StructStencil *stencil;
   int                  stencil_size;

   int                  ierr = 0;

   stencil = hypre_StructMatrixStencil(RAP);
   stencil_size = hypre_StructStencilSize(stencil);

   hypre_SetIndex(stridec, 1, 1, 1);

   cgrid = hypre_StructMatrixGrid(RAP);
   cgrid_boxes = hypre_StructGridBoxes(cgrid);

   if (hypre_IndexZ(hypre_StructGridPeriodic(cgrid)) == 1)
   {
      hypre_StructMatrixAssemble(RAP);

      hypre_ForBoxI(ci, cgrid_boxes)
         {
            cgrid_box = hypre_BoxArrayBox(cgrid_boxes, ci);

            cstart = hypre_BoxIMin(cgrid_box);

            RAP_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(RAP), ci);

            hypre_SetIndex(index,1,0,0);
            xOffset = hypre_BoxOffsetDistance(RAP_dbox,index);
            hypre_SetIndex(index,0,1,0);
            yOffset = hypre_BoxOffsetDistance(RAP_dbox,index);


            /*-----------------------------------------------------------------
             * Extract pointers for 15-point coarse grid operator:
             *-----------------------------------------------------------------*/

            hypre_SetIndex(index,0,0,-1);
            rap_bc = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,-1,0,-1);
            rap_bw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,1,0,-1);
            rap_be = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,-1,-1);
            rap_bs = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,1,-1);
            rap_bn = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,0,0);
            rap_cc = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,-1,0,0);
            rap_cw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,-1,0);
            rap_cs = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            /*-----------------------------------------------------------------
             * Extract additional pointers for 27-point coarse grid operator:
             *-----------------------------------------------------------------*/

            if(stencil_size == 27)
            {
               hypre_SetIndex(index,-1,-1,-1);
               rap_bsw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,1,-1,-1);
               rap_bse = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,-1,1,-1);
               rap_bnw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,1,1,-1);
               rap_bne = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,-1,-1,0);
               rap_csw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,1,-1,0);
               rap_cse = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);
            }

            /*-----------------------------------------------------------------
             * Collapse 15 point operator.
             *-----------------------------------------------------------------*/

            hypre_BoxGetSize(cgrid_box, loop_size);

            hypre_BoxLoop1Begin(loop_size,
                                RAP_dbox,  cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iAc,iAcmx,iAcmy
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop1For(loopi, loopj, loopk, iAc)
               {
                  iAcmx = iAc - xOffset;
                  iAcmy = iAc - yOffset;
               
                  rap_cc[iAc] += (2.0 * rap_bc[iAc]);
                  rap_cw[iAc] += (rap_bw[iAc] + rap_be[iAcmx]);
                  rap_cs[iAc] += (rap_bs[iAc] + rap_bn[iAcmy]);
               }
            hypre_BoxLoop1End(iAc);

            hypre_BoxLoop1Begin(loop_size,
                                RAP_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)
               {
                  rap_bc[iAc]  = zero;
                  rap_bw[iAc]  = zero;
                  rap_be[iAc]  = zero;
                  rap_bs[iAc]  = zero;
                  rap_bn[iAc]  = zero;
               }
            hypre_BoxLoop1End(iAc);

            /*-----------------------------------------------------------------
             * Collapse additional entries for 27 point operator.
             *-----------------------------------------------------------------*/

            if (stencil_size == 27)
            {
               hypre_BoxGetSize(cgrid_box, loop_size);

               hypre_BoxLoop1Begin(loop_size,
                                   RAP_dbox,  cstart, stridec, iAc);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,iAc,iAcmxmy,iAcpxmy
#include "hypre_box_smp_forloop.h"
               hypre_BoxLoop1For(loopi, loopj, loopk, iAc)
                  {
                     iAcmxmy = iAc - xOffset - yOffset;
                     iAcpxmy = iAc + xOffset - yOffset;
                  
                     rap_csw[iAc] += (rap_bsw[iAc] + rap_bne[iAcmxmy]);
                  
                     rap_cse[iAc] += (rap_bse[iAc] + rap_bnw[iAcpxmy]);
                  
                  }
               hypre_BoxLoop1End(iAc);

               hypre_BoxLoop1Begin(loop_size,
                                   RAP_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)
                  {
                     rap_bsw[iAc]  = zero;
                     rap_bse[iAc]  = zero;
                     rap_bnw[iAc]  = zero;
                     rap_bne[iAc]  = zero;
                  }
               hypre_BoxLoop1End(iAc);
            }

         } /* end ForBoxI */

   }
   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMG3RAPPeriodicNoSym 
 *    Collapses stencil in periodic direction on coarsest grid.
 *--------------------------------------------------------------------------*/
 
int
hypre_SMG3RAPPeriodicNoSym( hypre_StructMatrix *RAP,
                            hypre_Index         cindex,
                            hypre_Index         cstride )

{

   hypre_Index             index;

   hypre_StructGrid       *cgrid;
   hypre_BoxArray         *cgrid_boxes;
   hypre_Box              *cgrid_box;
   hypre_IndexRef          cstart;
   hypre_Index             stridec;
   hypre_Index             loop_size;

   int                  ci;
   int                  loopi, loopj, loopk;

   hypre_Box           *RAP_dbox;

   double               *rap_bc, *rap_bw, *rap_be, *rap_bs, *rap_bn;
   double               *rap_cc, *rap_cw, *rap_ce, *rap_cs, *rap_cn;
   double               *rap_ac, *rap_aw, *rap_ae, *rap_as, *rap_an;
   double               *rap_bsw, *rap_bse, *rap_bnw, *rap_bne;
   double               *rap_csw, *rap_cse, *rap_cnw, *rap_cne;
   double               *rap_asw, *rap_ase, *rap_anw, *rap_ane;

   int                  iAc;

   double               zero = 0.0;

   hypre_StructStencil *stencil;
   int                  stencil_size;

   int                  ierr = 0;

   stencil = hypre_StructMatrixStencil(RAP);
   stencil_size = hypre_StructStencilSize(stencil);

   hypre_SetIndex(stridec, 1, 1, 1);

   cgrid = hypre_StructMatrixGrid(RAP);
   cgrid_boxes = hypre_StructGridBoxes(cgrid);

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

            cstart = hypre_BoxIMin(cgrid_box);

            RAP_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(RAP), ci);

            /*-----------------------------------------------------------------
             * Extract pointers for 15-point coarse grid operator:
             *-----------------------------------------------------------------*/

            hypre_SetIndex(index,0,0,-1);
            rap_bc = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,-1,0,-1);
            rap_bw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,1,0,-1);
            rap_be = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,-1,-1);
            rap_bs = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,1,-1);
            rap_bn = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,0,0);
            rap_cc = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,-1,0,0);
            rap_cw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,1,0,0);
            rap_ce = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,-1,0);
            rap_cs = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,1,0);
            rap_cn = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,0,1);
            rap_ac = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,-1,0,1);
            rap_aw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,1,0,1);
            rap_ae = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,-1,1);
            rap_as = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            hypre_SetIndex(index,0,1,1);
            rap_an = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            /*-----------------------------------------------------------------
             * Extract additional pointers for 27-point coarse grid operator:
             *-----------------------------------------------------------------*/

            if(stencil_size == 27)
            {

               hypre_SetIndex(index,-1,-1,-1);
               rap_bsw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,1,-1,-1);
               rap_bse = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,-1,1,-1);
               rap_bnw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,1,1,-1);
               rap_bne = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,-1,-1,0);
               rap_csw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,1,-1,0);
               rap_cse = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,-1,1,0);
               rap_cnw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,1,1,0);
               rap_cne = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,-1,-1,1);
               rap_asw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,1,-1,1);
               rap_ase = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,-1,1,1);
               rap_anw = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

               hypre_SetIndex(index,1,1,1);
               rap_ane = hypre_StructMatrixExtractPointerByIndex(RAP, ci, index);

            }

            /*-----------------------------------------------------------------
             * Collapse 15 point operator.
             *-----------------------------------------------------------------*/

            hypre_BoxGetSize(cgrid_box, loop_size);

            hypre_BoxLoop1Begin(loop_size,
                                RAP_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)
               {
                  rap_cc[iAc] += (rap_bc[iAc] + rap_ac[iAc]);
                  rap_bc[iAc]  = zero;
                  rap_ac[iAc]  = zero;
               
                  rap_cw[iAc] += (rap_bw[iAc] + rap_aw[iAc]);
                  rap_bw[iAc]  = zero;
                  rap_aw[iAc]  = zero;
               
                  rap_ce[iAc] += (rap_be[iAc] + rap_ae[iAc]);
                  rap_be[iAc]  = zero;
                  rap_ae[iAc]  = zero;
               
                  rap_cs[iAc] += (rap_bs[iAc] + rap_as[iAc]);
                  rap_bs[iAc]  = zero;
                  rap_as[iAc]  = zero;
               
                  rap_cn[iAc] += (rap_bn[iAc] + rap_an[iAc]);
                  rap_bn[iAc]  = zero;
                  rap_an[iAc]  = zero;
               }
            hypre_BoxLoop1End(iAc);

            /*-----------------------------------------------------------------
             * Collapse additional entries for 27 point operator.
             *-----------------------------------------------------------------*/

            if (stencil_size == 27)
            {
               hypre_BoxGetSize(cgrid_box, loop_size);

               hypre_BoxLoop1Begin(loop_size,
                                   RAP_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)
                  {
                     rap_csw[iAc] += (rap_bsw[iAc] + rap_asw[iAc]);
                     rap_bsw[iAc]  = zero;
                     rap_asw[iAc]  = zero;
                  
                     rap_cse[iAc] += (rap_bse[iAc] + rap_ase[iAc]);
                     rap_bse[iAc]  = zero;
                     rap_ase[iAc]  = zero;
                  
                     rap_cnw[iAc] += (rap_bnw[iAc] + rap_anw[iAc]);
                     rap_bnw[iAc]  = zero;
                     rap_anw[iAc]  = zero;
                  
                     rap_cne[iAc] += (rap_bne[iAc] + rap_ane[iAc]);
                     rap_bne[iAc]  = zero;
                     rap_ane[iAc]  = zero;            
                  }
               hypre_BoxLoop1End(iAc);
            }

         } /* end ForBoxI */

   }
   return ierr;
}