/******************************************************************************
 * Copyright (c) 2024 Lawrence Livermore National Security, LLC and other
 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
 *
 * SPDX-License-Identifier: MIT
 ******************************************************************************/

#include "internal/amg.h"
#include <stdlib.h>
#include <string.h>
#include "HYPRE_parcsr_mv.h"
#include "_hypre_IJ_mv.h"     // For hypre_IJVectorGlobalNumRows
#include "_hypre_parcsr_mv.h" // For hypre_ParVectorComm, hypre_ParVectorInitialize_v2
#include "_hypre_utilities.h" // For hypre_TAlloc, hypre_TMemcpy, hypre_TFree
#include "internal/gen_macros.h"
#include "logging.h"

/*-----------------------------------------------------------------------------
 * Define Field/Offset/Setter mappings
 *-----------------------------------------------------------------------------*/

/* AMG's interpolation fields */
#define AMGint_FIELDS(_prefix)                                                      \
   ADD_FIELD_OFFSET_ENTRY(_prefix, prolongation_type, hypredrv_FieldTypeIntSet)     \
   ADD_FIELD_OFFSET_ENTRY(_prefix, restriction_type, hypredrv_FieldTypeIntSet)      \
   ADD_FIELD_OFFSET_ENTRY(_prefix, max_nnz_row, hypredrv_FieldTypeIntSet)           \
   ADD_FIELD_OFFSET_ENTRY(_prefix, trunc_factor, hypredrv_FieldTypeDoubleSet)       \
   ADD_FIELD_OFFSET_ENTRY(_prefix, restrict_strong_th, hypredrv_FieldTypeDoubleSet) \
   ADD_FIELD_OFFSET_ENTRY(_prefix, restrict_filter_th, hypredrv_FieldTypeDoubleSet)

/* AMG's coarsening fields */
#define AMGcsn_FIELDS(_prefix)                                                          \
   ADD_FIELD_OFFSET_ENTRY(_prefix, type, hypredrv_FieldTypeIntSet)                      \
   ADD_FIELD_OFFSET_ENTRY(_prefix, rap2, hypredrv_FieldTypeIntSet)                      \
   ADD_FIELD_OFFSET_ENTRY(_prefix, mod_rap2, hypredrv_FieldTypeIntSet)                  \
   ADD_FIELD_OFFSET_ENTRY(_prefix, keep_transpose, hypredrv_FieldTypeIntSet)            \
   ADD_FIELD_OFFSET_ENTRY(_prefix, sabs, hypredrv_FieldTypeIntSet)                      \
   ADD_FIELD_OFFSET_ENTRY(_prefix, num_functions, hypredrv_FieldTypeIntSet)             \
   ADD_FIELD_OFFSET_ENTRY(_prefix, nodal_type, hypredrv_FieldTypeIntSet)                \
   ADD_FIELD_OFFSET_ENTRY(_prefix, interp_vec_variant, hypredrv_FieldTypeIntSet)        \
   ADD_FIELD_OFFSET_ENTRY(_prefix, smooth_interp_vecs, hypredrv_FieldTypeIntSet)        \
   ADD_FIELD_OFFSET_ENTRY(_prefix, interp_vec_qmax, hypredrv_FieldTypeIntSet)           \
   ADD_FIELD_OFFSET_ENTRY(_prefix, interp_vec_abs_q_trunc, hypredrv_FieldTypeDoubleSet) \
   ADD_FIELD_OFFSET_ENTRY(_prefix, filter_functions, hypredrv_FieldTypeIntSet)          \
   ADD_FIELD_OFFSET_ENTRY(_prefix, nodal, hypredrv_FieldTypeIntSet)                     \
   ADD_FIELD_OFFSET_ENTRY(_prefix, seq_amg_th, hypredrv_FieldTypeIntSet)                \
   ADD_FIELD_OFFSET_ENTRY(_prefix, min_coarse_size, hypredrv_FieldTypeIntSet)           \
   ADD_FIELD_OFFSET_ENTRY(_prefix, max_coarse_size, hypredrv_FieldTypeIntSet)           \
   ADD_FIELD_OFFSET_ENTRY(_prefix, max_levels, hypredrv_FieldTypeIntSet)                \
   ADD_FIELD_OFFSET_ENTRY(_prefix, max_row_sum, hypredrv_FieldTypeDoubleSet)            \
   ADD_FIELD_OFFSET_ENTRY(_prefix, strong_th, hypredrv_FieldTypeDoubleSet)

/* AMG's aggressive coarsening fields */
#define AMGagg_FIELDS(_prefix)                                                    \
   ADD_FIELD_OFFSET_ENTRY(_prefix, num_levels, hypredrv_FieldTypeIntSet)          \
   ADD_FIELD_OFFSET_ENTRY(_prefix, num_paths, hypredrv_FieldTypeIntSet)           \
   ADD_FIELD_OFFSET_ENTRY(_prefix, prolongation_type, hypredrv_FieldTypeIntSet)   \
   ADD_FIELD_OFFSET_ENTRY(_prefix, max_nnz_row, hypredrv_FieldTypeIntSet)         \
   ADD_FIELD_OFFSET_ENTRY(_prefix, trunc_factor, hypredrv_FieldTypeDoubleSet)     \
   ADD_FIELD_OFFSET_ENTRY(_prefix, P12_max_elements, hypredrv_FieldTypeDoubleSet) \
   ADD_FIELD_OFFSET_ENTRY(_prefix, P12_trunc_factor, hypredrv_FieldTypeDoubleSet)

/* AMG's relaxation fields */
#define AMGrlx_FIELDS(_prefix)                                                \
   ADD_FIELD_OFFSET_ENTRY(_prefix, down_type, hypredrv_FieldTypeIntSet)       \
   ADD_FIELD_OFFSET_ENTRY(_prefix, up_type, hypredrv_FieldTypeIntSet)         \
   ADD_FIELD_OFFSET_ENTRY(_prefix, coarse_type, hypredrv_FieldTypeIntSet)     \
   ADD_FIELD_OFFSET_ENTRY(_prefix, down_sweeps, hypredrv_FieldTypeIntSet)     \
   ADD_FIELD_OFFSET_ENTRY(_prefix, up_sweeps, hypredrv_FieldTypeIntSet)       \
   ADD_FIELD_OFFSET_ENTRY(_prefix, coarse_sweeps, hypredrv_FieldTypeIntSet)   \
   ADD_FIELD_OFFSET_ENTRY(_prefix, num_sweeps, hypredrv_FieldTypeIntSet)      \
   ADD_FIELD_OFFSET_ENTRY(_prefix, order, hypredrv_FieldTypeIntSet)           \
   ADD_FIELD_OFFSET_ENTRY(_prefix, points, hypredrv_FieldTypeIntSet)          \
   ADD_FIELD_OFFSET_ENTRY(_prefix, weight, hypredrv_FieldTypeDoubleSet)       \
   ADD_FIELD_OFFSET_ENTRY(_prefix, outer_weight, hypredrv_FieldTypeDoubleSet) \
   ADD_FIELD_OFFSET_ENTRY(_prefix, chebyshev, hypredrv_ChebySetArgs)

