Add specific intersection function for sections/stripes.

src/Makefile.am

@@ -105,6 +105,7 @@ libyaxt_c_la_SOURCES =						\
 	xt_idxmod.c						\
 	xt_idxsection_internal.h				\
 	xt_idxsection.c						\
+	xt_idxsection_stripe_intersection.h			\
 	xt_idxstripes_internal.h				\
 	xt_idxstripes_pos_ext_map.h				\
 	xt_idxstripes.c						\

src/xt_idxlist_intersection.c

@@ -82,7 +82,7 @@ intersection_get_matrix[num_idxlist_classes][num_idxlist_classes] = {
   { empty_isect, xt_idxsection_get_intersection_with_other_idxlist,
     xt_idxsection_get_intersection_with_other_idxlist,
     xt_idxsection_get_intersection,
-    xt_idxsection_get_intersection_with_other_idxlist },
+    xt_idxsection_get_idxstripes_intersection },
   { empty_isect, xt_idxstripes_get_idxvec_intersection, xt_default_isect, xt_default_isect,
     xt_idxstripes_get_intersection },
 };

src/xt_idxsection.c

@@ -69,6 +69,7 @@
 #include "xt_config_internal.h"
 #include "xt_stripe_util.h"
 #include "xt_idxvec_internal.h"
+#include "xt_idxstripes_internal.h"
 
 static void
 idxsection_delete(Xt_idxlist data);
@@ -495,6 +496,41 @@ get_pos:;
   return result;
 }
 
