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

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxData data structure
 *--------------------------------------------------------------------------*/

typedef struct
{
   int                     setup_temp_vec;
   int                     setup_a_rem;
   int                     setup_a_sol;
                       
   MPI_Comm                comm;
                       
   int                     memory_use;
   double                  tol;
   int                     max_iter;
   int                     zero_guess;
                         
   int                     num_spaces;
   int                    *space_indices;
   int                    *space_strides;
                       
   int                     num_pre_spaces;
   int                     num_reg_spaces;
   int                    *pre_space_ranks;
   int                    *reg_space_ranks;

   hypre_Index             base_index;
   hypre_Index             base_stride;
   hypre_BoxArray         *base_box_array;

   int                     stencil_dim;
                       
   hypre_StructMatrix     *A;
   hypre_StructVector     *b;
   hypre_StructVector     *x;

   hypre_StructVector     *temp_vec;
   hypre_StructMatrix     *A_sol;  /* Coefficients of A that make up
                                      the (sol)ve part of the relaxation */
   hypre_StructMatrix     *A_rem;  /* Coefficients of A (rem)aining:
                                      A_rem = A - A_sol                  */
   void                  **residual_data;  /* Array of size `num_spaces' */
   void                  **solve_data;     /* Array of size `num_spaces' */

   /* log info (always logged) */
   int                     num_iterations;
   int                     time_index;
                         
   int                     num_pre_relax;
   int                     num_post_relax;

} hypre_SMGRelaxData;

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxCreate
 *--------------------------------------------------------------------------*/