/* AMG's complex smoother fields */
#define AMGsmt_FIELDS(_prefix)                                           \
   ADD_FIELD_OFFSET_ENTRY(_prefix, type, hypredrv_FieldTypeIntSet)       \
   ADD_FIELD_OFFSET_ENTRY(_prefix, num_levels, hypredrv_FieldTypeIntSet) \
   ADD_FIELD_OFFSET_ENTRY(_prefix, num_sweeps, hypredrv_FieldTypeIntSet) \
   ADD_FIELD_OFFSET_ENTRY(_prefix, fsai, hypredrv_FSAISetArgs)           \
   ADD_FIELD_OFFSET_ENTRY(_prefix, ilu, hypredrv_ILUSetArgs)

/* AMG */
#define AMG_FIELDS(_prefix)                                                \
   ADD_FIELD_OFFSET_ENTRY(_prefix, max_iter, hypredrv_FieldTypeIntSet)     \
   ADD_FIELD_OFFSET_ENTRY(_prefix, print_level, hypredrv_FieldTypeIntSet)  \
   ADD_FIELD_OFFSET_ENTRY(_prefix, tolerance, hypredrv_FieldTypeDoubleSet) \
   ADD_FIELD_OFFSET_ENTRY(_prefix, interpolation, hypredrv_AMGintSetArgs)  \
   ADD_FIELD_OFFSET_ENTRY(_prefix, aggressive, hypredrv_AMGaggSetArgs)     \
   ADD_FIELD_OFFSET_ENTRY(_prefix, coarsening, hypredrv_AMGcsnSetArgs)     \
   ADD_FIELD_OFFSET_ENTRY(_prefix, relaxation, hypredrv_AMGrlxSetArgs)     \
   ADD_FIELD_OFFSET_ENTRY(_prefix, smoother, hypredrv_AMGsmtSetArgs)

/* Define the prefix list */
#define GENERATE_PREFIXED_LIST_AMG      \
   GENERATE_PREFIXED_COMPONENTS(AMGint) \
   GENERATE_PREFIXED_COMPONENTS(AMGcsn) \
   GENERATE_PREFIXED_COMPONENTS(AMGagg) \
   GENERATE_PREFIXED_COMPONENTS(AMGrlx) \
   GENERATE_PREFIXED_COMPONENTS(AMGsmt) \
   GENERATE_PREFIXED_COMPONENTS(AMG)

/* Define num_fields macros for each struct prefix */
// clang-format off
#define AMGint_NUM_FIELDS (sizeof(AMGint_field_offset_map) / sizeof(AMGint_field_offset_map[0]))
#define AMGcsn_NUM_FIELDS (sizeof(AMGcsn_field_offset_map) / sizeof(AMGcsn_field_offset_map[0]))
#define AMGagg_NUM_FIELDS (sizeof(AMGagg_field_offset_map) / sizeof(AMGagg_field_offset_map[0]))
#define AMGrlx_NUM_FIELDS (sizeof(AMGrlx_field_offset_map) / sizeof(AMGrlx_field_offset_map[0]))
#define AMGsmt_NUM_FIELDS (sizeof(AMGsmt_field_offset_map) / sizeof(AMGsmt_field_offset_map[0]))
#define AMG_NUM_FIELDS    (sizeof(AMG_field_offset_map)    / sizeof(AMG_field_offset_map[0]))
// clang-format on

/* Iterates over each prefix in the list and
   generates the various function declarations/definitions and field_offset_map object */
GENERATE_PREFIXED_LIST_AMG                               // LCOV_EXCL_LINE
hypredrv_DEFINE_VOID_GET_VALID_VALUES_FUNC(hypredrv_AMG) // LCOV_EXCL_LINE

   /*-----------------------------------------------------------------------------
    * hypredrv_AMGintSetDefaultArgs
    *-----------------------------------------------------------------------------*/

   void hypredrv_AMGintSetDefaultArgs(AMGint_args *args)
{
   args->prolongation_type  = 6;
   args->restriction_type   = 0;
   args->max_nnz_row        = 4;
   args->trunc_factor       = 0.0;
   args->restrict_strong_th = 0.25;
   args->restrict_filter_th = 0.0;
}

/*-----------------------------------------------------------------------------
 * AMGcsnSetDefaultArgs
 *-----------------------------------------------------------------------------*/

void
hypredrv_AMGcsnSetDefaultArgs(AMGcsn_args *args)
{
   args->rap2 = 0;
#ifdef HYPRE_USING_GPU
   args->mod_rap2       = 1;
   args->keep_transpose = 1;
   args->type           = 8;
#else
   args->mod_rap2       = 0;
   args->keep_transpose = 0;
   args->type           = 10;
#endif
   args->num_functions = 1;
   /* GM/LN knobs: "auto" (-1) reproduces the historical hard-coded setup, i.e.
      GM2 with QMax 4 and smoothed interpolation vectors for a full-system near
      null space, GM1 for one projected onto an MGR F-block. */
   args->nodal_type             = 4;
   args->interp_vec_variant     = -1;
   args->smooth_interp_vecs     = -1;
   args->interp_vec_qmax        = -1;
   args->interp_vec_abs_q_trunc = 0.0;
   args->sabs                   = 0;
   args->filter_functions       = 0;
   args->nodal                  = 0;
   args->seq_amg_th             = 0;
   args->min_coarse_size        = 0;
   args->max_coarse_size        = 64;
   args->max_levels             = 25;
   args->max_row_sum            = 0.9;
   args->strong_th              = 0.25;
}

/*-----------------------------------------------------------------------------
 * AMGaggSetDefaultArgs
 *-----------------------------------------------------------------------------*/

void
hypredrv_AMGaggSetDefaultArgs(AMGagg_args *args)
{
   args->num_levels        = 0;
   args->num_paths         = 1;
   args->prolongation_type = 4;
   args->max_nnz_row       = 0;
   args->P12_max_elements  = 0;
   args->P12_trunc_factor  = 0.0;
   args->trunc_factor      = 0.0;
}

/*-----------------------------------------------------------------------------
 * AMGrlxSetDefaultArgs
 *-----------------------------------------------------------------------------*/

void
hypredrv_AMGrlxSetDefaultArgs(AMGrlx_args *args)
{
   args->type = -1;
#ifdef HYPRE_USING_GPU
   args->down_type = 18;
   args->up_type   = 18;
#else
   args->down_type = 13;
   args->up_type   = 14;
#endif
   args->coarse_type   = 9;
   args->down_sweeps   = -1;
   args->up_sweeps     = -1;
   args->coarse_sweeps = 1;
   args->num_sweeps    = 1;
   args->order         = 0;
   args->points        = 0;
   args->weight        = 1.0;
   args->outer_weight  = 1.0;

   hypredrv_ChebySetDefaultArgs(&args->chebyshev);
}

