test_interp_operator_sum_mvp_at_src.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_src.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_src(struct yac_collection_selection * sel,
                                       Xt_redist * src_halo_redists,
                                       Xt_redist tgt_result_redist,
                                       size_t tgt_count,
                                       size_t * num_src_points_per_field,
                                       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 = 4, .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 local target points that are to be interpolated
        //----------------------------------------------------------------------
 
        struct {
 
          Xt_int * ids;  
          int * offsets; 
          size_t count;  
        } tgt_receive;
 
        //----------------------------------------------------------------------
        // Information about the target points that are to be interpolated
        // (these target points can either be owned by the local process or
        //  remote process
        //----------------------------------------------------------------------
 
        struct {
 
          //--------------------------------------------------------------------
          // Information about the target points that are to be interpolated by
          // the local process
          //--------------------------------------------------------------------
 
          Xt_int * ids; 
          size_t count; 
 
          //--------------------------------------------------------------------
          // Information about the source points from remote processes required
          // the interpolation of the 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 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
 
        } tgt_send;
 
        //----------------------------------------------------------------------
        // 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 =
          {{.tgt_receive =
              {.ids = NULL,
               .offsets = NULL,
               .count = 0},
            .tgt_send = {
              .ids = NULL,
              .count = 0,
              .src_global_ids_per_field = {NULL},
              .num_src_points_per_field = {0},
              .num_src_per_tgt = NULL,
              .weights = NULL,
              .src_field_idx = NULL,
              .src_idx = NULL},
            .is_source = 0,
            .is_target = 0},
           {.tgt_receive =
              {.ids = NULL,
               .offsets = NULL,
               .count = 0},
            .tgt_send = {
              .ids = NULL,
              .count = 0,
              .src_global_ids_per_field = {NULL},
              .num_src_points_per_field = {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 on rank 1 contains the source points of rank 0 of
       // source field 1
       // dimensions (for documentation only):
       //   tgt[rank][local_idx]
       //   src[rank][src_field_idx][local_idx]
       //
       // compute rank 0:
       //   tgt[1][0] = src[0][1][0]
       //   tgt[1][1] = src[0][1][1]
       //   tgt[1][2] = src[0][1][2]
       //   tgt[1][3] = src[0][1][3]
       {.interp_data =
          {{.tgt_receive =
              {.ids = NULL,
               .offsets = NULL,
               .count = 0},
            .tgt_send = {
              .ids = (Xt_int[]){4,5,6,7},
              .count = 4,
              .src_global_ids_per_field = {NULL,NULL,NULL},
              .num_src_points_per_field = {0,0,0},
              .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[]){1,1,1,1},
              .src_idx = (size_t[]){0,1,2,3}},
            .is_source = 1,
            .is_target = 0},
           {.tgt_receive =
              {.ids = (Xt_int[]){4,5,6,7},
               .offsets = (int[]){0,1,2,3},
               .count = 4},
            .tgt_send = {
              .ids = NULL,
              .count = 0,
              .src_global_ids_per_field = {NULL,NULL,NULL},
              .num_src_points_per_field = {0,0,0},
              .num_src_per_tgt = NULL,
              .weights = NULL,
              .src_field_idx = NULL,
              .src_idx = NULL},
            .is_source = 0,
            .is_target = 1}},
        .num_src_fields = 3},
       // some target data on rank 1 contains a weighted sum of source points
       // of rank 0 of all source fields
       // dimensions (for documentation only):
       //   tgt[rank][local_idx]
       //   src[rank][src_field_idx][local_idx]
       //
       // compute rank 0:
       //   tgt[1][1] = 0.4 * src[0][0][1] + 0.3 * src[0][1][1] + 0.3 * src[0][2][1]
       //   tgt[1][3] = 0.4 * src[0][0][3] + 0.3 * src[0][1][3] + 0.3 * src[0][2][3]
       {.interp_data =
          {{.tgt_receive =
              {.ids = NULL,
               .offsets = NULL,
               .count = 0},
            .tgt_send = {
              .ids = (Xt_int[]){5,7},
              .count = 2,
              .src_global_ids_per_field = {NULL,NULL,NULL},
              .num_src_points_per_field = {0,0,0},
              .num_src_per_tgt = (size_t[]){3,3},
              .weights = (double[]){0.4,0.3,0.3, 0.4,0.3,0.3},
              .src_field_idx = (size_t[]){0,1,2, 0,1,2},
              .src_idx = (size_t[]){1,1,1, 3,3,3}},
            .is_source = 1,
            .is_target = 0},
           {.tgt_receive =
              {.ids = (Xt_int[]){5,7},
               .offsets = (int[]){1,3},
               .count = 2},
            .tgt_send = {
              .ids = NULL,
              .count = 0,
              .src_global_ids_per_field = {NULL,NULL,NULL},
              .num_src_points_per_field = {0,0,0},
              .num_src_per_tgt = NULL,
              .weights = NULL,
              .src_field_idx = NULL,
              .src_idx = NULL},
            .is_source = 0,
            .is_target = 1}},
        .num_src_fields = 3},
       // some target data on rank 0 contains a weighted sum of source points
       // of rank 0 of all source fields
       // dimensions (for documentation only):
       //   tgt[rank][local_idx]
       //   src[rank][src_field_idx][local_idx]
       //
       // compute rank 0:
       //   tgt[0][1] = 0.4 * src[0][0][1] + 0.3 * src[0][1][1] + 0.3 * src[0][2][1]
       //   tgt[0][3] = 0.4 * src[0][0][3] + 0.3 * src[0][1][3] + 0.3 * src[0][2][3]
       {.interp_data =
          {{.tgt_receive =
              {.ids = (Xt_int[]){1,3},
               .offsets = (int[]){1,3},
               .count = 2},
            .tgt_send = {
              .ids = (Xt_int[]){1,3},
              .count = 2,
              .src_global_ids_per_field = {NULL,NULL,NULL},
              .num_src_points_per_field = {0,0,0},
              .num_src_per_tgt = (size_t[]){3,3},
              .weights = (double[]){0.4,0.3,0.3, 0.4,0.3,0.3},
              .src_field_idx = (size_t[]){0,1,2, 0,1,2},
              .src_idx = (size_t[]){1,1,1, 3,3,3}},
            .is_source = 1,
            .is_target = 1},
           {.tgt_receive =
              {.ids = NULL,
               .offsets = NULL,
               .count = 0},
            .tgt_send = {
              .ids = NULL,
              .count = 0,
              .src_global_ids_per_field = {NULL,NULL,NULL},
              .num_src_points_per_field = {0,0,0},
              .num_src_per_tgt = NULL,
              .weights = NULL,
              .src_field_idx = NULL,
              .src_idx = NULL},
            .is_source = 0,
            .is_target = 0}},
        .num_src_fields = 3},
       // some target data on rank 1 contains a weighted sum of various source
       // points from both ranks
       // dimensions (for documentation only):
       //   tgt[rank][local_idx]
       //   src[rank][src_field_idx][local_idx]
       //
       // compute rank 0:
       //   tgt[1][0] = 0.1 * src[0][1][0] + 0.1 * src[0][1][1] +
       //               0.1 * src[1][1][0] + 0.1 * src[1][1][1]
       //   tgt[1][2] = 0.2 * src[0][0][2] + 0.2 * src[0][2][2] +
       //               0.2 * src[1][2][2]
       {.interp_data =
          {{.tgt_receive =
              {.ids = NULL,
               .offsets = NULL,
               .count = 0},
            .tgt_send = {
              .ids = (Xt_int[]){4,6},
              .count = 2,
              .src_global_ids_per_field = {NULL,
                                           (Xt_int[]){4,5},
                                           (Xt_int[]){6}},
              .num_src_points_per_field = {0,2,1},
              .num_src_per_tgt = (size_t[]){4,3},
              .weights = (double[]){0.1,0.1,0.1,0.1, 0.2,0.2,0.2},
              .src_field_idx = (size_t[]){1,1,SIZE_MAX,SIZE_MAX,
                                          0,2,SIZE_MAX},
              .src_idx = (size_t[]){0,1,0,1, 2,2,2}},
            .is_source = 1,
            .is_target = 0},
           {.tgt_receive =
              {.ids = (Xt_int[]){4,6},
               .offsets = (int[]){0,2},
               .count = 2},
            .tgt_send = {
              .ids = NULL,
              .count = 0,
              .src_global_ids_per_field = {NULL,NULL,NULL},
              .num_src_points_per_field = {0,0,0},
              .num_src_per_tgt = NULL,
              .weights = NULL,
              .src_field_idx = NULL,
              .src_idx = NULL},
            .is_source = 1,
            .is_target = 1}},
        .num_src_fields = 3},
       // complex interpolation involving all processes
       // dimensions (for documentation only):
       //   tgt[rank][local_idx]
       //   src[rank][src_field_idx][local_idx]
       //
       // compute rank 0:
       //   tgt[0][1] = 0.1 * src[0][0][0] + 0.2 * src[0][0][1] +
       //               0.3 * src[1][0][2] + 0.4 * src[1][0][3]
       //   tgt[0][2] = 0.1 * src[0][1][0] + 0.2 * src[0][1][1]
       //   tgt[0][0] = 0.2 * src[1][2][2]
       //   tgt[1][0] = 0.1 * src[0][1][0] + 0.2 * src[0][1][1] +
       //               0.3 * src[1][1][2] + 0.4 * src[1][1][3]
       //   tgt[1][1] = 0.1 * src[0][1][0] + 0.2 * src[0][1][1]
       // compute rank 1:
       //   tgt[0][3] = 0.25 * src[1][2][0] + 0.25 * src[1][2][1] +
       //               0.25 * src[1][2][2] + 0.25 * src[1][2][3]
       //   tgt[1][2] = 0.1 * src[0][0][0] + 0.2 * src[0][0][1] +
       //               0.3 * src[1][0][2] + 0.4 * src[1][0][3]
       //   tgt[1][3] = 0.1 * src[0][1][0] + 0.2 * src[0][1][1]
       //
       {.interp_data =
          {{.tgt_receive =
              {.ids = (Xt_int[]){0,1,2,3},
               .offsets = (int[]){0,1,2,3},
               .count = 4},
            .tgt_send = {
              .ids = (Xt_int[]){1,2,0,4,5},
              .count = 5,
              .src_global_ids_per_field = {(Xt_int[]){6,7},
                                           (Xt_int[]){6,7},
                                           (Xt_int[]){6}},
              .num_src_points_per_field = {2,2,1},
              .num_src_per_tgt = (size_t[]){4,2,1,4,2},
              .weights = (double[]){0.1,0.2,0.3,0.4,
                                    0.1,0.2,
                                    0.2,
                                    0.1,0.2,0.3,0.4,
                                    0.1,0.2},
              .src_field_idx = (size_t[]){0,0,SIZE_MAX,SIZE_MAX,
                                          1,1,
                                          SIZE_MAX,
                                          1,1,SIZE_MAX,SIZE_MAX,
                                          1,1},
              .src_idx = (size_t[]){0,1,0,1,
                                    0,1,
                                    4,
                                    0,1,2,3,
                                    0,1}},
            .is_source = 1,
            .is_target = 1},
           {.tgt_receive =
              {.ids = (Xt_int[]){4,5,6,7},
               .offsets = (int[]){0,1,2,3},
               .count = 4},
            .tgt_send = {
              .ids = (Xt_int[]){3,6,7},
              .count = 3,
              .src_global_ids_per_field = {(Xt_int[]){0,1},
                                           (Xt_int[]){0,1},
                                           NULL},
              .num_src_points_per_field = {2,2,0},
              .num_src_per_tgt = (size_t[]){4,4,2},
              .weights = (double[]){0.25,0.25,0.25,0.25,
                                    0.1,0.2,0.3,0.4,
                                    0.1,0.2},
              .src_field_idx = (size_t[]){2,2,2,2,
                                          SIZE_MAX,SIZE_MAX,0,0,
                                          SIZE_MAX,SIZE_MAX},
              .src_idx = (size_t[]){0,1,2,3,
                                    0,1,2,3,
                                    2,3}},
            .is_source = 1,
            .is_target = 1}},
        .num_src_fields = 3}};
    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_interp_operator_sum_mvp_at_src_new from list
      // of links
      size_t num_src_fields =
        interp_configs[interp_config_idx].num_src_fields;
      Xt_redist src_halo_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 target points, which are to be
        // interpolated by the local process
        Xt_idxlist required_source_idxlist =
          (interp_configs[interp_config_idx].
            interp_data[comm_rank].
              tgt_send.num_src_points_per_field[src_field_idx] > 0)?
            xt_idxvec_new(
              interp_configs[interp_config_idx].
                interp_data[comm_rank].
                  tgt_send.src_global_ids_per_field[src_field_idx],
              interp_configs[interp_config_idx].
                interp_data[comm_rank].
                  tgt_send.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_halo_redists[src_field_idx] = xt_redist_p2p_new(xmap, MPI_DOUBLE);
        xt_xmap_delete(xmap);
        xt_idxlist_delete(required_source_idxlist);
 
      } // src_field_idx
 
      // create redist for receiving local target points that were
      // interpolated on the source processes
      Xt_redist tgt_result_redist;
      {
        Xt_idxlist interpolated_target_idxlist =
          (interp_configs[interp_config_idx].
            interp_data[comm_rank].tgt_send.count > 0)?
            xt_idxvec_new(
              interp_configs[interp_config_idx].
                interp_data[comm_rank].tgt_send.ids,
              interp_configs[interp_config_idx].
                interp_data[comm_rank].tgt_send.count):
            xt_idxempty_new();
        int * interpolated_target_offsets =
          xmalloc(
            interp_configs[interp_config_idx].
              interp_data[comm_rank].tgt_send.count *
            sizeof(*interpolated_target_offsets));
        for (size_t t = 0;
             t < interp_configs[interp_config_idx].
                   interp_data[comm_rank].tgt_send.count; ++t) {
          interpolated_target_offsets[t] = (int)t;
        }
        Xt_idxlist required_target_idxlist =
          (interp_configs[interp_config_idx].
            interp_data[comm_rank].tgt_receive.count > 0)?
            xt_idxvec_new(
              interp_configs[interp_config_idx].
                interp_data[comm_rank].tgt_receive.ids,
              interp_configs[interp_config_idx].
                interp_data[comm_rank].tgt_receive.count):
            xt_idxempty_new();
 
        Xt_xmap xmap =
          xt_xmap_all2all_new(
            interpolated_target_idxlist, required_target_idxlist,
            MPI_COMM_WORLD);
        tgt_result_redist =
          xt_redist_p2p_off_new(
            xmap, interpolated_target_offsets,
            interp_configs[interp_config_idx].
                interp_data[comm_rank].tgt_receive.offsets, MPI_DOUBLE);
 
        xt_xmap_delete(xmap);
        xt_idxlist_delete(required_target_idxlist);
        xt_idxlist_delete(interpolated_target_idxlist);
        free(interpolated_target_offsets);
      }
 
      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].tgt_send.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_send.count; ++i) {
              total_num_receive_src_points +=
                interp_configs[interp_config_idx].
                  interp_data[comm_rank].tgt_send.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].
                    tgt_send.src_field_idx[i] == SIZE_MAX) ||
                (interp_configs[interp_config_idx].
                   interp_data[comm_rank].
                     tgt_send.src_field_idx[i] < num_src_fields),
                "ERROR in src_field_idx");
              YAC_ASSERT(
                interp_configs[interp_config_idx].
                  interp_data[comm_rank].tgt_send.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_src_new(
              sel, src_halo_redists,
              interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.count,
              interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_per_tgt,
              weights,
              interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.src_field_idx,
              interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.src_idx,
              num_src_fields, tgt_result_redist, 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_src(
                sel, src_halo_redists, tgt_result_redist,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.count,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.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_src(
                sel, src_halo_redists, tgt_result_redist,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.count,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.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_src(
                sel, src_halo_redists, tgt_result_redist,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.count,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.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_src(
                sel, src_halo_redists, tgt_result_redist,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.count,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.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_src(
                sel, src_halo_redists, tgt_result_redist,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.count,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_points_per_field,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.num_src_per_tgt,
                weights,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.src_field_idx,
                interp_configs[interp_config_idx].interp_data[comm_rank].tgt_send.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
 
      xt_redist_delete(tgt_result_redist);
      for (size_t f = 0; f < num_src_fields; ++f)
        xt_redist_delete(src_halo_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_src(struct yac_collection_selection * sel,
                                       Xt_redist * src_halo_redists,
                                       Xt_redist tgt_result_redist,
                                       size_t tgt_count,
                                       size_t * num_src_points_per_field,
                                       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_wgt
  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 temporary target field collection buffer, which will hold
  // interpolation results on the source process. These results are then
  // redistributed to the target processes using tgt_result_redist.
  double ** temp_tgt_field_collection =
    xmalloc(collection_size * sizeof(*temp_tgt_field_collection));
  for (size_t c = 0; c < collection_size; ++c) {
    temp_tgt_field_collection[c] =
      xmalloc(tgt_count * sizeof(**temp_tgt_field_collection));
    for (size_t i = 0; i < tgt_count; ++i) {
      temp_tgt_field_collection[c][i] = TGT_UNSET_VALUE;
    }
  }
 
  // On the source processes the interpolated target points are written into
  // a contiguous.
  size_t * tgt_pos = xmalloc(tgt_count * sizeof(*tgt_pos));
  for (size_t t = 0; t < tgt_count; ++t) {tgt_pos[t] = t;}
 
  // 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_wgt 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_halo_redists[f],
        (const void *)(src_fields_collection[c_idx][f]),
        (void *)(remote_src_fields[c][f]));
      xt_redist_s_exchange1(
        src_halo_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_wgt(
    (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,
    temp_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);
 
  // for each collection entry
  for (size_t c = 0; c < collection_size; ++c) {
 
    // redistribute the interpolated target points from the source processes to
    // the target ones
    xt_redist_s_exchange1(
      tgt_result_redist,
      (const void *)(temp_tgt_field_collection[c]),
      (void *)(ref_tgt_field_collection[c]));
  }
 
  // 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_src 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);
  for (size_t c = 0; c < collection_size; ++c) {
    free(temp_tgt_field_collection[c]);
  }
  free(temp_tgt_field_collection);
  free(tgt_pos);
}