/*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*/
/******************************************************************************
 *
 * Structured matrix-vector multiply routine
 *
 *****************************************************************************/

#include "headers.h"

/* this currently cannot be greater than 7 */
#ifdef MAX_DEPTH
#undef MAX_DEPTH
#endif
#define MAX_DEPTH 7

/*--------------------------------------------------------------------------
 * hypre_StructMatvecData data structure
 *--------------------------------------------------------------------------*/

typedef struct
{
   hypre_StructMatrix  *A;
   hypre_StructVector  *x;
   hypre_ComputePkg    *compute_pkg;

} hypre_StructMatvecData;

/*--------------------------------------------------------------------------
 * hypre_StructMatvecCreate
 *--------------------------------------------------------------------------*/

void *
hypre_StructMatvecCreate( )
{
   hypre_StructMatvecData *matvec_data;

   matvec_data = hypre_CTAlloc(hypre_StructMatvecData, 1);

   return (void *) matvec_data;
}

/*--------------------------------------------------------------------------
 * hypre_StructMatvecSetup
 *--------------------------------------------------------------------------*/

int
hypre_StructMatvecSetup( void               *matvec_vdata,
                         hypre_StructMatrix *A,
                         hypre_StructVector *x            )
{
   int ierr = 0;

   hypre_StructMatvecData  *matvec_data = matvec_vdata;
                          
   hypre_StructGrid        *grid;
   hypre_StructStencil     *stencil;
                          
   hypre_BoxArrayArray     *send_boxes;
   hypre_BoxArrayArray     *recv_boxes;
   int                    **send_processes;
   int                    **recv_processes;
   hypre_BoxArrayArray     *indt_boxes;
   hypre_BoxArrayArray     *dept_boxes;

   hypre_Index              unit_stride;
                       
   hypre_ComputePkg        *compute_pkg;

   /*----------------------------------------------------------
    * Set up the compute package
    *----------------------------------------------------------*/

   grid    = hypre_StructMatrixGrid(A);
   stencil = hypre_StructMatrixStencil(A);

   hypre_CreateComputeInfo(grid, stencil,
                        &send_boxes, &recv_boxes,
                        &send_processes, &recv_processes,
                        &indt_boxes, &dept_boxes);

   hypre_SetIndex(unit_stride, 1, 1, 1);
   hypre_ComputePkgCreate(send_boxes, recv_boxes,
                          unit_stride, unit_stride,
                          send_processes, recv_processes,
                          indt_boxes, dept_boxes,
                          unit_stride,
                          grid, hypre_StructVectorDataSpace(x), 1,
                          &compute_pkg);

   /*----------------------------------------------------------
    * Set up the matvec data structure
    *----------------------------------------------------------*/

   (matvec_data -> A)           = hypre_StructMatrixRef(A);
   (matvec_data -> x)           = hypre_StructVectorRef(x);
   (matvec_data -> compute_pkg) = compute_pkg;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_StructMatvecCompute
 *--------------------------------------------------------------------------*/

int
hypre_StructMatvecCompute( void               *matvec_vdata,
                           double              alpha,
                           hypre_StructMatrix *A,
                           hypre_StructVector *x,
                           double              beta,
                           hypre_StructVector *y            )
{
   int ierr = 0;

   hypre_StructMatvecData  *matvec_data = matvec_vdata;
                          
   hypre_ComputePkg        *compute_pkg;
                          
   hypre_CommHandle        *comm_handle;
                          
   hypre_BoxArrayArray     *compute_box_aa;
   hypre_BoxArray          *compute_box_a;
   hypre_Box               *compute_box;
                          
   hypre_Box               *A_data_box;
   hypre_Box               *x_data_box;
   hypre_Box               *y_data_box;
                          
   int                      Ai;
   int                      xi;
   int                      xoff0;
   int                      xoff1;
   int                      xoff2;
   int                      xoff3;
   int                      xoff4;
   int                      xoff5;
   int                      xoff6;
   int                      yi;
                          
   double                  *Ap0;
   double                  *Ap1;
   double                  *Ap2;
   double                  *Ap3;
   double                  *Ap4;
   double                  *Ap5;
   double                  *Ap6;
   double                  *xp;
   double                  *yp;
                          
   hypre_BoxArray          *boxes;
   hypre_Box               *box;
   hypre_Index              loop_size;
   hypre_IndexRef           start;
   hypre_IndexRef           stride;
                          
   hypre_StructStencil     *stencil;
   hypre_Index             *stencil_shape;
   int                      stencil_size;
   int                      depth;
                          
   double                   temp;
   int                      compute_i, i, j, si;
   int                      loopi, loopj, loopk;

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

   compute_pkg = (matvec_data -> compute_pkg);

   stride = hypre_ComputePkgStride(compute_pkg);

   /*-----------------------------------------------------------------------
    * Do (alpha == 0.0) computation
    *-----------------------------------------------------------------------*/

   if (alpha == 0.0)
   {
      boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(A));
      hypre_ForBoxI(i, boxes)
         {
            box   = hypre_BoxArrayBox(boxes, i);
            start = hypre_BoxIMin(box);

            y_data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(y), i);
            yp = hypre_StructVectorBoxData(y, i);

            hypre_BoxGetSize(box, loop_size);

            hypre_BoxLoop1Begin(loop_size,
                                y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi
#include "hypre_box_smp_forloop.h"
            hypre_BoxLoop1For(loopi, loopj, loopk, yi)
               {
                  yp[yi] *= beta;
               }
            hypre_BoxLoop1End(yi);
         }

      return ierr;
   }

   /*-----------------------------------------------------------------------
    * Do (alpha != 0.0) computation
    *-----------------------------------------------------------------------*/

   stencil       = hypre_StructMatrixStencil(A);
   stencil_shape = hypre_StructStencilShape(stencil);
   stencil_size  = hypre_StructStencilSize(stencil);

   for (compute_i = 0; compute_i < 2; compute_i++)
   {
      switch(compute_i)
      {
         case 0:
         {
            xp = hypre_StructVectorData(x);
            hypre_InitializeIndtComputations(compute_pkg, xp, &comm_handle);
            compute_box_aa = hypre_ComputePkgIndtBoxes(compute_pkg);

            /*--------------------------------------------------------------
             * initialize y= (beta/alpha)*y
             *--------------------------------------------------------------*/

            temp = beta / alpha;
            if (temp != 1.0)
            {
               boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(A));
               hypre_ForBoxI(i, boxes)
                  {
                     box   = hypre_BoxArrayBox(boxes, i);
                     start = hypre_BoxIMin(box);

                     y_data_box =
                        hypre_BoxArrayBox(hypre_StructVectorDataSpace(y), i);
                     yp = hypre_StructVectorBoxData(y, i);

                     if (temp == 0.0)
                     {
                        hypre_BoxGetSize(box, loop_size);

                        hypre_BoxLoop1Begin(loop_size,
                                            y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi
#include "hypre_box_smp_forloop.h"
                        hypre_BoxLoop1For(loopi, loopj, loopk, yi)
                           {
                              yp[yi] = 0.0;
                           }
                        hypre_BoxLoop1End(yi);
                     }
                     else
                     {
                        hypre_BoxGetSize(box, loop_size);

                        hypre_BoxLoop1Begin(loop_size,
                                            y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi
#include "hypre_box_smp_forloop.h"
                        hypre_BoxLoop1For(loopi, loopj, loopk, yi)
                           {
                              yp[yi] *= temp;
                           }
                        hypre_BoxLoop1End(yi);
                     }
                  }
            }
         }
         break;

         case 1:
         {
            hypre_FinalizeIndtComputations(comm_handle);
            compute_box_aa = hypre_ComputePkgDeptBoxes(compute_pkg);
         }
         break;
      }

      /*--------------------------------------------------------------------
       * y += A*x
       *--------------------------------------------------------------------*/

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

            A_data_box = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(A), i);
            x_data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(x), i);
            y_data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(y), i);

            xp = hypre_StructVectorBoxData(x, i);
            yp = hypre_StructVectorBoxData(y, i);

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

                  hypre_BoxGetSize(compute_box, loop_size);
                  start  = hypre_BoxIMin(compute_box);

                  /* unroll up to depth MAX_DEPTH */
                  for (si = 0; si < stencil_size; si+= MAX_DEPTH)
                  {
                     depth = hypre_min(MAX_DEPTH, (stencil_size -si));
                     switch(depth)
                     {
                        case 7:
                        Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
                        Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
                        Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
                        Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
                        Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
                        Ap5 = hypre_StructMatrixBoxData(A, i, si+5);
                        Ap6 = hypre_StructMatrixBoxData(A, i, si+6);

                        xoff0 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+0]);
                        xoff1 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+1]);
                        xoff2 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+2]);
                        xoff3 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+3]);
                        xoff4 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+4]);
                        xoff5 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+5]);
                        xoff6 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+6]);

                        hypre_BoxLoop3Begin(loop_size,
                                            A_data_box, start, stride, Ai,
                                            x_data_box, start, stride, xi,
                                            y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi,xi,Ai
#include "hypre_box_smp_forloop.h"
                        hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
                           {
                              yp[yi] +=
                                 Ap0[Ai] * xp[xi + xoff0] +
                                 Ap1[Ai] * xp[xi + xoff1] +
                                 Ap2[Ai] * xp[xi + xoff2] +
                                 Ap3[Ai] * xp[xi + xoff3] +
                                 Ap4[Ai] * xp[xi + xoff4] +
                                 Ap5[Ai] * xp[xi + xoff5] +
                                 Ap6[Ai] * xp[xi + xoff6];
                           }
                        hypre_BoxLoop3End(Ai, xi, yi);

                        break;

                        case 6:
                        Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
                        Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
                        Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
                        Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
                        Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
                        Ap5 = hypre_StructMatrixBoxData(A, i, si+5);

                        xoff0 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+0]);
                        xoff1 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+1]);
                        xoff2 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+2]);
                        xoff3 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+3]);
                        xoff4 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+4]);
                        xoff5 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+5]);

                        hypre_BoxLoop3Begin(loop_size,
                                            A_data_box, start, stride, Ai,
                                            x_data_box, start, stride, xi,
                                            y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi,xi,Ai
#include "hypre_box_smp_forloop.h"
                        hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
                           {
                              yp[yi] +=
                                 Ap0[Ai] * xp[xi + xoff0] +
                                 Ap1[Ai] * xp[xi + xoff1] +
                                 Ap2[Ai] * xp[xi + xoff2] +
                                 Ap3[Ai] * xp[xi + xoff3] +
                                 Ap4[Ai] * xp[xi + xoff4] +
                                 Ap5[Ai] * xp[xi + xoff5];
                           }
                        hypre_BoxLoop3End(Ai, xi, yi);
                        
                        break;

                        case 5:
                        Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
                        Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
                        Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
                        Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
                        Ap4 = hypre_StructMatrixBoxData(A, i, si+4);

                        xoff0 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+0]);
                        xoff1 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+1]);
                        xoff2 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+2]);
                        xoff3 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+3]);
                        xoff4 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+4]);

                        hypre_BoxLoop3Begin(loop_size,
                                            A_data_box, start, stride, Ai,
                                            x_data_box, start, stride, xi,
                                            y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi,xi,Ai
#include "hypre_box_smp_forloop.h"
                        hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
                           {
                              yp[yi] +=
                                 Ap0[Ai] * xp[xi + xoff0] +
                                 Ap1[Ai] * xp[xi + xoff1] +
                                 Ap2[Ai] * xp[xi + xoff2] +
                                 Ap3[Ai] * xp[xi + xoff3] +
                                 Ap4[Ai] * xp[xi + xoff4]; 
                           }
                        hypre_BoxLoop3End(Ai, xi, yi);

                        break;

                        case 4:
                        Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
                        Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
                        Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
                        Ap3 = hypre_StructMatrixBoxData(A, i, si+3);

                        xoff0 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+0]);
                        xoff1 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+1]);
                        xoff2 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+2]);
                        xoff3 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+3]);

                        hypre_BoxLoop3Begin(loop_size,
                                            A_data_box, start, stride, Ai,
                                            x_data_box, start, stride, xi,
                                            y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi,xi,Ai
#include "hypre_box_smp_forloop.h"
                        hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
                           {
                              yp[yi] +=
                                 Ap0[Ai] * xp[xi + xoff0] +
                                 Ap1[Ai] * xp[xi + xoff1] +
                                 Ap2[Ai] * xp[xi + xoff2] +
                                 Ap3[Ai] * xp[xi + xoff3];
                           }
                        hypre_BoxLoop3End(Ai, xi, yi);

                        break;

                        case 3:
                        Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
                        Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
                        Ap2 = hypre_StructMatrixBoxData(A, i, si+2);

                        xoff0 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+0]);
                        xoff1 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+1]);
                        xoff2 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+2]);

                        hypre_BoxLoop3Begin(loop_size,
                                            A_data_box, start, stride, Ai,
                                            x_data_box, start, stride, xi,
                                            y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi,xi,Ai
#include "hypre_box_smp_forloop.h"
                        hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
                           {
                              yp[yi] +=
                                 Ap0[Ai] * xp[xi + xoff0] +
                                 Ap1[Ai] * xp[xi + xoff1] +
                                 Ap2[Ai] * xp[xi + xoff2]; 
                           }
                        hypre_BoxLoop3End(Ai, xi, yi);

                        break;

                        case 2:
                        Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
                        Ap1 = hypre_StructMatrixBoxData(A, i, si+1);

                        xoff0 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+0]);
                        xoff1 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+1]);

                        hypre_BoxLoop3Begin(loop_size,
                                            A_data_box, start, stride, Ai,
                                            x_data_box, start, stride, xi,
                                            y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi,xi,Ai
#include "hypre_box_smp_forloop.h"
                        hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
                           {
                              yp[yi] +=
                                 Ap0[Ai] * xp[xi + xoff0] +
                                 Ap1[Ai] * xp[xi + xoff1];
                           }
                        hypre_BoxLoop3End(Ai, xi, yi);

                        break;

                        case 1:
                        Ap0 = hypre_StructMatrixBoxData(A, i, si+0);

                        xoff0 = hypre_BoxOffsetDistance(x_data_box,
                                                        stencil_shape[si+0]);

                        hypre_BoxLoop3Begin(loop_size,
                                            A_data_box, start, stride, Ai,
                                            x_data_box, start, stride, xi,
                                            y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi,xi,Ai
#include "hypre_box_smp_forloop.h"
                        hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
                           {
                              yp[yi] +=
                                 Ap0[Ai] * xp[xi + xoff0];
                           }
                        hypre_BoxLoop3End(Ai, xi, yi);

                        break;
                     }
                  }

                  if (alpha != 1.0)
                  {
                     hypre_BoxLoop1Begin(loop_size,
                                         y_data_box, start, stride, yi);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,yi
#include "hypre_box_smp_forloop.h"
                     hypre_BoxLoop1For(loopi, loopj, loopk, yi)
                        {
                           yp[yi] *= alpha;
                        }
                     hypre_BoxLoop1End(yi);
                  }
               }
         }
   }
   
   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_StructMatvecDestroy
 *--------------------------------------------------------------------------*/

int
hypre_StructMatvecDestroy( void *matvec_vdata )
{
   int ierr = 0;

   hypre_StructMatvecData *matvec_data = matvec_vdata;

   if (matvec_data)
   {
      hypre_StructMatrixDestroy(matvec_data -> A);
      hypre_StructVectorDestroy(matvec_data -> x);
      hypre_ComputePkgDestroy(matvec_data -> compute_pkg );
      hypre_TFree(matvec_data);
   }

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_StructMatvec
 *--------------------------------------------------------------------------*/

int
hypre_StructMatvec( double              alpha,
                    hypre_StructMatrix *A,
                    hypre_StructVector *x,
                    double              beta,
                    hypre_StructVector *y     )
{
   int ierr = 0;

   void *matvec_data;

   matvec_data = hypre_StructMatvecCreate();
   ierr = hypre_StructMatvecSetup(matvec_data, A, x);
   ierr = hypre_StructMatvecCompute(matvec_data, alpha, A, x, beta, y);
   ierr = hypre_StructMatvecDestroy(matvec_data);

   return ierr;
}