/*-----------------------------------------------------------------------------
 * hypredrv_AMGsmtSetDefaultArgs
 *-----------------------------------------------------------------------------*/

void
hypredrv_AMGsmtSetDefaultArgs(AMGsmt_args *args)
{
   args->type       = 5;
   args->num_levels = 0;
   args->num_sweeps = 1;
   hypredrv_FSAISetDefaultArgs(&args->fsai);
   hypredrv_ILUSetDefaultArgs(&args->ilu);
}

/*-----------------------------------------------------------------------------
 * AMGSetDefaultArgs
 *-----------------------------------------------------------------------------*/

void
hypredrv_AMGSetDefaultArgs(AMG_args *args)
{
   args->max_iter               = 1;
   args->print_level            = 0;
   args->tolerance              = 0.0;
   args->interp_vec_variant     = 2;
   args->num_rbms               = 0;
   args->rbms[0]                = NULL;
   args->rbms[1]                = NULL;
   args->rbms[2]                = NULL;
   args->rbm_source_generation  = 0;
   args->rbm_source_labels_size = 0;
   args->rbm_source_labels_hash = 0;

   hypredrv_AMGintSetDefaultArgs(&args->interpolation);
   hypredrv_AMGaggSetDefaultArgs(&args->aggressive);
   hypredrv_AMGcsnSetDefaultArgs(&args->coarsening);
   hypredrv_AMGrlxSetDefaultArgs(&args->relaxation);
   hypredrv_AMGsmtSetDefaultArgs(&args->smoother);
}

/*-----------------------------------------------------------------------------
 * AMGintGetValidValues
 *-----------------------------------------------------------------------------*/

// clang-format off
StrIntMapArray
hypredrv_AMGintGetValidValues(const char *key)
{
   if (!strcmp(key, "prolongation_type"))
   {
      static StrIntMap map[] = {{"mod_classical",         0},
                                {"least_squares",         1},
                                {"mod_classical_he",      2},
                                {"direct_sep_weights",    3},
                                {"multipass",             4},
                                {"multipass_sep_weights", 5},
                                {"extended+i",            6},
                                {"extended+i_c",          7},
                                {"standard",              8},
                                {"standard_sep_weights",  9},
                                {"blk_classical",        10},
                                {"blk_classical_diag",   11},
                                {"f_f",                  12},
                                {"f_f1",                 13},
                                {"extended",             14},
                                {"direct_sep_weights",   15},
                                {"mm_extended",          16},
                                {"mm_extended+i",        17},
                                {"mm-ext+i",             17},
                                {"mm_extended+e",        18},
                                {"mm-ext+e",             18},
                                {"blk_direct",           24},
                                {"one_point",           100},};

      return STR_INT_MAP_ARRAY_CREATE(map);
   }
   if (!strcmp(key, "restriction_type"))
   {
      static StrIntMap map[] = {{"p_transpose",    0},
                                {"air_1",          1},
                                {"air_2",          2},
                                {"neumann_air_0",  3},
                                {"neumann_air_1",  4},
                                {"neumann_air_2",  5},
                                {"air_1.5",       15},};

      return STR_INT_MAP_ARRAY_CREATE(map);
   }
   else
   {
      return STR_INT_MAP_ARRAY_VOID();
   }
}

/*-----------------------------------------------------------------------------
 * AMGcsnGetValidValues
 *-----------------------------------------------------------------------------*/

StrIntMapArray
hypredrv_AMGcsnGetValidValues(const char *key)
{
   if (!strcmp(key, "type"))
   {
      static StrIntMap map[] = {{"cljp",    0},
                                {"rs",      1},
                                {"rs3",     3},
                                {"falgout", 6},
                                {"pmis",    8},
                                {"hmis",   10},};

      return STR_INT_MAP_ARRAY_CREATE(map);
   }
   if (!strcmp(key, "filter_functions") ||
       !strcmp(key, "nodal") ||
       !strcmp(key, "rap2") ||
       !strcmp(key, "mod_rap2") ||
       !strcmp(key, "keep_transpose"))
   {
      return STR_INT_MAP_ARRAY_CREATE_ON_OFF();
   }
   else
   {
      return STR_INT_MAP_ARRAY_VOID();
   }
}

/*-----------------------------------------------------------------------------
 * AMGaggGetValidValues
 *-----------------------------------------------------------------------------*/

StrIntMapArray
hypredrv_AMGaggGetValidValues(const char *key)
{
   if (!strcmp(key, "prolongation_type"))
   {
      static StrIntMap map[] = {{"2_stage_extended+i", 1},
                                {"2_stage_standard",   2},
                                {"2_stage_extended",   3},
                                {"multipass",          4},
                                {"mm_extended",        5},
                                {"mm_extended+i",      6},
                                {"mm_extended+e",      7},};

      return STR_INT_MAP_ARRAY_CREATE(map);
   }
   else
   {
      return STR_INT_MAP_ARRAY_VOID();
   }
}

/*-----------------------------------------------------------------------------
 * AMGrlxGetValidValues
 *-----------------------------------------------------------------------------*/

StrIntMapArray
hypredrv_AMGrlxGetValidValues(const char *key)
{
   if (!strcmp(key, "down_type"))
   {
      static StrIntMap map[] = {{"jacobi_non_mv",  0},
                                {"forward-hgs",    3},
                                {"backward-hgs",   4},
                                {"chaotic-hgs",    5},
                                {"hsgs",           6},
                                {"jacobi",         7},
                                {"l1-hsgs",        8},
                                {"forward-solve", 10},
                                {"2gs-it1",       11},
                                {"2gs-it2",       12},
                                {"forward-hl1gs", 13},
                                {"backward-hl1gs", 14},
                                {"cg",            15},
                                {"chebyshev",     16},
                                {"l1-jacobi",     18},
                                {"l1sym-hgs",     89},};

      return STR_INT_MAP_ARRAY_CREATE(map);
   }
   if (!strcmp(key, "up_type"))
   {
      static StrIntMap map[] = {{"jacobi_non_mv",   0},
                                {"forward-hgs",     3},
                                {"backward-hgs",    4},
                                {"chaotic-hgs",     5},
                                {"hsgs",            6},
                                {"jacobi",          7},
                                {"l1-hsgs",         8},
                                {"forward-solve",  10},
                                {"2gs-it1",        11},
                                {"2gs-it2",        12},
                                {"forward-hl1gs",  13},
                                {"backward-hl1gs", 14},
                                {"cg",             15},
                                {"chebyshev",      16},
                                {"l1-jacobi",      18},
                                {"l1sym-hgs",      89},};

      return STR_INT_MAP_ARRAY_CREATE(map);
   }
   if (!strcmp(key, "points"))
   {
      /* Which grid points each sweep of the cycle relaxes. "all" leaves hypre's
         default schedule untouched; "air" applies the F/C schedule that the
         approximate ideal restriction algorithm expects. */
      static StrIntMap map[] = {{"all", 0},
                                {"air", 1},};

      return STR_INT_MAP_ARRAY_CREATE(map);
   }
   if (!strcmp(key, "coarse_type"))
   {
      static StrIntMap map[] = {{"jacobi_non_mv",   0},
                                {"hsgs",            6},
                                {"jacobi",          7},
                                {"l1-hsgs",         8},
                                {"ge",              9},
                                {"2gs-it1",        11},
                                {"2gs-it2",        12},
                                {"forward-hl1gs",  13},
                                {"backward-hl1gs", 14},
                                {"cg",             15},
                                {"chebyshev",      16},
                                {"l1-jacobi",      18},
                                {"l1sym-hgs",      89},
                                {"lu_piv",         99},
                                {"lu_inv",        199},};

      return STR_INT_MAP_ARRAY_CREATE(map);
   }
   else
   {
      return STR_INT_MAP_ARRAY_VOID();
   }
}