void *
hypre_SMGRelaxCreate( MPI_Comm  comm )
{
   hypre_SMGRelaxData *relax_data;

   relax_data = hypre_CTAlloc(hypre_SMGRelaxData, 1);
   (relax_data -> setup_temp_vec) = 1;
   (relax_data -> setup_a_rem)    = 1;
   (relax_data -> setup_a_sol)    = 1;
   (relax_data -> comm)           = comm;
   (relax_data -> base_box_array) = NULL;
   (relax_data -> time_index)     = hypre_InitializeTiming("SMGRelax");
   /* set defaults */
   (relax_data -> memory_use)         = 0;
   (relax_data -> tol)                = 1.0e-06;
   (relax_data -> max_iter)           = 1000;
   (relax_data -> zero_guess)         = 0;
   (relax_data -> num_spaces)         = 1;
   (relax_data -> space_indices)      = hypre_TAlloc(int, 1);
   (relax_data -> space_strides)      = hypre_TAlloc(int, 1);
   (relax_data -> space_indices[0])   = 0;
   (relax_data -> space_strides[0])   = 1;
   (relax_data -> num_pre_spaces)     = 0;
   (relax_data -> num_reg_spaces)     = 1;
   (relax_data -> pre_space_ranks)    = NULL;
   (relax_data -> reg_space_ranks)    = hypre_TAlloc(int, 1);
   (relax_data -> reg_space_ranks[0]) = 0;
   hypre_SetIndex((relax_data -> base_index), 0, 0, 0);
   hypre_SetIndex((relax_data -> base_stride), 1, 1, 1);
   (relax_data -> A)                  = NULL;
   (relax_data -> b)                  = NULL;
   (relax_data -> x)                  = NULL;
   (relax_data -> temp_vec)           = NULL;

   (relax_data -> num_pre_relax) = 1;
   (relax_data -> num_post_relax) = 1;

   return (void *) relax_data;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxDestroyTempVec
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxDestroyTempVec( void *relax_vdata )
{
   hypre_SMGRelaxData  *relax_data = relax_vdata;
   int                  ierr = 0;

   hypre_StructVectorDestroy(relax_data -> temp_vec);
   (relax_data -> setup_temp_vec) = 1;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxDestroyARem
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxDestroyARem( void *relax_vdata )
{
   hypre_SMGRelaxData  *relax_data = relax_vdata;
   int                  i;
   int                  ierr = 0;

   if (relax_data -> A_rem)
   {
      for (i = 0; i < (relax_data -> num_spaces); i++)
      {
         hypre_SMGResidualDestroy(relax_data -> residual_data[i]);
      }
      hypre_TFree(relax_data -> residual_data);
      hypre_StructMatrixDestroy(relax_data -> A_rem);
      (relax_data -> A_rem) = NULL;
   }
   (relax_data -> setup_a_rem) = 1;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxDestroyASol
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxDestroyASol( void *relax_vdata )
{
   hypre_SMGRelaxData  *relax_data = relax_vdata;
   int                  stencil_dim;
   int                  i;
   int                  ierr = 0;

   if (relax_data -> A_sol)
   {
      stencil_dim = (relax_data -> stencil_dim);
      for (i = 0; i < (relax_data -> num_spaces); i++)
      {
         if (stencil_dim > 2)
            hypre_SMGDestroy(relax_data -> solve_data[i]);
         else
            hypre_CyclicReductionDestroy(relax_data -> solve_data[i]);
      }
      hypre_TFree(relax_data -> solve_data);
      hypre_StructMatrixDestroy(relax_data -> A_sol);
      (relax_data -> A_sol) = NULL;
   }
   (relax_data -> setup_a_sol) = 1;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxDestroy
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxDestroy( void *relax_vdata )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 ierr = 0;

   if (relax_data)
   {
      hypre_TFree(relax_data -> space_indices);
      hypre_TFree(relax_data -> space_strides);
      hypre_TFree(relax_data -> pre_space_ranks);
      hypre_TFree(relax_data -> reg_space_ranks);
      hypre_BoxArrayDestroy(relax_data -> base_box_array);

      hypre_StructMatrixDestroy(relax_data -> A);
      hypre_StructVectorDestroy(relax_data -> b);
      hypre_StructVectorDestroy(relax_data -> x);

      hypre_SMGRelaxDestroyTempVec(relax_vdata);
      hypre_SMGRelaxDestroyARem(relax_vdata);
      hypre_SMGRelaxDestroyASol(relax_vdata);

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

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelax
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelax( void               *relax_vdata,
                hypre_StructMatrix *A,
                hypre_StructVector *b,
                hypre_StructVector *x           )
{
   hypre_SMGRelaxData   *relax_data = relax_vdata;

   int                   zero_guess;
   int                   stencil_dim;
   hypre_StructVector   *temp_vec;
   hypre_StructMatrix   *A_sol;
   hypre_StructMatrix   *A_rem;
   void                **residual_data;
   void                **solve_data;

   hypre_IndexRef        base_stride;
   hypre_BoxArray       *base_box_a;
   double                zero = 0.0;

   int                   max_iter;
   int                   num_spaces;
   int                  *space_ranks;
                    
   int                   i, j, k, is;
                    
   int                   ierr = 0;

   /*----------------------------------------------------------
    * Note: The zero_guess stuff is not handled correctly
    * for general relaxation parameters.  It is correct when
    * the spaces are independent sets in the direction of
    * relaxation.
    *----------------------------------------------------------*/

   hypre_BeginTiming(relax_data -> time_index);

   /*----------------------------------------------------------
    * Set up the solver
    *----------------------------------------------------------*/

   /* insure that the solver memory gets fully set up */
   if ((relax_data -> setup_a_sol) > 0)
   {
      (relax_data -> setup_a_sol) = 2;
   }

   hypre_SMGRelaxSetup(relax_vdata, A, b, x);

   zero_guess      = (relax_data -> zero_guess);
   stencil_dim     = (relax_data -> stencil_dim);
   temp_vec        = (relax_data -> temp_vec);
   A_sol           = (relax_data -> A_sol);
   A_rem           = (relax_data -> A_rem);
   residual_data   = (relax_data -> residual_data);
   solve_data      = (relax_data -> solve_data);


   /*----------------------------------------------------------
    * Set zero values
    *----------------------------------------------------------*/

   if (zero_guess)
   {
      base_stride = (relax_data -> base_stride);
      base_box_a = (relax_data -> base_box_array);
      ierr = hypre_SMGSetStructVectorConstantValues(x, zero, base_box_a,
                                                    base_stride); 
   }

   /*----------------------------------------------------------
    * Iterate
    *----------------------------------------------------------*/

   for (k = 0; k < 2; k++)
   {
      switch(k)
      {
         /* Do pre-relaxation iterations */
         case 0:
         max_iter    = 1;
         num_spaces  = (relax_data -> num_pre_spaces);
         space_ranks = (relax_data -> pre_space_ranks);
         break;

         /* Do regular relaxation iterations */
         case 1:
         max_iter    = (relax_data -> max_iter);
         num_spaces  = (relax_data -> num_reg_spaces);
         space_ranks = (relax_data -> reg_space_ranks);
         break;
      }

      for (i = 0; i < max_iter; i++)
      {
         for (j = 0; j < num_spaces; j++)
         {
            is = space_ranks[j];

            hypre_SMGResidual(residual_data[is], A_rem, x, b, temp_vec);

            if (stencil_dim > 2)
               hypre_SMGSolve(solve_data[is], A_sol, temp_vec, x);
            else
               hypre_CyclicReduction(solve_data[is], A_sol, temp_vec, x);
         }

         (relax_data -> num_iterations) = (i + 1);
      }
   }

   /*----------------------------------------------------------
    * Free up memory according to memory_use parameter
    *----------------------------------------------------------*/

   if ((stencil_dim - 1) <= (relax_data -> memory_use))
   {
      hypre_SMGRelaxDestroyASol(relax_vdata);
   }

   hypre_EndTiming(relax_data -> time_index);

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetup
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetup( void               *relax_vdata,
                     hypre_StructMatrix *A,
                     hypre_StructVector *b,
                     hypre_StructVector *x           )
{
   hypre_SMGRelaxData  *relax_data = relax_vdata;
   int                  stencil_dim;
   int                  a_sol_test;
   int                  ierr = 0;

   stencil_dim = hypre_StructStencilDim(hypre_StructMatrixStencil(A));
   (relax_data -> stencil_dim) = stencil_dim;
   hypre_StructMatrixDestroy(relax_data -> A);
   hypre_StructVectorDestroy(relax_data -> b);
   hypre_StructVectorDestroy(relax_data -> x);
   (relax_data -> A) = hypre_StructMatrixRef(A);
   (relax_data -> b) = hypre_StructVectorRef(b);
   (relax_data -> x) = hypre_StructVectorRef(x);

   /*----------------------------------------------------------
    * Set up memory according to memory_use parameter.
    *
    * If a subset of the solver memory is not to be set up
    * until the solve is actually done, it's "setup" tag
    * should have a value greater than 1.
    *----------------------------------------------------------*/

   if ((stencil_dim - 1) <= (relax_data -> memory_use))
   {
      a_sol_test = 1;
   }
   else
   {
      a_sol_test = 0;
   }

   /*----------------------------------------------------------
    * Set up the solver
    *----------------------------------------------------------*/

   if ((relax_data -> setup_temp_vec) > 0)
   {
      ierr = hypre_SMGRelaxSetupTempVec(relax_vdata, A, b, x);
   }

   if ((relax_data -> setup_a_rem) > 0)
   {
      ierr = hypre_SMGRelaxSetupARem(relax_vdata, A, b, x);
   }

   if ((relax_data -> setup_a_sol) > a_sol_test)
   {
      ierr = hypre_SMGRelaxSetupASol(relax_vdata, A, b, x);
   }

   if ((relax_data -> base_box_array) == NULL)
   {
      ierr = hypre_SMGRelaxSetupBaseBoxArray(relax_vdata, A, b, x);
   }
   

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetupTempVec
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetupTempVec( void               *relax_vdata,
                            hypre_StructMatrix *A,
                            hypre_StructVector *b,
                            hypre_StructVector *x           )
{
   hypre_SMGRelaxData  *relax_data = relax_vdata;
   hypre_StructVector  *temp_vec   = (relax_data -> temp_vec);
   int                  ierr = 0;

   /*----------------------------------------------------------
    * Set up data
    *----------------------------------------------------------*/

   if ((relax_data -> temp_vec) == NULL)
   {
      temp_vec = hypre_StructVectorCreate(hypre_StructVectorComm(b),
                                          hypre_StructVectorGrid(b));
      hypre_StructVectorSetNumGhost(temp_vec, hypre_StructVectorNumGhost(b));
      hypre_StructVectorInitialize(temp_vec);
      hypre_StructVectorAssemble(temp_vec);
      (relax_data -> temp_vec) = temp_vec;
   }
   (relax_data -> setup_temp_vec) = 0;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetupARem
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetupARem( void               *relax_vdata,
                         hypre_StructMatrix *A,
                         hypre_StructVector *b,
                         hypre_StructVector *x           )
{
   hypre_SMGRelaxData   *relax_data = relax_vdata;

   int                   num_spaces    = (relax_data -> num_spaces);
   int                  *space_indices = (relax_data -> space_indices);
   int                  *space_strides = (relax_data -> space_strides);
   hypre_StructVector   *temp_vec      = (relax_data -> temp_vec);

   hypre_StructStencil  *stencil       = hypre_StructMatrixStencil(A);     
   hypre_Index          *stencil_shape = hypre_StructStencilShape(stencil);
   int                   stencil_size  = hypre_StructStencilSize(stencil); 
   int                   stencil_dim   = hypre_StructStencilDim(stencil);
                       
   hypre_StructMatrix   *A_rem;
   void                **residual_data;

   hypre_Index           base_index;
   hypre_Index           base_stride;

   int                   num_stencil_indices;
   int                  *stencil_indices;
                       
   int                   i;

   int                   ierr = 0;

   /*----------------------------------------------------------
    * Free up old data before putting new data into structure
    *----------------------------------------------------------*/

   hypre_SMGRelaxDestroyARem(relax_vdata);

   /*----------------------------------------------------------
    * Set up data
    *----------------------------------------------------------*/

   hypre_CopyIndex((relax_data -> base_index),  base_index);
   hypre_CopyIndex((relax_data -> base_stride), base_stride);

   stencil_indices = hypre_TAlloc(int, stencil_size);
   num_stencil_indices = 0;
   for (i = 0; i < stencil_size; i++)
   {
      if (hypre_IndexD(stencil_shape[i], (stencil_dim - 1)) != 0)
      {
         stencil_indices[num_stencil_indices] = i;
         num_stencil_indices++;
      }
   }
   A_rem = hypre_StructMatrixCreateMask(A, num_stencil_indices, stencil_indices);
   hypre_TFree(stencil_indices);

   /* Set up residual_data */
   residual_data = hypre_TAlloc(void *, num_spaces);

   for (i = 0; i < num_spaces; i++)
   {
      hypre_IndexD(base_index,  (stencil_dim - 1)) = space_indices[i];
      hypre_IndexD(base_stride, (stencil_dim - 1)) = space_strides[i];

      residual_data[i] = hypre_SMGResidualCreate();
      hypre_SMGResidualSetBase(residual_data[i], base_index, base_stride);
      hypre_SMGResidualSetup(residual_data[i], A_rem, x, b, temp_vec);
   }

   (relax_data -> A_rem)         = A_rem;
   (relax_data -> residual_data) = residual_data;

   (relax_data -> setup_a_rem) = 0;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetupASol
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetupASol( void               *relax_vdata,
                         hypre_StructMatrix *A,
                         hypre_StructVector *b,
                         hypre_StructVector *x           )
{
   hypre_SMGRelaxData   *relax_data = relax_vdata;

   int                   num_spaces    = (relax_data -> num_spaces);
   int                  *space_indices = (relax_data -> space_indices);
   int                  *space_strides = (relax_data -> space_strides);
   hypre_StructVector   *temp_vec      = (relax_data -> temp_vec);

   int                   num_pre_relax   = (relax_data -> num_pre_relax);
   int                   num_post_relax  = (relax_data -> num_post_relax);

   hypre_StructStencil  *stencil       = hypre_StructMatrixStencil(A);     
   hypre_Index          *stencil_shape = hypre_StructStencilShape(stencil);
   int                   stencil_size  = hypre_StructStencilSize(stencil); 
   int                   stencil_dim   = hypre_StructStencilDim(stencil);
                       
   hypre_StructMatrix   *A_sol;
   void                **solve_data;

   hypre_Index           base_index;
   hypre_Index           base_stride;

   int                   num_stencil_indices;
   int                  *stencil_indices;
                       
   int                   i;

   int                   ierr = 0;

   /*----------------------------------------------------------
    * Free up old data before putting new data into structure
    *----------------------------------------------------------*/

   hypre_SMGRelaxDestroyASol(relax_vdata);

   /*----------------------------------------------------------
    * Set up data
    *----------------------------------------------------------*/

   hypre_CopyIndex((relax_data -> base_index),  base_index);
   hypre_CopyIndex((relax_data -> base_stride), base_stride);

   stencil_indices = hypre_TAlloc(int, stencil_size);
   num_stencil_indices = 0;
   for (i = 0; i < stencil_size; i++)
   {
      if (hypre_IndexD(stencil_shape[i], (stencil_dim - 1)) == 0)
      {
         stencil_indices[num_stencil_indices] = i;
         num_stencil_indices++;
      }
   }
   A_sol = hypre_StructMatrixCreateMask(A, num_stencil_indices, stencil_indices);
   hypre_StructStencilDim(hypre_StructMatrixStencil(A_sol)) = stencil_dim - 1;
   hypre_TFree(stencil_indices);

   /* Set up solve_data */
   solve_data    = hypre_TAlloc(void *, num_spaces);

   for (i = 0; i < num_spaces; i++)
   {
      hypre_IndexD(base_index,  (stencil_dim - 1)) = space_indices[i];
      hypre_IndexD(base_stride, (stencil_dim - 1)) = space_strides[i];

      if (stencil_dim > 2)
      {
         solve_data[i] = hypre_SMGCreate(relax_data -> comm);
         hypre_SMGSetNumPreRelax( solve_data[i], num_pre_relax);
         hypre_SMGSetNumPostRelax( solve_data[i], num_post_relax);
         hypre_SMGSetBase(solve_data[i], base_index, base_stride);
         hypre_SMGSetMemoryUse(solve_data[i], (relax_data -> memory_use));
         hypre_SMGSetTol(solve_data[i], 0.0);
         hypre_SMGSetMaxIter(solve_data[i], 1);
         hypre_SMGSetup(solve_data[i], A_sol, temp_vec, x);
      }
      else
      {
         solve_data[i] = hypre_CyclicReductionCreate(relax_data -> comm);
         hypre_CyclicReductionSetBase(solve_data[i], base_index, base_stride);
         hypre_CyclicReductionSetup(solve_data[i], A_sol, temp_vec, x);
      }
   }

   (relax_data -> A_sol)      = A_sol;
   (relax_data -> solve_data) = solve_data;

   (relax_data -> setup_a_sol) = 0;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetTempVec
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetTempVec( void               *relax_vdata,
                          hypre_StructVector *temp_vec    )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 ierr = 0;

   hypre_SMGRelaxDestroyTempVec(relax_vdata);
   (relax_data -> temp_vec) = hypre_StructVectorRef(temp_vec);

   (relax_data -> setup_temp_vec) = 1;
   (relax_data -> setup_a_rem)    = 1;
   (relax_data -> setup_a_sol)    = 1;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetMemoryUse
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetMemoryUse( void *relax_vdata,
                            int   memory_use  )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int               ierr = 0;

   (relax_data -> memory_use) = memory_use;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetTol
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetTol( void   *relax_vdata,
                      double  tol         )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 ierr = 0;

   (relax_data -> tol) = tol;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetMaxIter
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetMaxIter( void *relax_vdata,
                          int   max_iter    )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 ierr = 0;

   (relax_data -> max_iter) = max_iter;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetZeroGuess
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetZeroGuess( void *relax_vdata,
                            int   zero_guess  )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 ierr = 0;

   (relax_data -> zero_guess) = zero_guess;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetNumSpaces
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetNumSpaces( void *relax_vdata,
                            int   num_spaces      )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 i;
   int                 ierr = 0;

   (relax_data -> num_spaces) = num_spaces;

   hypre_TFree(relax_data -> space_indices);
   hypre_TFree(relax_data -> space_strides);
   hypre_TFree(relax_data -> pre_space_ranks);
   hypre_TFree(relax_data -> reg_space_ranks);
   (relax_data -> space_indices)   = hypre_TAlloc(int, num_spaces);
   (relax_data -> space_strides)   = hypre_TAlloc(int, num_spaces);
   (relax_data -> num_pre_spaces)  = 0;
   (relax_data -> num_reg_spaces)  = num_spaces;
   (relax_data -> pre_space_ranks) = NULL;
   (relax_data -> reg_space_ranks) = hypre_TAlloc(int, num_spaces);

   for (i = 0; i < num_spaces; i++)
   {
      (relax_data -> space_indices[i]) = 0;
      (relax_data -> space_strides[i]) = 1;
      (relax_data -> reg_space_ranks[i]) = i;
   }

   (relax_data -> setup_temp_vec) = 1;
   (relax_data -> setup_a_rem)    = 1;
   (relax_data -> setup_a_sol)    = 1;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetNumPreSpaces
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetNumPreSpaces( void *relax_vdata,
                               int   num_pre_spaces )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 i;
   int                 ierr = 0;

   (relax_data -> num_pre_spaces) = num_pre_spaces;

   hypre_TFree(relax_data -> pre_space_ranks);
   (relax_data -> pre_space_ranks) = hypre_TAlloc(int, num_pre_spaces);

   for (i = 0; i < num_pre_spaces; i++)
      (relax_data -> pre_space_ranks[i]) = 0;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetNumRegSpaces
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetNumRegSpaces( void *relax_vdata,
                               int   num_reg_spaces )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 i;
   int                 ierr = 0;

   (relax_data -> num_reg_spaces) = num_reg_spaces;

   hypre_TFree(relax_data -> reg_space_ranks);
   (relax_data -> reg_space_ranks) = hypre_TAlloc(int, num_reg_spaces);

   for (i = 0; i < num_reg_spaces; i++)
      (relax_data -> reg_space_ranks[i]) = 0;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetSpace
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetSpace( void *relax_vdata,
                        int   i,
                        int   space_index,
                        int   space_stride )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 ierr = 0;

   (relax_data -> space_indices[i]) = space_index;
   (relax_data -> space_strides[i]) = space_stride;

   (relax_data -> setup_temp_vec) = 1;
   (relax_data -> setup_a_rem)    = 1;
   (relax_data -> setup_a_sol)    = 1;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetRegSpaceRank
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetRegSpaceRank( void *relax_vdata,
                               int   i,
                               int   reg_space_rank )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 ierr = 0;

   (relax_data -> reg_space_ranks[i]) = reg_space_rank;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetPreSpaceRank
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetPreSpaceRank( void *relax_vdata,
                               int   i,
                               int   pre_space_rank  )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 ierr = 0;

   (relax_data -> pre_space_ranks[i]) = pre_space_rank;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetBase
 *--------------------------------------------------------------------------*/
 
int
hypre_SMGRelaxSetBase( void        *relax_vdata,
                       hypre_Index  base_index,
                       hypre_Index  base_stride )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 d;
   int                 ierr = 0;
 
   for (d = 0; d < 3; d++)
   {
      hypre_IndexD((relax_data -> base_index),  d) =
         hypre_IndexD(base_index,  d);
      hypre_IndexD((relax_data -> base_stride), d) =
         hypre_IndexD(base_stride, d);
   }
 
   if ((relax_data -> base_box_array) != NULL)
   {
      hypre_BoxArrayDestroy((relax_data -> base_box_array));
      (relax_data -> base_box_array) = NULL;
   }

   (relax_data -> setup_temp_vec) = 1;
   (relax_data -> setup_a_rem)    = 1;
   (relax_data -> setup_a_sol)    = 1;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetNumPreRelax
 * Note that we require at least 1 pre-relax sweep.
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetNumPreRelax( void *relax_vdata,
                              int   num_pre_relax )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 ierr = 0;

   (relax_data -> num_pre_relax) = hypre_max(num_pre_relax,1);

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetNumPostRelax
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetNumPostRelax( void *relax_vdata,
                               int   num_post_relax )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;
   int                 ierr = 0;

   (relax_data -> num_post_relax) = num_post_relax;

   return ierr;
}

/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetNewMatrixStencil
 *--------------------------------------------------------------------------*/
 
int
hypre_SMGRelaxSetNewMatrixStencil( void                *relax_vdata,
                                   hypre_StructStencil *diff_stencil )
{
   hypre_SMGRelaxData *relax_data = relax_vdata;

   hypre_Index        *stencil_shape = hypre_StructStencilShape(diff_stencil);
   int                 stencil_size  = hypre_StructStencilSize(diff_stencil); 
   int                 stencil_dim   = hypre_StructStencilDim(diff_stencil);
                         
   int                 i;
                     
   int                 ierr = 0;

   for (i = 0; i < stencil_size; i++)
   {
      if (hypre_IndexD(stencil_shape[i], (stencil_dim - 1)) != 0)
      {
         (relax_data -> setup_a_rem) = 1;
      }
      else
      {
         (relax_data -> setup_a_sol) = 1;
      }
   }

   return ierr;
}


/*--------------------------------------------------------------------------
 * hypre_SMGRelaxSetupBaseBoxArray
 *--------------------------------------------------------------------------*/

int
hypre_SMGRelaxSetupBaseBoxArray( void               *relax_vdata,
                                 hypre_StructMatrix *A,
                                 hypre_StructVector *b,
                                 hypre_StructVector *x           )
{
   hypre_SMGRelaxData  *relax_data = relax_vdata;
                       
   hypre_StructGrid    *grid;
   hypre_BoxArray      *boxes;
   hypre_BoxArray      *base_box_array;
                       
   int                  ierr = 0;

   grid  = hypre_StructVectorGrid(x);
   boxes = hypre_StructGridBoxes(grid);

   base_box_array = hypre_BoxArrayDuplicate(boxes);
   hypre_ProjectBoxArray(base_box_array, 
                         (relax_data -> base_index),
                         (relax_data -> base_stride));

   (relax_data -> base_box_array) = base_box_array;

   return ierr;
}