adds alternative implementation for routine yac_proc_sphere_part_do_point_search

* this should espectially work well on vector machines
* it is enabled by compiling with "-DYAC_NEC_EXPERIMENTAL"

src/core/proc_sphere_part.c

@@ -632,6 +632,51 @@ static int is_serial_node(struct proc_sphere_part_node * node) {
          (node->U.data.rank == 0) && (node->T.data.rank == 0);
 }
 
+// #define YAC_NEC_EXPERIMENTAL
+#ifdef YAC_NEC_EXPERIMENTAL
+// the following code may be better for a vector machine
+static void yac_proc_sphere_part_do_point_search_recursive(
+  struct proc_sphere_part_node * node, yac_coordinate_pointer search_coords,
+  size_t * search_idx, size_t * temp_search_idx, int * flag,
+  size_t count, int * ranks) {
+
+  for (size_t i = 0; i < count; ++i) {
+
+    // compute cos angle between current search point and norm vector
+    double dot = search_coords[search_idx[i]][0] * node->gc_norm_vector[0] +
+                 search_coords[search_idx[i]][1] * node->gc_norm_vector[1] +
+                 search_coords[search_idx[i]][2] * node->gc_norm_vector[2];
+
+    // if (angle >= M_PI_2)
+    flag[i] = (dot <= 0.0);
+  }
+
+  size_t u_size = 0, t_size = 0;
+  for (size_t i = 0; i < count; ++i) {
+    if (flag[i]) search_idx[u_size++] = search_idx[i];
+    else         temp_search_idx[t_size++] = search_idx[i];
+  }
+
+  if (node->U.is_leaf) {
+    size_t rank = node->U.data.rank;;
+    for (size_t i = 0; i < u_size; ++i) ranks[search_idx[i]] = rank;
+  } else {
+    yac_proc_sphere_part_do_point_search_recursive(
+      node->U.data.node, search_coords, search_idx, temp_search_idx + t_size,
+      flag, u_size, ranks);
+  }
+
+  if (node->T.is_leaf) {
+    size_t rank = node->T.data.rank;
+    for (size_t i = 0; i < t_size; ++i) ranks[temp_search_idx[i]] = rank;
+  } else {
+    yac_proc_sphere_part_do_point_search_recursive(
+      node->T.data.node, search_coords, temp_search_idx, search_idx + u_size,
+      flag, t_size, ranks);
+  }
+}
+#endif // YAC_NEC_EXPERIMENTAL
+
 void yac_proc_sphere_part_do_point_search(
   struct proc_sphere_part_node * node, yac_coordinate_pointer search_coords,
   size_t count, int * ranks) {
@@ -641,7 +686,17 @@ void yac_proc_sphere_part_do_point_search(
     return;
   }
 
+#ifdef YAC_NEC_EXPERIMENTAL
+
+  size_t * search_idx = xmalloc(2 * count * sizeof(*search_idx));
+  for (size_t i = 0; i < count; ++i) search_idx[i] = i;
+  int * flag = xmalloc(count * sizeof(*flag));
+  yac_proc_sphere_part_do_point_search_recursive(
+    node, search_coords, search_idx, search_idx + count, flag, count, ranks);
+  free(flag);
+  free(search_idx);
 
+#else
 
   YAC_OMP_PARALLEL
   {
@@ -681,6 +736,7 @@ void yac_proc_sphere_part_do_point_search(
       }
     }
   }
+#endif // YAC_NEC_EXPERIMENTAL
 }
 
 static void bnd_circle_search(