/*-----------------------------------------------------------------------------
 * AMGsmtGetValidValues
 *-----------------------------------------------------------------------------*/

StrIntMapArray
hypredrv_AMGsmtGetValidValues(const char *key)
{
   if (!strcmp(key, "type"))
   {
      static StrIntMap map[] = {{"fsai",      4},
                                {"ilu",       5},
                                {"schwarz",   6},
                                {"pilut",     7},
                                {"parasails", 8},
                                {"euclid",    9},};

      return STR_INT_MAP_ARRAY_CREATE(map);
   }
   else
   {
      return STR_INT_MAP_ARRAY_VOID();
   }
}
// clang-format on

/*-----------------------------------------------------------------------------
 * AMGSetRBMs
 *-----------------------------------------------------------------------------*/

void
hypredrv_AMGDestroyRBMs(AMG_args *args)
{
   if (!args)
   {
      return;
   }

   for (HYPRE_Int i = 0; i < args->num_rbms; i++)
   {
      HYPRE_ParVectorDestroy(args->rbms[i]);
      args->rbms[i] = NULL;
   }
   args->interp_vec_variant     = 2;
   args->num_rbms               = 0;
   args->rbm_source_generation  = 0;
   args->rbm_source_labels_size = 0;
   args->rbm_source_labels_hash = 0;
}

#if HYPRE_CHECK_MIN_VERSION(22600, 0)
enum
{
   AMG_NUM_RBMS            = 3,
   AMG_FIRST_RBM_COMPONENT = 3,
};

static int
AMGIntCompare(const void *lhs, const void *rhs)
{
   int left  = *(const int *)lhs;
   int right = *(const int *)rhs;
   return (left > right) - (left < right);
}

/* Copy the complete multivector to host memory. HYPRE_IJVectorGetValues()
 * expects its output buffer in the vector's execution memory on GPU builds,
 * so passing a malloc'ed buffer to it is not portable. Keeping the original
 * strides also handles both supported multivector storage layouts. */
static HYPRE_Complex *
AMGNearNullDataToHost(HYPRE_IJVector vec_nn, HYPRE_Int num_entries,
                      HYPRE_Int *vector_stride, HYPRE_Int *index_stride)
{
   void            *object     = NULL;
   hypre_ParVector *par_vector = NULL;
   hypre_Vector    *local      = NULL;

   HYPRE_IJVectorGetObject(vec_nn, &object);
   par_vector = (hypre_ParVector *)object;
   if (!par_vector || !(local = hypre_ParVectorLocalVector(par_vector)) ||
       hypre_VectorSize(local) != num_entries ||
       hypre_VectorNumVectors(local) < AMG_FIRST_RBM_COMPONENT + AMG_NUM_RBMS)
   {
      return NULL;
   }

   size_t total_entries = (size_t)num_entries * (size_t)hypre_VectorNumVectors(local);
   HYPRE_Complex *host_values =
      hypre_TAlloc(HYPRE_Complex, total_entries, HYPRE_MEMORY_HOST);
   if (total_entries > 0 && !host_values)
   {
      return NULL;
   }

   if (total_entries > 0)
   {
      HYPRE_Complex *source = hypre_VectorData(local);
      if (!source)
      {
         hypre_TFree(host_values, HYPRE_MEMORY_HOST);
         return NULL;
      }
      hypre_TMemcpy(host_values, source, HYPRE_Complex, total_entries, HYPRE_MEMORY_HOST,
                    hypre_VectorMemoryLocation(local));
   }

   *vector_stride = hypre_VectorVectorStride(local);
   *index_stride  = hypre_VectorIndexStride(local);
   return host_values;
}

/* Build the three displacement RBMs from a host copy of the near-null-space
 * multivector. A NULL selection mask copies every local row. */
static int
AMGCreateRBMsFromHost(AMG_args *args, MPI_Comm comm, HYPRE_BigInt global_size,
                      HYPRE_BigInt local_start, HYPRE_Int local_size,
                      HYPRE_Int source_size, const unsigned char *selected,
                      const HYPRE_Complex *host_values, HYPRE_Int vector_stride,
                      HYPRE_Int index_stride)
{
   hypredrv_AMGDestroyRBMs(args);

   for (HYPRE_Int mode = 0; mode < AMG_NUM_RBMS; mode++)
   {
      /* hypre_ParVectorCreate copies these two entries. */
      HYPRE_BigInt partitioning[2] = {local_start, local_start + local_size};
      HYPRE_ParVectorCreate(comm, global_size, partitioning, &args->rbms[mode]);
      int local_ok = args->rbms[mode] != NULL;
      if (local_ok)
      {
         hypre_ParVectorInitialize_v2(args->rbms[mode], HYPRE_MEMORY_HOST);
         hypre_Vector *local = hypre_ParVectorLocalVector(args->rbms[mode]);
         local_ok            = local && (local_size == 0 || hypre_VectorData(local));
      }

      int global_ok = 0;
      MPI_Allreduce(&local_ok, &global_ok, 1, MPI_INT, MPI_MIN, comm);
      if (!global_ok)
      {
         if (args->rbms[mode])
         {
            HYPRE_ParVectorDestroy(args->rbms[mode]);
            args->rbms[mode] = NULL;
         }
         hypredrv_AMGDestroyRBMs(args);
         hypredrv_ErrorCodeSet(ERROR_ALLOCATION);
         hypredrv_ErrorMsgAdd("Failed to create rigid-body-mode vectors");
         return 0;
      }

      args->num_rbms++;
      HYPRE_Complex *data =
         hypre_VectorData(hypre_ParVectorLocalVector(args->rbms[mode]));
      HYPRE_Int projected_i = 0;
      for (HYPRE_Int i = 0; i < source_size; i++)
      {
         if (!selected || selected[i])
         {
            size_t source_index =
               ((size_t)(AMG_FIRST_RBM_COMPONENT + mode) * (size_t)vector_stride) +
               ((size_t)i * (size_t)index_stride);
            data[projected_i++] = host_values[source_index];
         }
      }
   }

   return 1;
}
#endif

