test_interp_operator_sum_mvp_at_tgt.c

Go to the documentation of this file.

// Copyright (c) 2025 The YAC Authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <math.h>
 
#include <mpi.h>
#include <yaxt.h>
 
#include "tests.h"
#include "interpolation/interpolation.h"
#include "yac_mpi.h"
#include "test_common.h"
#include "interpolation/interpolation_utils.h"
#include "interpolation/operators/interp_operator_sum_mvp_at_tgt.h"
#include "collection_selection_internal.h"
 
static void utest_init_data(
  struct yac_collection_selection * sel, double **** src_fields_collection,
  double **** src_fields_frac_mask_collection, double *** tgt_field_collection);
static void utest_free_data(
  struct yac_collection_selection * sel, double *** src_fields_collection,
  double *** src_fields_frac_mask_collection, double ** tgt_field_collection);
static void utest_check_sum_mvp_at_tgt(struct yac_collection_selection * sel,
                                       Xt_redist * src_redists,
                                       size_t * num_src_points_per_field,
                                       size_t * tgt_pos, size_t tgt_count,
                                       size_t * num_src_per_tgt, double * weights,
                                       size_t * src_field_idx, size_t * src_idx,
                                       size_t num_src_fields,
                                       double frac_mask_fallback_value,
                                       double *** src_fields_collection,
                                       double *** src_fields_frac_mask_collection,
                                       double ** tgt_field_collection);
 
enum {
  MAX_COLLECTION_SIZE = 4,
  MAX_NUM_SRC_FIELDS = 3,
  NUM_SRC_POINTS = 4,
  NUM_TGT_POINTS = 4,
  NUM_PROCS = 2,
};
#define TGT_UNSET_VALUE (-1.0)
 