+static void
+idxsection_init_sorted_copy(Xt_idxsection orig, Xt_idxsection copy);
+
+Xt_idxlist
+xt_idxsection_get_idxstripes_intersection(Xt_idxlist src_idxlist,
+                                          Xt_idxlist dst_idxlist,
+                                          Xt_config config)
+{
+  INSTR_DEF(instr,"idxsection_get_idxstripes_intersection")
+  INSTR_START(instr);
+  (void)config;
+  assert(dst_idxlist->vtable->idxlist_pack_code == STRIPES);
+  Xt_idxsection idxsection = (Xt_idxsection)src_idxlist;
+  Xt_idxlist idxstripes_list = dst_idxlist;
+#include "xt_idxsection_stripe_intersection.h"
+  INSTR_STOP(instr);
+  return result;
+}
+
+Xt_idxlist
+xt_idxsection_get_idxstripes_r_intersection(Xt_idxlist src_idxlist,
+                                            Xt_idxlist dst_idxlist,
+                                            Xt_config config)
+{
+  INSTR_DEF(instr,"idxsection_get_idxstripes_r_intersection")
+  INSTR_START(instr);
+  (void)config;
+  assert(src_idxlist->vtable->idxlist_pack_code == STRIPES);
+  Xt_idxsection idxsection = (Xt_idxsection)dst_idxlist;
+  Xt_idxlist idxstripes_list = src_idxlist;
+#include "xt_idxsection_stripe_intersection.h"
+  INSTR_STOP(instr);
+  return result;
+}
+
 Xt_idxlist
 xt_idxsection_get_intersection(Xt_idxlist idxlist_src,
                                Xt_idxlist idxlist_dst,

src/xt_idxsection_internal.h

@@ -67,6 +67,11 @@ PPM_DSO_INTERNAL Xt_idxlist
 xt_idxsection_get_intersection_with_other_idxlist(
   Xt_idxlist src_idxsection, Xt_idxlist dst_idxlist, Xt_config config);
 
+PPM_DSO_INTERNAL Xt_idxlist
+xt_idxsection_get_idxstripes_intersection(Xt_idxlist src_idxsection,
+                                          Xt_idxlist dst_idxlist,
+                                          Xt_config config);
+
 #endif
 
 /*

src/xt_idxsection_stripe_intersection.h

+/**
+ * @file xt_idxsection_stripe_intersection.h
+ * Contains code parts independent of whether section
+ * is destination or source side of intersection computation.
+ *
+ * @copyright Copyright  (C)  2016 Jörg Behrens <behrens@dkrz.de>
+ *                                 Moritz Hanke <hanke@dkrz.de>
+ *                                 Thomas Jahns <jahns@dkrz.de>
+ *
+ * @author Jörg Behrens <behrens@dkrz.de>
+ *         Moritz Hanke <hanke@dkrz.de>
+ *         Thomas Jahns <jahns@dkrz.de>
+ */
+/*
+ * Keywords:
+ * Maintainer: Jörg Behrens <behrens@dkrz.de>
+ *             Moritz Hanke <hanke@dkrz.de>
+ *             Thomas Jahns <jahns@dkrz.de>
+ * URL: https://dkrz-sw.gitlab-pages.dkrz.de/yaxt/
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are  permitted provided that the following conditions are
+ * met:
+ *
+ * Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ *
+ * Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * Neither the name of the DKRZ GmbH nor the names of its contributors
+ * may be used to endorse or promote products derived from this software
+ * without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+ * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+// intersection between an idxsection and a set of stripes
+size_t num_stripes = xt_idxstripes_get_num_index_stripes(idxstripes_list);
+
+const struct Xt_stripe *stripes
+  = xt_idxstripes_get_index_stripes_const(idxstripes_list);
+size_t num_dimensions = (size_t)idxsection->ndim;
+#ifndef __PGI
+struct {
+  struct Xt_idxsection_ h;
+  struct dim_desc d[num_dimensions];
+} sss;                        /* sorted section stack copy object */
+idxsection_init_sorted_copy(idxsection, &sss.h);
+struct dim_desc *dims = sss.d;
+#else
+struct Xt_idxsection_ *sssp
+= xmalloc(sizeof (struct Xt_idxsection_)
+          + sizeof (struct dim_desc) * num_dimensions);
+idxsection_init_sorted_copy(idxsection, sssp);
+struct dim_desc *dims = sssp->dims;
+#endif
+size_t num_result_stripes_total = 0;
+{
+  /* start with non-appendable terminator */
+  Xt_int accum_start = XT_INT_MIN;
+  int buf_nstrides = 0;
+  Xt_int buf_stride = -1;
+  Xt_int index_continuation = XT_INT_MIN-1;
+  for (size_t k = 0; k < num_stripes; ++k) {
+    struct Xt_stripe query = stripes[k];
+    size_t num_result_stripes = 0;
+    for (int j = 0; j < query.nstrides; ++j) {
+      Xt_int insert_index = (Xt_int)(query.start + query.stride * j),
+        running_index = insert_index;
+      bool out_of_bounds = false;
+      for (size_t i = 0; i < num_dimensions; ++i) {
+        Xt_int abs_global_stride = dims[i].global_stride;
+        Xt_int curr_global_position
+          = (Xt_int)(running_index / abs_global_stride);
+        running_index = (Xt_int)(running_index % abs_global_stride);
+        out_of_bounds |= (curr_global_position < dims[i].local_start)
+          | (curr_global_position >= dims[i].local_start + dims[i].local_size);
+      }
+      if (!out_of_bounds) {
+        if (insert_index == index_continuation || buf_nstrides == 1) {
+          /* already existing stripe can be continued */
+          ++buf_nstrides;
+          if (insert_index != index_continuation)
+            buf_stride = (Xt_int)(insert_index - accum_start);
+        } else {
+          buf_stride = 1;
+          buf_nstrides = 1;
+          accum_start = insert_index;
+          ++num_result_stripes;
+        }
+        index_continuation = (Xt_int)(insert_index + buf_stride);
+      }
+    }
+    num_result_stripes_total += num_result_stripes;
+  }
+}
+Xt_idxlist result;
+if (num_result_stripes_total) {
+  struct Xt_stripes_alloc result_alloc
+    = xt_idxstripes_alloc(num_result_stripes_total);
+  struct Xt_stripe *restrict result_stripes = result_alloc.stripes;
+  size_t num_result_stripes = (size_t)-1;
+  Xt_int index_continuation = XT_INT_MIN-1;
+  int buf_nstrides = 0;
+  Xt_int buf_stride = -1;
+  for (size_t k = 0; k < num_stripes; ++k) {
+    struct Xt_stripe query = stripes[k];
+    for (int j = 0; j < query.nstrides; ++j) {
+      Xt_int insert_index = (Xt_int)(query.start + query.stride * j),
+        running_index = insert_index;
+      bool out_of_bounds = false;
+      for (size_t i = 0; i < num_dimensions; ++i) {
+        Xt_int abs_global_stride = dims[i].global_stride;
+        Xt_int curr_global_position
+          = (Xt_int)(running_index / abs_global_stride);
+        running_index = (Xt_int)(running_index % abs_global_stride);
+        out_of_bounds |= (curr_global_position < dims[i].local_start)
+          | (curr_global_position >= dims[i].local_start + dims[i].local_size);
+      }
+      if (!out_of_bounds) {
+        if (insert_index == index_continuation || buf_nstrides == 1) {
+          /* already existing stripe can be expanded */
+          ++buf_nstrides;
+          if (insert_index != index_continuation)
+            buf_stride = (Xt_int)(insert_index
+                                  - result_stripes[num_result_stripes].start);
+        } else {
+          result_stripes[num_result_stripes].nstrides = buf_nstrides;
+          result_stripes[num_result_stripes].stride = buf_stride;
+          buf_stride = 1;
+          buf_nstrides = 1;
+          result_stripes[++num_result_stripes].start = insert_index;
+        }
+        index_continuation = (Xt_int)(insert_index + buf_stride);
+      }
+    }
+  }
+  result_stripes[num_result_stripes].stride = buf_stride;
+  result_stripes[num_result_stripes].nstrides = buf_nstrides;
+  assert(num_result_stripes_total == num_result_stripes+1);
+  Xt_idxlist initial_result = xt_idxstripes_congeal(result_alloc);
+  if (XT_CONFIG_GET_FORCE_NOSORT(config)
+      || xt_idxlist_get_sorting(initial_result) == 1)
+    result = initial_result;
+  else {
+    result = xt_idxlist_sorted_copy(initial_result);
+    xt_idxlist_delete(initial_result);
+  }
+} else
+  result = xt_idxempty_new();
+#ifdef __PGI
+free(sssp);
+#endif
+
+/*
+ * Local Variables:
+ * c-basic-offset: 2
+ * coding: utf-8
+ * indent-tabs-mode: nil
+ * show-trailing-whitespace: t
+ * require-trailing-newline: t
+ * End:
+ */

tests/test_idxsection.c

@@ -558,6 +558,44 @@ int main(int argc, char **argv) {
     xt_idxlist_delete(idxsection);
   }
 
+  {
+    // check idxsection_get_idxstripes_intersection
+
+    Xt_int start = 0;
+    enum { num_dimensions = 2 };
+    static const Xt_int global_size[num_dimensions] = {4,4};
+    static const int local_size [num_dimensions] = {2,2};
+    static const Xt_int local_start[num_dimensions] = {0,2};
+
+    Xt_idxlist idxsection
+      = xt_idxsection_new(start, num_dimensions, global_size,
+                          local_size, local_start);
+
+    // we have indices = {2,3,6,7}
+    static const struct Xt_stripe stripes[]
+      = { { .start = 2, .stride = -1, .nstrides = 2 },
+          { .start = 5, .stride = 2, .nstrides = 2 },
+          { .start = 6, .stride = 1, .nstrides = 2 },
+          { .start = 7, .stride = -1, .nstrides = 2 },
+          { .start = 8, .stride = 1, .nstrides = 1 },
+    };
+    enum { num_stripes = sizeof(stripes) / sizeof(stripes[0]) };
+
+    Xt_idxlist sel_idxlist = xt_idxstripes_new(stripes, num_stripes);
+
+    Xt_idxlist inter_idxlist
+      = xt_idxlist_get_intersection(idxsection, sel_idxlist);
+
+    static const Xt_int ref_inter_idx[] = {2,6,6,7,7,7};
+    enum { num_ref_inter_idx = sizeof(ref_inter_idx) / sizeof(*ref_inter_idx) };
+
+    check_idxlist(inter_idxlist, ref_inter_idx, num_ref_inter_idx);
+
+    xt_idxlist_delete(inter_idxlist);
+    xt_idxlist_delete(sel_idxlist);
+    xt_idxlist_delete(idxsection);
+  }
+
   { // test 2D section with negative global size
     // iterate through all sign combinations of -/+ for local and global
     // and for x and y, giving 2^2^2 combinations