void
hypredrv_AMGSetRBMs(AMG_args *args, HYPRE_IJVector vec_nn)
{
   HYPRE_BigInt jlower = 0, jupper = 0;
   HYPRE_Int    num_entries = 0;

   if (!args)
   {
      return;
   }
   args->interp_vec_variant = 2;
   hypredrv_AMGDestroyRBMs(args);

   /* Sanity: check if the near null space vector is set
      We do not error out when NOT using nodal coarsening. */
   if (!vec_nn || !args->coarsening.nodal)
   {
      if (args->coarsening.nodal)
      {
         hypredrv_ErrorCodeSet(ERROR_UNKNOWN);
         hypredrv_ErrorMsgAdd("Near null space vectors (RBMs) required"
                              " for nodal coarsening, but not set");
      }
      return;
   }

#if HYPRE_CHECK_MIN_VERSION(22600, 0)
   HYPRE_IJVectorGetLocalRange(vec_nn, &jlower, &jupper);
   num_entries                  = (HYPRE_Int)(jupper - jlower + 1);
   HYPRE_Int      vector_stride = 0, index_stride = 0;
   HYPRE_Complex *host_values =
      AMGNearNullDataToHost(vec_nn, num_entries, &vector_stride, &index_stride);
   int      local_ok  = num_entries == 0 || host_values != NULL;
   int      global_ok = 0;
   MPI_Comm comm      = hypre_IJVectorComm((hypre_IJVector *)vec_nn);
   MPI_Allreduce(&local_ok, &global_ok, 1, MPI_INT, MPI_MIN, comm);
   if (!global_ok)
   {
      hypredrv_ErrorCodeSet(ERROR_ALLOCATION);
      hypredrv_ErrorMsgAdd("Failed to access near-null-space components on host");
      hypre_TFree(host_values, HYPRE_MEMORY_HOST);
      return;
   }

   AMGCreateRBMsFromHost(args, comm, hypre_IJVectorGlobalNumRows(vec_nn), jlower,
                         num_entries, num_entries, NULL, host_values, vector_stride,
                         index_stride);
   hypre_TFree(host_values, HYPRE_MEMORY_HOST);
#else
   (void)vec_nn;
   return;
#endif
}

/*-----------------------------------------------------------------------------
 * hypredrv_AMGSetProjectedRBMs
 *
 * Project full-system rigid-body modes onto selected MGR labels while preserving
 * each rank's row order, which is also the ordering of hypre's extracted block.
 *-----------------------------------------------------------------------------*/

static uint64_t
AMGSelectedDofsHash(const int *selected_dofs, size_t num_selected_dofs)
{
   uint64_t hash = UINT64_C(1469598103934665603);
   for (size_t i = 0; i < num_selected_dofs; i++)
   {
      hash ^= (uint64_t)(uint32_t)selected_dofs[i];
      hash *= UINT64_C(1099511628211);
   }
   return hash;
}

void
hypredrv_AMGSetProjectedRBMs(AMG_args *args, HYPRE_IJVector vec_nn,
                             const IntArray *dofmap, uint64_t input_generation,
                             const int *selected_dofs, size_t num_selected_dofs)
{
#if HYPRE_CHECK_MIN_VERSION(22600, 0)
   HYPRE_BigInt   jlower = 0, jupper = -1;
   HYPRE_Int      num_entries      = 0;
   HYPRE_Int      num_components   = 0;
   unsigned char *selected         = NULL;
   uint64_t       projected_local  = 0;
   uint64_t       projected_global = 0;
   uint64_t       projected_scan   = 0;
   int           *selected_labels  = NULL;
   MPI_Comm       comm             = MPI_COMM_NULL;

   if (!args)
   {
      return;
   }
   if (num_selected_dofs && !selected_dofs)
   {
      hypredrv_ErrorCodeSet(ERROR_INVALID_VAL);
      return;
   }
   uint64_t labels_hash = AMGSelectedDofsHash(selected_dofs, num_selected_dofs);
   if (!args->coarsening.nodal)
   {
      hypredrv_AMGDestroyRBMs(args);
      args->interp_vec_variant = 1;
      return;
   }
   if (!vec_nn || !dofmap)
   {
      /* Nodal AMG is valid without user interpolation vectors. Passing zero
       * vectors below preserves hypre's standard nodal-coarsening behavior. */
      hypredrv_AMGDestroyRBMs(args);
      args->interp_vec_variant = 1;
      return;
   }

   comm                = hypre_IJVectorComm((hypre_IJVector *)vec_nn);
   int local_cache_hit = args->interp_vec_variant == 1 &&
                         args->rbm_source_generation == input_generation &&
                         args->rbm_source_labels_size == num_selected_dofs &&
                         args->rbm_source_labels_hash == labels_hash;
   int global_cache_hit = 0;
   MPI_Allreduce(&local_cache_hit, &global_cache_hit, 1, MPI_INT, MPI_MIN, comm);
   if (global_cache_hit)
   {
      return;
   }

   hypredrv_AMGDestroyRBMs(args);
   args->interp_vec_variant = 1;
   num_components           = hypre_IJVectorNumComponents(vec_nn);
   HYPRE_IJVectorGetLocalRange(vec_nn, &jlower, &jupper);
   num_entries     = (jupper >= jlower) ? (HYPRE_Int)(jupper - jlower + 1) : 0;
   int local_valid = num_components >= AMG_FIRST_RBM_COMPONENT + AMG_NUM_RBMS &&
                     (size_t)num_entries == dofmap->size;
   int global_valid = 0;
   MPI_Allreduce(&local_valid, &global_valid, 1, MPI_INT, MPI_MIN, comm);
   if (!global_valid)
   {
      hypredrv_ErrorCodeSet(ERROR_INVALID_VAL);
      hypredrv_ErrorMsgAdd(
         "Elasticity near null space requires 6 components and a local size matching "
         "the MGR dofmap on every rank");
      return;
   }

   selected = (unsigned char *)calloc((size_t)num_entries, sizeof(unsigned char));
   if (num_selected_dofs > 0)
   {
      selected_labels = (int *)malloc(num_selected_dofs * sizeof(*selected_labels));
   }
   HYPRE_Int      vector_stride = 0, index_stride = 0;
   HYPRE_Complex *host_values =
      AMGNearNullDataToHost(vec_nn, num_entries, &vector_stride, &index_stride);
   int local_ok = (num_selected_dofs == 0 || selected_labels) &&
                  (num_entries == 0 || (selected && host_values));
   int global_ok = 0;
   MPI_Allreduce(&local_ok, &global_ok, 1, MPI_INT, MPI_MIN, comm);
   if (!global_ok || (num_selected_dofs > 0 && !selected_labels))
   {
      hypredrv_ErrorCodeSet(ERROR_ALLOCATION);
      hypredrv_ErrorMsgAdd(
         "Failed to allocate projected near-null-space buffers (%d entries)",
         (int)num_entries);
      free(selected);
      free(selected_labels);
      hypre_TFree(host_values, HYPRE_MEMORY_HOST);
      return;
   }

   if (num_selected_dofs > 0)
   {
      memcpy(selected_labels, selected_dofs,
             num_selected_dofs * sizeof(*selected_labels));
      qsort(selected_labels, num_selected_dofs, sizeof(*selected_labels), AMGIntCompare);
   }
   for (HYPRE_Int i = 0; i < num_entries; i++)
   {
      if (num_selected_dofs > 0 &&
          bsearch(&dofmap->data[i], selected_labels, num_selected_dofs,
                  sizeof(*selected_labels), AMGIntCompare))
      {
         selected[i] = 1;
         projected_local++;
      }
   }

   MPI_Allreduce(&projected_local, &projected_global, 1, MPI_UINT64_T, MPI_SUM, comm);
   MPI_Scan(&projected_local, &projected_scan, 1, MPI_UINT64_T, MPI_SUM, comm);

   if (projected_global == 0)
   {
      hypredrv_ErrorCodeSet(ERROR_INVALID_VAL);
      hypredrv_ErrorMsgAdd("MGR label set selects no near-null-space entries");
   }
   else
   {
      int local_nonzero = 0;
      for (HYPRE_Int i = 0; i < num_entries && !local_nonzero; i++)
      {
         if (!selected[i])
         {
            continue;
         }
         for (HYPRE_Int mode = 0; mode < AMG_NUM_RBMS; mode++)
         {
            size_t source_index =
               ((size_t)(AMG_FIRST_RBM_COMPONENT + mode) * (size_t)vector_stride) +
               ((size_t)i * (size_t)index_stride);
            if (hypre_cabs(host_values[source_index]) > 0.0)
            {
               local_nonzero = 1;
               break;
            }
         }
      }
      int global_nonzero = 0;
      MPI_Allreduce(&local_nonzero, &global_nonzero, 1, MPI_INT, MPI_MAX, comm);
      if (global_nonzero)
      {
         AMGCreateRBMsFromHost(args, comm, (HYPRE_BigInt)projected_global,
                               (HYPRE_BigInt)(projected_scan - projected_local),
                               (HYPRE_Int)projected_local, num_entries, selected,
                               host_values, vector_stride, index_stride);
      }
      args->interp_vec_variant = 1;
      if (!hypredrv_ErrorCodeActive())
      {
         args->rbm_source_generation  = input_generation;
         args->rbm_source_labels_size = num_selected_dofs;
         args->rbm_source_labels_hash = labels_hash;
      }
   }

   free(selected);
   free(selected_labels);
   hypre_TFree(host_values, HYPRE_MEMORY_HOST);
#else
   (void)args;
   (void)vec_nn;
   (void)dofmap;
   (void)input_generation;
   (void)selected_dofs;
   (void)num_selected_dofs;
#endif
}