int main(void) {
 
  MPI_Init(NULL, NULL);
  xt_initialize(MPI_COMM_WORLD);
 
  int comm_rank, comm_size;
  MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
 
  if (comm_size != NUM_PROCS) {
    PUT_ERR("ERROR: test requires 2 processes");
    xt_finalize();
    MPI_Finalize();
    return TEST_EXIT_CODE;
  }
 
  // collection selection configs (only applied to source field data,
  // target field data is expected to have a contiguous collection of size
  // collection_size)
  struct {
    size_t N;
    size_t * indices;
  } collection_selection_configs[] =
    {// selects three entries from the collection in random order
     {.N = 3, .indices = (size_t[]){0, 3, 2}},
     // selects all entries from the collection in ascending order
     {.N = 4, .indices = (size_t[]){0, 1, 2, 3}},
     // selects all entries from the collection in default (ascending) order
     {.N = MAX_COLLECTION_SIZE, .indices = NULL}};
  enum { NUM_COLL_SEL_CONFIGS =
           sizeof(collection_selection_configs) /
           sizeof(collection_selection_configs[0]) };
 
  // loop over collection selections
  for (size_t sel_idx = 0; sel_idx < NUM_COLL_SEL_CONFIGS; ++sel_idx) {
 
    struct yac_collection_selection * sel =
      yac_collection_selection_new(
        collection_selection_configs[sel_idx].N,
        collection_selection_configs[sel_idx].indices);
 
    // interpolation configurations
    // (basically contains entries of the sparsely populated interpolation
    //  weight matrix)
    // We assume that the global ids of the source source points of each field
    // are:
    // * for rank 0: [0, 1, 2, 3]
    // * for rank 1: [4, 5, 6, 7]
    struct {
      struct {
 
        //----------------------------------------------------------------------
        // Information about the target points that are to be interpolated
        //----------------------------------------------------------------------
 
        size_t * tgt_pos;               // indices of all interpolated target
                                        // points
        size_t tgt_count;               // number of interpolated target points
 
        //----------------------------------------------------------------------
        // Information about the source points from remote processes required
        // the interpolation of the local target points
        //----------------------------------------------------------------------
 
        Xt_int * src_global_ids_per_field[MAX_NUM_SRC_FIELDS];
                                        // global ids of all source points of
                                        // each source field from remote
                                        // processes required for the
                                        // interpolation of the local target
                                        // points
        size_t num_src_points_per_field[MAX_NUM_SRC_FIELDS];
                                        // number of source points per source
                                        // field required from remote processes
 
        //----------------------------------------------------------------------
        // Information on how the target points are to be interpolated
        //
        // Information on how to interpolate the i'th target
        // (with i being in the range [0;tgt_count-1]) is stored at:
        //   * num_src_per_tgt[i]
        //   * weights[SUM(num_src_per_tgt[0..i])..SUM(num_src_per_tgt[0..i+1]))
        //   * src_field_idx[SUM(num_src_per_tgt[0..i])..SUM(num_src_per_tgt[0..i+1]))
        //   * src_idx[SUM(num_src_per_tgt[0..i])..SUM(num_src_per_tgt[0..i+1]))
        //
        // A value of SIZE_MAX in the src_field_idx array indicates a source
        // field point from a remote process. In that case, the respective
        // src_idx in the contiguous receive buffer containing received
        // source points for all fields of a collection entry.
        //----------------------------------------------------------------------
 
        size_t * num_src_per_tgt;       // number of source points per target
        double * weights;               // weights for the interpolation of
                                        // the target points
        size_t * src_field_idx;         // source field indices of source points
                                        // required for the interpolation of the
                                        // target points.
        size_t * src_idx;               // source indices / receive buffer index
 
        //----------------------------------------------------------------------
        // reference data
        //----------------------------------------------------------------------
        int is_source;
        int is_target;
 
      } interp_data[NUM_PROCS];
      size_t num_src_fields;            // number of source fields
    } interp_configs[] =
      // no interpolation
      // (both processes are neither source nor target)
      {{.interp_data =
          {{.src_global_ids_per_field = {NULL},
            .num_src_points_per_field = {0},
            .tgt_pos = NULL,
            .tgt_count = 0,
            .num_src_per_tgt = NULL,
            .weights = NULL,
            .src_field_idx = NULL,
            .src_idx = NULL,
            .is_source = 0,
            .is_target = 0},
           {.src_global_ids_per_field = {NULL},
            .num_src_points_per_field = {0},
            .tgt_pos = NULL,
            .tgt_count = 0,
            .num_src_per_tgt = NULL,
            .weights = NULL,
            .src_field_idx = NULL,
            .src_idx = NULL,
            .is_source = 0,
            .is_target = 0}},
        .num_src_fields = 1},
       // target data contains the source points of the respective other rank
       // and from the local one
       // (rank 0 and rank 1 are source and rank 1 is only target)
       // dimensions (for documentation only):
       //   tgt[rank][local_idx]
       //   src[rank][src_field_idx][local_idx]
       //
       // tgt[1][1] = 0.5 * src[0][1][0] + 0.5 * src[1][1][0]
       // tgt[1][3] = 0.5 * src[0][1][1] + 0.5 * src[1][1][1]
       {.interp_data =
          {{.src_global_ids_per_field = {NULL, NULL, NULL},
            .num_src_points_per_field = {0,0,0},
            .tgt_pos = NULL,
            .tgt_count = 0,
            .num_src_per_tgt = NULL,
            .weights = NULL,
            .src_field_idx = NULL,
            .src_idx = NULL,
            .is_source = 1,
            .is_target = 0},
           {.src_global_ids_per_field = {NULL, (Xt_int[]){0,1}, NULL},
            .num_src_points_per_field = {0,2,0},
            .tgt_pos = (size_t[]){1,3},
            .tgt_count = 2,
            .num_src_per_tgt = (size_t[]){2,2},
            .weights = (double[]){0.5,0.5, 0.5,0.5},
            .src_field_idx = (size_t[]){SIZE_MAX, 1, SIZE_MAX, 1},
            .src_idx = (size_t[]){0,0,1,1},
            .is_source = 1,
            .is_target = 1}},
        .num_src_fields = 3},
       // target data contains the source points of the respective other rank
       // (both processes are source and target;
       //  source points are only from remote process)
       {.interp_data =
          {{.src_global_ids_per_field = {(Xt_int[]){4,5,6,7}},
            .num_src_points_per_field = {4},
            .tgt_pos = (size_t[]){0,1,2,3},
            .tgt_count = 4,
            .num_src_per_tgt = (size_t[]){1,1,1,1},
            .weights = (double[]){1.0, 1.0, 1.0, 1.0},
            .src_field_idx = (size_t[]){SIZE_MAX, SIZE_MAX, SIZE_MAX, SIZE_MAX},
            .src_idx = (size_t[]){0,1,2,3},
            .is_source = 1,
            .is_target = 1},
           {.src_global_ids_per_field = {(Xt_int[]){0,1,2,3}},
            .num_src_points_per_field = {4},
            .tgt_pos = (size_t[]){0,1,2,3},
            .tgt_count = 4,
            .num_src_per_tgt = (size_t[]){1,1,1,1},
            .weights = (double[]){1.0, 1.0, 1.0, 1.0},
            .src_field_idx = (size_t[]){SIZE_MAX, SIZE_MAX, SIZE_MAX, SIZE_MAX},
            .src_idx = (size_t[]){0,1,2,3},
            .is_source = 1,
            .is_target = 1}},
        .num_src_fields = 1},
       // more complex combination of the source data from multiple fields
       // dimensions (for documentation only):
       //   tgt[rank][local_idx]
       //   src[rank][src_field_idx][local_idx]
       //
       // tgt[0][0] = 0.2 * src[0][0][1] + 0.4 src[1][1][0] + 0.4 * src[1][2][3]
       // tgt[0][2] = sum(src[0][1][:]) + sum(src[1][2][:])
       // tgt[1][1] = 0.5 * src[0][2][1] + 0.5 * src[0][2][3]
       // tgt[1][3] = 0.1 * src[0][0][0] +
       //             0.1 * src[0][0][2] +
       //             0.1 * src[0][1][2] +
       //             0.1 * src[0][2][2] +
       //             0.1 * src[1][0][2] +
       //             0.1 * src[1][1][2] +
       //             0.1 * src[1][2][2]
       // (both processes are source and target;
       //  source points are both from local and remote process)
       {.interp_data =
          {{.src_global_ids_per_field = {NULL,
                                         (Xt_int[]){4},
                                         (Xt_int[]){4,5,6,7}},
            .num_src_points_per_field = {0,1,4},
            .tgt_pos = (size_t[]){0,2},
            .tgt_count = 2,
            .num_src_per_tgt = (size_t[]){3,8},
            .weights = (double[]){0.2,0.4,0.4,
                                  1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0},
            .src_field_idx = (size_t[]){0,SIZE_MAX,SIZE_MAX,
                                        1,1,1,1,
                                          SIZE_MAX,SIZE_MAX,SIZE_MAX,SIZE_MAX},
            .src_idx = (size_t[]){1,0,4,
                                  0,1,2,3,1,2,3,4},
            .is_source = 1,
            .is_target = 1},
           {.src_global_ids_per_field = {(Xt_int[]){0,2},
                                         (Xt_int[]){2},
                                         (Xt_int[]){1,2,3}},
            .num_src_points_per_field = {2,1,3},
            .tgt_pos = (size_t[]){1,3},
            .tgt_count = 2,
            .num_src_per_tgt = (size_t[]){2,7},
            .weights = (double[]){0.5,0.5,
                                  0.1,0.1,0.1,0.1,0.1,0.1,0.1},
            .src_field_idx = (size_t[]){SIZE_MAX,SIZE_MAX,
                                        SIZE_MAX,SIZE_MAX,SIZE_MAX,SIZE_MAX,
                                          0,1,2},
            .src_idx = (size_t[]){3,5,
                                  0,1,2,4,2,2,2},
            .is_source = 1,
            .is_target = 1}},
        .num_src_fields = 3},
       // more complex combination of the source data from multiple fields
       // dimensions (for documentation only):
       //   tgt[rank][local_idx]
       //   src[rank][src_field_idx][local_idx]
       //
       // tgt[0][0] = 0.5 * src[0][0][0] + 0.5 * src[0][0][0]
       // tgt[0][1] = 0.5 * src[0][0][1] + 0.5 * src[0][1][1]
       // tgt[0][2] = 0.5 * src[0][1][2] + 0.5 * src[0][0][2]
       // tgt[0][3] = 0.5 * src[0][1][3] + 0.5 * src[0][1][3]
       // tgt[1][0] = 0.5 * src[1][0][0] + 0.5 * src[1][1][0]
       // tgt[1][2] = 0.5 * src[1][0][2] + 0.5 * src[1][1][2]
       // (both processes are source and target;
       //  source points are only from local)
       {.interp_data =
          {{.src_global_ids_per_field = {NULL, NULL},
            .num_src_points_per_field = {0,0},
            .tgt_pos = (size_t[]){0,1,2,3},
            .tgt_count = 4,
            .num_src_per_tgt = (size_t[]){2,2,2,2},
            .weights = (double[]){0.5,0.5, 0.5,0.5, 0.5,0.5, 0.5,0.5},
            .src_field_idx = (size_t[]){0,0, 0,1, 1,0, 1,1},
            .src_idx = (size_t[]){0,0, 1,1, 2,2, 3,3},
            .is_source = 1,
            .is_target = 1},
           {.src_global_ids_per_field = {NULL, NULL},
            .num_src_points_per_field = {0,0},
            .tgt_pos = (size_t[]){0,2},
            .tgt_count = 2,
            .num_src_per_tgt = (size_t[]){2,2},
            .weights = (double[]){0.5,0.5, 0.5,0.5},
            .src_field_idx = (size_t[]){0,1, 0,1},
            .src_idx = (size_t[]){0,0, 2,2},
            .is_source = 1,
            .is_target = 1}},
        .num_src_fields = 2}};
    enum {
      NUM_INTERP_CONFIGS = sizeof(interp_configs) / sizeof(interp_configs[0])};
 
    // generate yaxt index list containing global ids of all locally available
    // source points (each source field contains the same global ids)
    // each process holds the following global ids:
    //  i + comm_rank * NUM_SRC_POINTS
    //  with i in [0..NUM_SRC_POINTS]
    Xt_int src_global_ids[NUM_SRC_POINTS];
    for (size_t src_idx = 0; src_idx < NUM_SRC_POINTS; ++src_idx)
      src_global_ids[src_idx] =
        (Xt_int)(src_idx + comm_rank * NUM_SRC_POINTS);
    Xt_idxlist src_idxlist =
      xt_idxvec_new(src_global_ids, NUM_SRC_POINTS);
 
    // loop over all interpolation configurations
    for (size_t interp_config_idx = 0; interp_config_idx < NUM_INTERP_CONFIGS;
         ++interp_config_idx) {
 
      // generate input data for yac_interpolation_sum_mvp_at_tgt_new from list
      // of links
      size_t num_src_fields =
        interp_configs[interp_config_idx].num_src_fields;
      Xt_redist src_redists[num_src_fields];
      for (size_t src_field_idx = 0; src_field_idx < num_src_fields;
           ++src_field_idx) {
 
        // create redist for receiving source data of the current source field
        // required for the interpolation of the local target points
        Xt_idxlist required_source_idxlist =
          (interp_configs[interp_config_idx].
            interp_data[comm_rank].
              num_src_points_per_field[src_field_idx] > 0)?
            xt_idxvec_new(
              interp_configs[interp_config_idx].
                interp_data[comm_rank].
                  src_global_ids_per_field[src_field_idx],
              interp_configs[interp_config_idx].
                interp_data[comm_rank].
                  num_src_points_per_field[src_field_idx]):
            xt_idxempty_new();
        Xt_xmap xmap =
          xt_xmap_all2all_new(
            src_idxlist, required_source_idxlist, MPI_COMM_WORLD);
        src_redists[src_field_idx] = xt_redist_p2p_new(xmap, MPI_DOUBLE);
        xt_xmap_delete(xmap);
        xt_idxlist_delete(required_source_idxlist);
 
      } // src_field_idx
 
      for (int with_frac_mask = 0; with_frac_mask <= 1; ++with_frac_mask) {
 
        // if without fractional masking, set fallback value to
        // YAC_FRAC_MASK_NO_VALUE, which deactivates it
        double frac_mask_fallback_value =
          with_frac_mask?-10.0:YAC_FRAC_MASK_NO_VALUE;
 
        for (int use_weights = 0; use_weights <= 1; ++use_weights) {
 
          // if without weights, set it to NULL, which deactivates it
          double * weights =
            use_weights?
              interp_configs[interp_config_idx].interp_data[comm_rank].weights:
              NULL;
 
          // some consistency checking of the input data
          {
            size_t total_num_receive_src_points = 0;
            for (size_t i = 0;
                 i < interp_configs[interp_config_idx].
                       interp_data[comm_rank].tgt_count; ++i) {
              total_num_receive_src_points +=
                interp_configs[interp_config_idx].
                  interp_data[comm_rank].num_src_per_tgt[i];
            }
            for (size_t i = 0; i < total_num_receive_src_points; ++i) {
              YAC_ASSERT(
                (interp_configs[interp_config_idx].
                   interp_data[comm_rank].src_field_idx[i] == SIZE_MAX) ||
                (interp_configs[interp_config_idx].
                   interp_data[comm_rank].src_field_idx[i] < num_src_fields),
                "ERROR in src_field_idx");
              YAC_ASSERT(
                interp_configs[interp_config_idx].
                  interp_data[comm_rank].src_idx[i] <
                total_num_receive_src_points,
                "ERROR in src_idx");
            }
          }
 
          // generate interpolation for the current configuration
          struct yac_interp_operator * interp =
            yac_interp_operator_sum_mvp_at_tgt_new(
              sel, src_redists,
              interp_configs[interp_config_idx].interp_data[comm_rank].tgt_pos,
              interp_configs[interp_config_idx].interp_data[comm_rank].tgt_count,
              interp_configs[interp_config_idx].interp_data[comm_rank].num_src_per_tgt,
              weights,
              interp_configs[interp_config_idx].interp_data[comm_rank].src_field_idx,
              interp_configs[interp_config_idx].interp_data[comm_rank].src_idx,
              num_src_fields, with_frac_mask);
 
          // test interpolation built above and a copy of it, to make sure that
          // the copy-function works properly
          for (int use_copy = 0; use_copy <= 1; ++use_copy) {
 
            if (use_copy) {
              struct yac_interp_operator * interp_copy =
                yac_interp_operator_copy(interp);
              yac_interp_operator_delete(interp);
              interp = interp_copy;
            }
 
            { // check is_source and is_target
 
              if (yac_interp_operator_is_source(interp) !=
                  interp_configs[interp_config_idx].
                    interp_data[comm_rank].is_source)
                PUT_ERR("ERROR in is_source");
 
              if (yac_interp_operator_is_target(interp) !=
                  interp_configs[interp_config_idx].
                    interp_data[comm_rank].is_target)
                PUT_ERR("ERROR in is_target");
            }
 
            { // synchronous test (put and get in a single exchange call)
 
              // initialise source and target fields
              double *** src_fields_collection;
              double *** src_fields_frac_mask_collection;
              double ** tgt_field_collection;
              utest_init_data(
                sel, &src_fields_collection, &src_fields_frac_mask_collection,
                &tgt_field_collection);
 
              yac_interp_operator_execute(
                interp, src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection, frac_mask_fallback_value, 1.0, 0.0);
 
              // there should be no open put or get operation
              if (yac_interp_operator_execute_put_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_put_test");
              if (yac_interp_operator_execute_get_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_get_test");
 
              // check results
              utest_check_sum_mvp_at_tgt(
                sel, src_redists,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_pos,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_count,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_idx,
                num_src_fields, frac_mask_fallback_value,
                src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
 
              // free source and target fields
              utest_free_data(
                sel, src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
            }
 
            { // independent put + get
 
              // initialise source and target fields
              double *** src_fields_collection;
              double *** src_fields_frac_mask_collection;
              double ** tgt_field_collection;
              utest_init_data(
                sel, &src_fields_collection, &src_fields_frac_mask_collection,
                &tgt_field_collection);
 
              yac_interp_operator_execute_put(
                interp, src_fields_collection,
                src_fields_frac_mask_collection,
                1, frac_mask_fallback_value, 1.0, 0.0);
              yac_interp_operator_execute_get(
                interp, tgt_field_collection,
                frac_mask_fallback_value, 1.0, 0.0);
 
              // there should be no open put or get operation
              if (yac_interp_operator_execute_put_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_put_test");
              if (yac_interp_operator_execute_get_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_get_test");
 
              // check results
              utest_check_sum_mvp_at_tgt(
                sel, src_redists,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_pos,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_count,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_idx,
                num_src_fields, frac_mask_fallback_value,
                src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
 
              // free source and target fields
              utest_free_data(
                sel, src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
            }
 
            { // independent put + async get + wait
 
              // initialise source and target fields
              double *** src_fields_collection;
              double *** src_fields_frac_mask_collection;
              double ** tgt_field_collection;
              utest_init_data(
                sel, &src_fields_collection, &src_fields_frac_mask_collection,
                &tgt_field_collection);
 
              yac_interp_operator_execute_put(
                interp, src_fields_collection,
                src_fields_frac_mask_collection,
                1, frac_mask_fallback_value, 1.0, 0.0);
              yac_interp_operator_execute_get_async(
                interp, tgt_field_collection,
                frac_mask_fallback_value, 1.0, 0.0);
              yac_interp_operator_execute_wait(interp);
 
              // there should be no open put or get operation
              if (yac_interp_operator_execute_put_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_put_test");
              if (yac_interp_operator_execute_get_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_get_test");
 
              // check results
              utest_check_sum_mvp_at_tgt(
                sel, src_redists,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_pos,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_count,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_idx,
                num_src_fields, frac_mask_fallback_value,
                src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
 
              // free source and target fields
              utest_free_data(
                sel, src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
            }
 
            { // independent put + async get + test_get-loop
 
              // initialise source and target fields
              double *** src_fields_collection;
              double *** src_fields_frac_mask_collection;
              double ** tgt_field_collection;
              utest_init_data(
                sel, &src_fields_collection, &src_fields_frac_mask_collection,
                &tgt_field_collection);
 
              yac_interp_operator_execute_put(
                interp, src_fields_collection,
                src_fields_frac_mask_collection,
                1, frac_mask_fallback_value, 1.0, 0.0);
              yac_interp_operator_execute_get_async(
                interp, tgt_field_collection,
                frac_mask_fallback_value, 1.0, 0.0);
              while(
                yac_interp_operator_execute_get_test(interp) ==
                YAC_INTERP_INCOMPLETE);
 
              // there should be no open put or get operation
              if (yac_interp_operator_execute_put_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_put_test");
              if (yac_interp_operator_execute_get_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_get_test");
 
              // check results
              utest_check_sum_mvp_at_tgt(
                sel, src_redists,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_pos,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_count,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_idx,
                num_src_fields, frac_mask_fallback_value,
                src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
 
              // free source and target fields
              utest_free_data(
                sel, src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
            }
 
            { // independent async get + put + wait
 
              // initialise source and target fields
              double *** src_fields_collection;
              double *** src_fields_frac_mask_collection;
              double ** tgt_field_collection;
              utest_init_data(
                sel, &src_fields_collection, &src_fields_frac_mask_collection,
                &tgt_field_collection);
 
              yac_interp_operator_execute_get_async(
                interp, tgt_field_collection,
                frac_mask_fallback_value, 1.0, 0.0);
              yac_interp_operator_execute_put(
                interp, src_fields_collection,
                src_fields_frac_mask_collection,
                1, frac_mask_fallback_value, 1.0, 0.0);
              yac_interp_operator_execute_wait(interp);
 
              // there should be no open put or get operation
              if (yac_interp_operator_execute_put_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_put_test");
              if (yac_interp_operator_execute_get_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_get_test");
 
              // check results
              utest_check_sum_mvp_at_tgt(
                sel, src_redists,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_pos,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_count,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_idx,
                num_src_fields, frac_mask_fallback_value,
                src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
 
              // free source and target fields
              utest_free_data(
                sel, src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
            }
 
            { // regression test: async get with freed pointer array
              // (mimics yac_cget_async_ which frees the pointer array
              // immediately after execute_get_async returns)
              //
              // get_async is called BEFORE put so that the exchange has not
              // yet started when the pointer array is freed; this guarantees
              // that compute_tgt runs after the free and reliably exposes the
              // use-after-free bug that was present before the fix.
 
              // initialise source and target fields
              double *** src_fields_collection;
              double *** src_fields_frac_mask_collection;
              double ** tgt_field_collection;
              utest_init_data(
                sel, &src_fields_collection, &src_fields_frac_mask_collection,
                &tgt_field_collection);
 
              // build a transient pointer array over tgt_field_collection,
              // as yac_get_ does in yac.c
              size_t collection_size =
                yac_collection_selection_get_collection_size(sel);
              double ** tgt_field_ptr =
                xmalloc(collection_size * sizeof(*tgt_field_ptr));
              for (size_t c = 0; c < collection_size; ++c)
                tgt_field_ptr[c] = tgt_field_collection[c];
 
              yac_interp_operator_execute_get_async(
                interp, tgt_field_ptr,
                frac_mask_fallback_value, 1.0, 0.0);
 
              // free the transient pointer array immediately, reproducing the
              // use-after-free that occurred before the fix
              free(tgt_field_ptr);
 
              yac_interp_operator_execute_put(
                interp, src_fields_collection,
                src_fields_frac_mask_collection,
                1, frac_mask_fallback_value, 1.0, 0.0);
              yac_interp_operator_execute_wait(interp);
 
              // there should be no open put or get operation
              if (yac_interp_operator_execute_put_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_put_test");
              if (yac_interp_operator_execute_get_test(interp) !=
                  YAC_INTERP_COMPLETE)
                PUT_ERR("ERROR in execute_get_test");
 
              // check results
              utest_check_sum_mvp_at_tgt(
                sel, src_redists,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_pos,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_count,
                interp_configs[interp_config_idx].interp_data[comm_rank].num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].src_idx,
                num_src_fields, frac_mask_fallback_value,
                src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
 
              // free source and target fields
              utest_free_data(
                sel, src_fields_collection, src_fields_frac_mask_collection,
                tgt_field_collection);
            }
 
          } // use_copy
 
          yac_interp_operator_delete(interp);
        } // use_weights
      } // with_frac_mask
 
      for (size_t f = 0; f < num_src_fields; ++f)
        xt_redist_delete(src_redists[f]);
 
    } // interp_config_idx
 
    xt_idxlist_delete(src_idxlist);
    yac_collection_selection_delete(sel);
 
  } // sel_idx
 
  xt_finalize();
  MPI_Finalize();
  return TEST_EXIT_CODE;
}
 
static void utest_init_data(
  struct yac_collection_selection * sel, double **** src_fields_collection,
  double **** src_fields_frac_mask_collection,
  double *** tgt_field_collection) {
 
  size_t collection_size = yac_collection_selection_get_collection_size(sel);
 
  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
  *src_fields_collection =
    xmalloc(MAX_COLLECTION_SIZE * sizeof(**src_fields_collection));
  *src_fields_frac_mask_collection =
    xmalloc(MAX_COLLECTION_SIZE * sizeof(**src_fields_frac_mask_collection));
 
  for (size_t c = 0; c < MAX_COLLECTION_SIZE; ++c) {
 
    (*src_fields_collection)[c] =
      xmalloc(MAX_NUM_SRC_FIELDS * sizeof(***src_fields_collection));
    (*src_fields_frac_mask_collection)[c] =
      xmalloc(MAX_NUM_SRC_FIELDS * sizeof(***src_fields_frac_mask_collection));
 
    for (size_t f = 0; f < MAX_NUM_SRC_FIELDS; ++f) {
 
      (*src_fields_collection)[c][f] =
        xmalloc(NUM_SRC_POINTS * sizeof(***src_fields_collection));
      (*src_fields_frac_mask_collection)[c][f] =
        xmalloc(NUM_SRC_POINTS * sizeof(***src_fields_frac_mask_collection));
 
      for (size_t i = 0; i < NUM_SRC_POINTS; ++i) {
        (*src_fields_collection)[c][f][i] =
          (double)(c*1000 + f*100 + rank*10 + i);
        (*src_fields_frac_mask_collection)[c][f][i] = 1.0;
      }
    }
  }
 
  *tgt_field_collection =
    xmalloc(collection_size * sizeof(**tgt_field_collection));
 
  for (size_t c = 0; c < collection_size; ++c) {
 
    (*tgt_field_collection)[c] =
      xmalloc(NUM_TGT_POINTS * sizeof(**tgt_field_collection));
 
    for (size_t i = 0; i < NUM_TGT_POINTS; ++i) {
      (*tgt_field_collection)[c][i] = TGT_UNSET_VALUE;
    }
  }
}
 
static void utest_free_data(
  struct yac_collection_selection * sel, double *** src_fields_collection,
  double *** src_fields_frac_mask_collection, double ** tgt_field_collection) {
 
  size_t collection_size = yac_collection_selection_get_collection_size(sel);
 
  for (size_t c = 0; c < MAX_COLLECTION_SIZE; ++c) {
 
    for (size_t f = 0; f < MAX_NUM_SRC_FIELDS; ++f) {
 
      free(src_fields_collection[c][f]);
      free(src_fields_frac_mask_collection[c][f]);
    }
 
    free(src_fields_collection[c]);
    free(src_fields_frac_mask_collection[c]);
  }
 
  free(src_fields_collection);
  free(src_fields_frac_mask_collection);
 
  for (size_t c = 0; c < collection_size; ++c) {
 
    free(tgt_field_collection[c]);
  }
 
  free(tgt_field_collection);
}
 
static void utest_check_sum_mvp_at_tgt(struct yac_collection_selection * sel,
                                       Xt_redist * src_redists,
                                       size_t * num_src_points_per_field,
                                       size_t * tgt_pos, size_t tgt_count,
                                       size_t * num_src_per_tgt, double * weights,
                                       size_t * src_field_idx, size_t * src_idx,
                                       size_t num_src_fields,
                                       double frac_mask_fallback_value,
                                       double *** src_fields_collection,
                                       double *** src_fields_frac_mask_collection,
                                       double ** tgt_field_collection) {
 
  size_t collection_size = yac_collection_selection_get_collection_size(sel);
  size_t const * collection_indices = yac_collection_selection_get_indices(sel);
 
  // initialise reference target field with dummy values
  // and generate array containg pointers to the field associated to
  // each collection
  double ** ref_tgt_field_collection =
    xmalloc(collection_size * sizeof(*ref_tgt_field_collection));
  for (size_t c = 0; c < collection_size; ++c) {
    ref_tgt_field_collection[c] =
      xmalloc(NUM_TGT_POINTS * sizeof(**ref_tgt_field_collection));
    for (size_t i = 0; i < NUM_TGT_POINTS; ++i) {
      ref_tgt_field_collection[c][i] = TGT_UNSET_VALUE;
    }
  }
 
  // compute prefix sum of num_src_per_tgt, which is required
  // by compute_tgt_field_stencil
  size_t * prefix_num_src_per_tgt =
    xmalloc((tgt_count+1) * sizeof(*prefix_num_src_per_tgt));
  prefix_num_src_per_tgt[0] = 0;
  for (size_t i = 0, accu = 0; i < tgt_count; ++i) {
    accu += num_src_per_tgt[i];
    prefix_num_src_per_tgt[i+1] = accu;
  }
 
  // Generate contiguous view of selected source collection selection
  // (This simplifies further steps.)
  double *** cont_src_fields_collection =
    xmalloc(collection_size * sizeof(*cont_src_fields_collection));
  double *** cont_src_fields_frac_mask_collection =
    xmalloc(collection_size * sizeof(*src_fields_frac_mask_collection));
  for (size_t c = 0; c < collection_size; ++c) {
    // get collection index (if non-contiguous) otherwise use c
    size_t c_idx = (collection_indices == NULL) ? c : collection_indices[c];
    cont_src_fields_collection[c] =
      src_fields_collection[c_idx];
    cont_src_fields_frac_mask_collection[c] =
      src_fields_frac_mask_collection[c_idx];
  }
 
  // setup pointers for receiving source data from other processes
  // (since these arrays are used during an intermediate step, contiguous
  //  collection indices can be assumed)
  // (compute_tgt_field_stencil assumes that received source points from all fields
  //  is in one contiguous buffer, which is why pointers stored in
  //  remote_src_fields and remote_src_frac_masks are computed based on
  //  num_src_points_per_field)
  enum {SRC_FIELD = 0, SRC_FRAC_MASK = 1};
  double * remote_src_buffer[2][collection_size];
  double * remote_src_fields[collection_size][num_src_fields];
  double * remote_src_frac_masks[collection_size][num_src_fields];
  for (size_t i = 0; i < 2; ++i) {
    for (size_t c = 0; c < collection_size; ++c) {
      // allocate buffer to the maximum possible buffer size
      remote_src_buffer[i][c] =
        xmalloc(
          num_src_fields * (NUM_PROCS - 1) * NUM_SRC_POINTS *
          sizeof(remote_src_buffer[i][c][0]));
      for (size_t k = 0;
           k < num_src_fields * (NUM_PROCS - 1) * NUM_SRC_POINTS; ++k) {
        // initialise buffer with NAN to make error detection easier
        remote_src_buffer[i][c][k] = NAN;
      }
    }
  }
  for (size_t c = 0; c < collection_size; ++c) {
    for (size_t f = 0, offset = 0; f < num_src_fields; ++f) {
      remote_src_fields[c][f] =
        &remote_src_buffer[SRC_FIELD][c][offset];
      remote_src_frac_masks[c][f] =
        &remote_src_buffer[SRC_FRAC_MASK][c][offset];
      offset += num_src_points_per_field[f];
    }
  }
 
  // for each collection entry
  for (size_t c = 0; c < collection_size; ++c) {
 
    // get collection index (if non-contiguous) otherwise use c
    size_t c_idx = (collection_indices == NULL) ? c : collection_indices[c];
 
    // for each source field
    for (size_t f = 0; f < num_src_fields; ++f) {
 
      // receive source point and fractional mask values required for the
      // interpolation of the local target points
      xt_redist_s_exchange1(
        src_redists[f],
        (const void *)(src_fields_collection[c_idx][f]),
        (void *)(remote_src_fields[c][f]));
      xt_redist_s_exchange1(
        src_redists[f],
        (const void *)(src_fields_frac_mask_collection[c_idx][f]),
        (void *)(remote_src_frac_masks[c][f]));
    }
  }
 
  // compute target field reference
  compute_tgt_field_stencil(
    (double const * restrict **)cont_src_fields_collection,
    (double const * restrict **)cont_src_fields_frac_mask_collection,
    (double const **)remote_src_fields,
    (double const **)remote_src_frac_masks,
    ref_tgt_field_collection,
    tgt_pos, tgt_count, prefix_num_src_per_tgt, weights,
    src_field_idx, src_idx, num_src_fields, collection_size,
    frac_mask_fallback_value, 1.0, 0.0);
 
  // check target field results generate by interpolation operation against
  // refernce (if a target point is not being interpolated by the
  // interpolation operation, the target field and the reference should both
  // contain the value TGT_UNSET_VALUE)
  for (size_t c = 0; c < collection_size; ++c) {
    for (size_t i = 0; i < NUM_TGT_POINTS; ++i) {
      if (fabs(tgt_field_collection[c][i] -
               ref_tgt_field_collection[c][i]) > 1e-9) {
        PUT_ERR("wrong data in sum_mvp_at_tgt interpolation");
      }
    }
  }
 
  // clean up
  free(prefix_num_src_per_tgt);
  free(cont_src_fields_collection);
  free(cont_src_fields_frac_mask_collection);
  for (size_t i = 0; i < 2; ++i) {
    for (size_t c = 0; c < collection_size; ++c) {
      free(remote_src_buffer[i][c]);
    }
  }
  for (size_t c = 0; c < collection_size; ++c) {
    free(ref_tgt_field_collection[c]);
  }
  free(ref_tgt_field_collection);
}