/*-----------------------------------------------------------------------------
 * hypredrv_AMGSetDofFunc
 *
 * Pass the user dofmap to BoomerAMG as the "dof_func" (unknown-to-function
 * map) when systems AMG is active (num_functions > 1). Without an explicit
 * dof_func, hypre assumes interleaved unknowns, i.e., dof_func[i] =
 * (i + offset) % num_functions, which is wrong for other unknown orderings
 * such as contiguous component blocks.
 *
 * hypre takes ownership of the array passed to HYPRE_BoomerAMGSetDofFunc.
 * From hypre 2.23.0, dof_func is wrapped in a hypre_IntArray whose memory
 * location is bound to the matrix's at setup time, so the copy of the dofmap
 * is allocated accordingly. Older versions store and free dof_func as a host
 * array, so the copy always lives on the host there.
 *
 * The dofmap is skipped (keeping hypre's default) when its labels do not form
 * a valid function map, i.e., when any global label falls outside
 * [0, num_functions). This can happen when a dofmap meant for MGR (more
 * labels than functions) is combined with a systems AMG preconditioner. The
 * check uses the global unique labels, which are identical on all ranks, so
 * the decision is rank-consistent without extra communication.
 *-----------------------------------------------------------------------------*/

void
hypredrv_AMGSetDofFunc(const AMG_args *args, const IntArray *dofmap, HYPRE_Solver precon,
                       HYPRE_IJMatrix ij_A)
{
   void     *vA            = NULL;
   HYPRE_Int num_functions = args->coarsening.num_functions;

   /* An empty local dofmap (a rank that owns no rows) is not a reason to skip this:
    * whether dof_func is attached changes the collective sequence of the AMG setup,
    * so the decision has to be the same on every rank. The remaining conditions below
    * are evaluated from globally reduced data and are rank-consistent by construction. */
   if (!precon || !ij_A || !dofmap || !dofmap->data || num_functions <= 1)
   {
      return;
   }

   if (dofmap->g_unique_size == 0 || !dofmap->g_unique_data ||
       dofmap->g_unique_data[0] < 0 ||
       dofmap->g_unique_data[dofmap->g_unique_size - 1] >= (int)num_functions)
   {
      HYPREDRV_LOG_COMMF(2, hypre_IJMatrixComm((hypre_IJMatrix *)ij_A), NULL, 0,
                         "dofmap labels do not fit in [0, num_functions=%d); "
                         "keeping hypre's default (interleaved) dof_func",
                         (int)num_functions);
      return;
   }

   HYPRE_IJMatrixGetObject(ij_A, &vA);
   hypre_ParCSRMatrix *par_A = (hypre_ParCSRMatrix *)vA;
   if (!par_A)
   {
      return;
   }

   HYPRE_Int num_local_rows = hypre_ParCSRMatrixNumRows(par_A);
   if ((size_t)num_local_rows != dofmap->size)
   {
      hypredrv_ErrorCodeSet(ERROR_INVALID_VAL);
      hypredrv_ErrorMsgAdd("Dofmap size (%d) does not match the number of local"
                           " matrix rows (%d)",
                           (int)dofmap->size, (int)num_local_rows);
      return;
   }

   /* Allocate at least one element so the attached dof_func pointer is non-NULL even
    * on a rank that owns no rows; hypre only tests the pointer, never reads past the
    * local row count. */
   HYPRE_Int *h_dof_func = hypre_TAlloc(
      HYPRE_Int, (num_local_rows > 0) ? num_local_rows : 1, HYPRE_MEMORY_HOST);
   for (HYPRE_Int i = 0; i < num_local_rows; i++)
   {
      h_dof_func[i] = (HYPRE_Int)dofmap->data[i];
   }

   HYPRE_Int *dof_func = h_dof_func;
#if defined(HYPRE_USING_GPU) && HYPRE_CHECK_MIN_VERSION(22300, 0)
   if (hypre_ParCSRMatrixMemoryLocation(par_A) == HYPRE_MEMORY_DEVICE)
   {
      dof_func = hypre_TAlloc(HYPRE_Int, num_local_rows, HYPRE_MEMORY_DEVICE);
      hypre_TMemcpy(dof_func, h_dof_func, HYPRE_Int, num_local_rows, HYPRE_MEMORY_DEVICE,
                    HYPRE_MEMORY_HOST);
      hypre_TFree(h_dof_func, HYPRE_MEMORY_HOST);
   }
#endif

   HYPRE_BoomerAMGSetDofFunc(precon, dof_func);
}

/*-----------------------------------------------------------------------------
 * AMGCreate
 *-----------------------------------------------------------------------------*/

void
hypredrv_AMGCreate(const AMG_args *args, HYPRE_Solver *precon_ptr)
{
   HYPRE_Solver precon = NULL;

   HYPRE_BoomerAMGCreate(&precon);
   HYPRE_BoomerAMGSetInterpType(precon, args->interpolation.prolongation_type);
   HYPRE_BoomerAMGSetRestriction(precon, args->interpolation.restriction_type);
#if HYPRE_CHECK_MIN_VERSION(21800, 0)
   /* Strength and filtering of the approximate ideal restriction operator */
   HYPRE_BoomerAMGSetStrongThresholdR(precon, args->interpolation.restrict_strong_th);
   HYPRE_BoomerAMGSetFilterThresholdR(precon, args->interpolation.restrict_filter_th);
#endif
   HYPRE_BoomerAMGSetCoarsenType(precon, args->coarsening.type);
   HYPRE_BoomerAMGSetSabs(precon, args->coarsening.sabs);
   HYPRE_BoomerAMGSetTol(precon, args->tolerance);
   HYPRE_BoomerAMGSetStrongThreshold(precon, args->coarsening.strong_th);
   HYPRE_BoomerAMGSetSeqThreshold(precon, args->coarsening.seq_amg_th);
   HYPRE_BoomerAMGSetMaxCoarseSize(precon, args->coarsening.max_coarse_size);
   HYPRE_BoomerAMGSetMinCoarseSize(precon, args->coarsening.min_coarse_size);
   HYPRE_BoomerAMGSetTruncFactor(precon, args->interpolation.trunc_factor);
   HYPRE_BoomerAMGSetPMaxElmts(precon, args->interpolation.max_nnz_row);
   HYPRE_BoomerAMGSetPrintLevel(precon, args->print_level);
   HYPRE_BoomerAMGSetChebyOrder(precon, args->relaxation.chebyshev.order);
   HYPRE_BoomerAMGSetChebyFraction(precon, args->relaxation.chebyshev.fraction);
   HYPRE_BoomerAMGSetChebyEigEst(precon, args->relaxation.chebyshev.eig_est);
   HYPRE_BoomerAMGSetChebyVariant(precon, args->relaxation.chebyshev.variant);
   HYPRE_BoomerAMGSetChebyScale(precon, args->relaxation.chebyshev.scale);
   if (args->relaxation.type >= 0)
   {
      HYPRE_BoomerAMGSetRelaxType(precon, args->relaxation.type);
   }
   HYPRE_BoomerAMGSetRelaxOrder(precon, args->relaxation.order);
   HYPRE_BoomerAMGSetRelaxWt(precon, args->relaxation.weight);
   HYPRE_BoomerAMGSetOuterWt(precon, args->relaxation.outer_weight);
   HYPRE_BoomerAMGSetMaxLevels(precon, args->coarsening.max_levels);
   HYPRE_BoomerAMGSetSmoothType(precon, args->smoother.type);
   HYPRE_BoomerAMGSetSmoothNumSweeps(precon, args->smoother.num_sweeps);
   HYPRE_BoomerAMGSetSmoothNumLevels(precon, args->smoother.num_levels);
   HYPRE_BoomerAMGSetMaxRowSum(precon, args->coarsening.max_row_sum);
#if HYPRE_CHECK_MIN_VERSION(22100, 0)
   HYPRE_BoomerAMGSetILUType(precon, args->smoother.ilu.type);
#endif
#if HYPRE_CHECK_MIN_VERSION(22600, 0)
#if HYPREDRV_HYPRE_RELEASE_NUMBER == 22900
   HYPRE_BoomerAMGSetILULocalReordering(
      precon, args->smoother.ilu.reordering ? args->smoother.ilu.reordering : 1);
#else
   HYPRE_BoomerAMGSetILULocalReordering(precon, args->smoother.ilu.reordering);
#endif
   HYPRE_BoomerAMGSetILUTriSolve(precon, args->smoother.ilu.tri_solve);
   HYPRE_BoomerAMGSetILULowerJacobiIters(precon, args->smoother.ilu.lower_jac_iters);
   HYPRE_BoomerAMGSetILUUpperJacobiIters(precon, args->smoother.ilu.upper_jac_iters);
#endif
#if HYPRE_CHECK_MIN_VERSION(22100, 0)
   HYPRE_BoomerAMGSetILULevel(precon, args->smoother.ilu.fill_level);
   HYPRE_BoomerAMGSetILUDroptol(precon, args->smoother.ilu.droptol);
   HYPRE_BoomerAMGSetILUMaxRowNnz(precon, args->smoother.ilu.max_row_nnz);
   HYPRE_BoomerAMGSetILUMaxIter(precon, args->smoother.num_sweeps);
#endif
#if HYPRE_CHECK_MIN_VERSION(23000, 0)
   HYPRE_BoomerAMGSetFSAIAlgoType(precon, args->smoother.fsai.algo_type);
   HYPRE_BoomerAMGSetFSAILocalSolveType(precon, args->smoother.fsai.ls_type);
   HYPRE_BoomerAMGSetFSAIMaxSteps(precon, args->smoother.fsai.max_steps);
   HYPRE_BoomerAMGSetFSAIMaxStepSize(precon, args->smoother.fsai.max_step_size);
   HYPRE_BoomerAMGSetFSAIMaxNnzRow(precon, args->smoother.fsai.max_nnz_row);
   HYPRE_BoomerAMGSetFSAINumLevels(precon, args->smoother.fsai.num_levels);
   HYPRE_BoomerAMGSetFSAIThreshold(precon, args->smoother.fsai.threshold);
   HYPRE_BoomerAMGSetFSAIEigMaxIters(precon, args->smoother.fsai.eig_max_iters);
   HYPRE_BoomerAMGSetFSAIKapTolerance(precon, args->smoother.fsai.kap_tolerance);
#endif
   HYPRE_BoomerAMGSetNumFunctions(precon, args->coarsening.num_functions);
#if HYPRE_CHECK_MIN_VERSION(23200, 0)
   HYPRE_BoomerAMGSetFilterFunctions(precon, args->coarsening.filter_functions);
#endif
   HYPRE_BoomerAMGSetAggNumLevels(precon, args->aggressive.num_levels);
   HYPRE_BoomerAMGSetAggInterpType(precon, args->aggressive.prolongation_type);
   HYPRE_BoomerAMGSetAggTruncFactor(precon, args->aggressive.trunc_factor);
   HYPRE_BoomerAMGSetAggP12TruncFactor(precon, args->aggressive.P12_trunc_factor);
   HYPRE_BoomerAMGSetAggPMaxElmts(precon, args->aggressive.max_nnz_row);
   HYPRE_BoomerAMGSetAggP12MaxElmts(precon, args->aggressive.P12_max_elements);
   HYPRE_BoomerAMGSetNumPaths(precon, args->aggressive.num_paths);
   HYPRE_BoomerAMGSetMaxIter(precon, args->max_iter);
   HYPRE_BoomerAMGSetRAP2(precon, args->coarsening.rap2);
#if HYPRE_CHECK_MIN_VERSION(21800, 0)
   HYPRE_BoomerAMGSetModuleRAP2(precon, args->coarsening.mod_rap2);
#endif
   HYPRE_BoomerAMGSetKeepTranspose(precon, args->coarsening.keep_transpose);

   /* Pre-smoothing */
   HYPRE_BoomerAMGSetCycleRelaxType(precon, args->relaxation.down_type, 1);
   if (args->relaxation.down_sweeps > -1)
   {
      HYPRE_BoomerAMGSetCycleNumSweeps(precon, args->relaxation.down_sweeps, 1);
   }
   else
   {
      HYPRE_BoomerAMGSetCycleNumSweeps(precon, args->relaxation.num_sweeps, 1);
   }

   /* Post-smoothing */
   HYPRE_BoomerAMGSetCycleRelaxType(precon, args->relaxation.up_type, 2);
   if (args->relaxation.up_sweeps > -1)
   {
      HYPRE_BoomerAMGSetCycleNumSweeps(precon, args->relaxation.up_sweeps, 2);
   }
   else
   {
      HYPRE_BoomerAMGSetCycleNumSweeps(precon, args->relaxation.num_sweeps, 2);
   }

   /* Coarsest level smoothing */
   HYPRE_BoomerAMGSetCycleRelaxType(precon, args->relaxation.coarse_type, 3);
   if (args->relaxation.coarse_sweeps > -1)
   {
      HYPRE_BoomerAMGSetCycleNumSweeps(precon, args->relaxation.coarse_sweeps, 3);
   }
   else
   {
      HYPRE_BoomerAMGSetCycleNumSweeps(precon, args->relaxation.num_sweeps, 3);
   }

   /* Per-sweep grid point schedule. The AIR algorithm skips the down cycle and
      relaxes F-points on the way up, finishing on the C-points when there are
      enough sweeps to do so. hypre takes ownership of the array. */
   if (args->relaxation.points == 1)
   {
      HYPRE_Int ns_down   = (args->relaxation.down_sweeps > -1)
                               ? args->relaxation.down_sweeps
                               : args->relaxation.num_sweeps;
      HYPRE_Int ns_up     = (args->relaxation.up_sweeps > -1) ? args->relaxation.up_sweeps
                                                              : args->relaxation.num_sweeps;
      HYPRE_Int ns_coarse = (args->relaxation.coarse_sweeps > -1)
                               ? args->relaxation.coarse_sweeps
                               : args->relaxation.num_sweeps;

      HYPRE_Int **grid_relax_points = hypre_CTAlloc(HYPRE_Int *, 4, HYPRE_MEMORY_HOST);

      grid_relax_points[0] = NULL;
      grid_relax_points[1] = hypre_CTAlloc(HYPRE_Int, ns_down, HYPRE_MEMORY_HOST);
      grid_relax_points[2] = hypre_CTAlloc(HYPRE_Int, ns_up, HYPRE_MEMORY_HOST);
      grid_relax_points[3] = hypre_CTAlloc(HYPRE_Int, ns_coarse, HYPRE_MEMORY_HOST);

      /* Down cycle and coarsest level relax every point */
      for (HYPRE_Int i = 0; i < ns_down; i++) grid_relax_points[1][i] = 0;
      for (HYPRE_Int i = 0; i < ns_coarse; i++) grid_relax_points[3][i] = 0;

      /* Up cycle: F-points, with a trailing C-point sweep when ns_up > 2 */
      for (HYPRE_Int i = 0; i < ns_up; i++) grid_relax_points[2][i] = -1;
      if (ns_up > 2) grid_relax_points[2][ns_up - 1] = 1;

      HYPRE_BoomerAMGSetGridRelaxPoints(precon, grid_relax_points);
   }

   if (args->coarsening.nodal)
   {
      /* Configurable GM/LN rigid-body-mode interpolation, used for both a
         top-level elasticity AMG and an MGR subblock AMG. The "auto" defaults
         reproduce the historical hard-coded setup; YAML can override them
         (coarsening.nodal_type / interp_vec_variant / smooth_interp_vecs /
         interp_vec_qmax / interp_vec_abs_q_trunc) to experiment toward an
         h-uniform elasticity AMG. */
      HYPRE_Int variant = (args->coarsening.interp_vec_variant > -1)
                             ? args->coarsening.interp_vec_variant
                             : args->interp_vec_variant;
      /* GM2 expands interpolation with the near-null-space vectors, so it wants a
         nonzero QMax and smoothed vectors; the other variants keep hypre's defaults. */
      HYPRE_Int qmax   = (variant == 2) ? 4 : 0;
      HYPRE_Int smooth = (variant == 2) ? 1 : 0;

      if (args->coarsening.interp_vec_qmax > -1)
      {
         qmax = args->coarsening.interp_vec_qmax;
      }
      if (args->coarsening.smooth_interp_vecs > -1)
      {
         smooth = args->coarsening.smooth_interp_vecs;
      }

      HYPRE_BoomerAMGSetNumFunctions(precon, 3);
      HYPRE_BoomerAMGSetNodal(precon, args->coarsening.nodal_type);
      if (args->coarsening.nodal_type == 4)
      {
         /* Nodal coarsening based on the row-sum norm */
         HYPRE_BoomerAMGSetNodalDiag(precon, 1);
      }
      HYPRE_BoomerAMGSetInterpVecVariant(precon, variant);
      HYPRE_BoomerAMGSetInterpVecQMax(precon, qmax);
      HYPRE_BoomerAMGSetInterpVecAbsQTrunc(precon,
                                           args->coarsening.interp_vec_abs_q_trunc);
#if HYPRE_CHECK_MIN_VERSION(30000, 0)
      HYPRE_BoomerAMGSetSmoothInterpVectors(precon, smooth);
#endif
      HYPRE_BoomerAMGSetInterpVectors(precon, args->num_rbms,
                                      (HYPRE_ParVector *)args->rbms);
   }

   *precon_ptr = precon;
}
