Changed the cell-multipole interactions. The task now takes two cells and...

Changed the cell-multipole interactions. The task now takes two cells and interacts all the relevant siblings at once rather than in many separated tasks. The number of multipole tasks has been reduced dramatically, removing a lot of the overhead time associted with them.

examples/test_fmm.c

@@ -379,11 +379,44 @@ void cell_split(struct cell *c, struct qsched *s) {
 /* New tree walk */
 /* -------------------------------------------------------------------------- */
 
+
+/**
+ * @brief Checks whether the cells are direct neighbours ot not. Both cells have
+ * to be of the same size
+ *
+ */
+static inline int are_neighbours(struct cell *ci, struct cell *cj) {
+
+#ifdef SANITY_CHECKS
+  if (ci->h != cj->h)
+    error(" Cells of different size in distance calculation.");
+#endif
+
+  /* Maximum allowed distance */
+  float min_dist = ci->h;
+
+  /* (Manhattan) Distance between the cells */
+  for (int k = 0; k < 3; k++) {
+    float center_i = ci->loc[k];
+    float center_j = cj->loc[k];
+    if (fabsf(center_i - center_j) > min_dist) return 0;
+  }
+
+  return 1;
+}
+
+
+
+
 /**
- * @brief Interacts all particles in ci with the monopole in cj
+ * @brief Interacts all count particles in parts
+ * with the monopole in cj
+ *
+ * @param parts An array of particles
+ * @param count The number of particles in the array
+ * @param cj The cell with which the particles will interact
  */
-static inline void make_interact_pc(struct part *parts, int count,
-                                    double *loc, struct cell *cj) {
+static inline void make_interact_pc(struct part *parts, int count, struct cell *cj) {
 
   int j, k;
   float com[3] = {0.0, 0.0, 0.0};
@@ -429,42 +462,51 @@ static inline void make_interact_pc(struct part *parts, int count,
 }
 
 
+
+
 /**
- * @brief Checks whether the cells are direct neighbours ot not. Both cells have
- * to be of the same size
+ * @brief Loops all siblings of cells ci and cj and launches all
+ * the cell-monopole interactions that can be done between the siblings
+ *
+ * @param ci The first cell
+ * @param cj The second cell
  */
-static inline int are_neighbours(struct cell *ci, struct cell *cj) {
+static inline void iact_pair_pc(struct cell *ci, struct cell *cj) {
 
+  struct cell *cp, *cps;
+  
 #ifdef SANITY_CHECKS
   if (ci->h != cj->h)
-    error(" Cells of different size in distance calculation.");
+    error(" Cells of different size !.");
+
+  if ( ! are_neighbours(ci, cj) ) 
+     error(" Cells that are not neighbours !");
 #endif
 
-  /* Maximum allowed distance */
-  float min_dist = ci->h;
+  /* Let's split both cells and build all possible non-neighbouring pairs */
+  for (cp = ci->firstchild; cp != ci->sibling; cp = cp->sibling) {
+    for (cps = cj->firstchild; cps != cj->sibling; cps = cps->sibling) {
 
-  /* (Manhattan) Distance between the cells */
-  for (int k = 0; k < 3; k++) {
-    float center_i = ci->loc[k];
-    float center_j = cj->loc[k];
-    if (fabsf(center_i - center_j) > min_dist) return 0;
-  }
+      if ( ! are_neighbours(cp, cps) ) {
 
-  return 1;
-}
+	make_interact_pc(cp->parts,  cp->count, cps);
+	make_interact_pc(cps->parts, cps->count, cp );
 
+      }
+    }
+  }
+}
 
-static inline void iact_pair_pc(struct cell *ci, struct cell *cj) {
 
-#ifdef SANITY_CHECKS
-  if (ci->h != cj->h)
-    error(" Cells of different size in distance calculation.");
-#endif
 
-  make_interact_pc(ci->parts, ci->count, ci->loc, cj);
-  make_interact_pc(cj->parts, cj->count, cj->loc, ci);
 
-}
+/**
+ * @brief Interacts all particles in a cell ci with all particles in another
+ * cell cj using a simple double for-loop.
+ *
+ * @param ci The first cell
+ * @param cj The second cell
+ */
 
 static inline void iact_pair_direct(struct cell *ci, struct cell *cj) {
 
@@ -556,6 +598,8 @@ static inline void iact_pair_direct(struct cell *ci, struct cell *cj) {
   }
 }
 
+
+
 /**
  * @brief Decides whether two cells use the direct summation interaction or the
  * multipole interactions
@@ -599,17 +643,25 @@ void iact_pair(struct cell *ci, struct cell *cj) {
 	if (are_neighbours(cp, cps)) {
 	  iact_pair(cp, cps);
 	}
-	else { /* Otherwise do cell-mulitpole interactions */
-	  iact_pair_pc(cp, cps);
-	}
       }
     }
+    
+    /* Build all possible cell-multipole interaction pairs */
+    iact_pair_pc(ci, cj);
+
   } else { /* Otherwise, compute the interactions at this level directly. */
     iact_pair_direct(ci, cj);
   }
 }
 
 
+
+/**
+ * @brief Interacts all particles in cell c with themselves
+ * either by recursing or by direct summation
+ *
+ * @param c The #cell containing the particles
+ */
 void iact_self_direct(struct cell *c) {
   int i, j, k, count = c->count;
   float xi[3] = {0.0, 0.0, 0.0};
@@ -702,6 +754,7 @@ void iact_self_direct(struct cell *c) {
   } /* otherwise, compute interactions directly. */
 }
 
+
 /**
  * @brief Create the tasks for the cell pair/self.
  *
@@ -746,14 +799,13 @@ void create_tasks(struct qsched *s, struct cell *ci, struct cell *cj) {
       data[1] = NULL;
 
       /* Create the task. */
-      tid =
-          qsched_addtask(s, task_type_self, task_flag_none, data,
-                         sizeof(struct cell *) * 2, ci->count * ci->count / 2);
+      tid =  qsched_addtask(s, task_type_self, task_flag_none, data,
+			    sizeof(struct cell *) * 2, ci->count * ci->count / 2);
 
       /* Add the resource (i.e. the cell) to the new task. */
       qsched_addlock(s, tid, ci->res);
 
-/* If this call might recurse, add a dependency on the cell's COM
+/* Since this call might recurse, add a dependency on the cell's COM
    task. */
 #ifdef COM_AS_TASK
       if (ci->split) qsched_addunlock(s, ci->com_tid, tid);
@@ -763,55 +815,55 @@ void create_tasks(struct qsched *s, struct cell *ci, struct cell *cj) {
     /* Otherwise, it's a pair. */
   } else {
 
-    /* Are the cells NOT neighbours ? */
-    if (!are_neighbours(ci, cj)) {
-
-      /* We can do particle-monopole tasks */
+    /* Are both cells split and we are above the task limit ? */
+    if (ci->split && cj->split && ci->count > task_limit / cj->count) {
       
+      /* Let's split both cells and build all possible pairs */
+      for (cp = ci->firstchild; cp != ci->sibling; cp = cp->sibling) {
+	for (cps = cj->firstchild; cps != cj->sibling; cps = cps->sibling) {
+	  
+	  /* Recurse */
+	  if (are_neighbours(cp, cps)) {
+	    create_tasks(s, cp, cps);
+	  }
+	}
+      }
+
+      /* Let's also build a particle-monopole task */
+	
       /* Create the task. */
       data[0] = ci;
       data[1] = cj;
       tid = qsched_addtask(s, task_type_pair_pc, task_flag_none, data,
 			   sizeof(struct cell *) * 2, ci->count + cj->count);
-
+      
       /* Add the resource and dependance */
       qsched_addlock(s, tid, ci->res);
       qsched_addlock(s, tid, cj->res);
+
 #ifdef COM_AS_TASK
       qsched_addunlock(s, cj->com_tid, tid);
       qsched_addunlock(s, ci->com_tid, tid);
 #endif
 
+     
+    } else {  /* Otherwise, at least one of the cells is not split, build a direct
+	       * interaction. */
 
-    } else {/* Cells are direct neighbours */
-
-      /* Are both cells split and we are above the task limit ? */
-      if (ci->split && cj->split && ci->count > task_limit / cj->count) {
-
-        /* Let's split both cells and build all possible pairs */
-        for (cp = ci->firstchild; cp != ci->sibling; cp = cp->sibling) {
-          for (cps = cj->firstchild; cps != cj->sibling; cps = cps->sibling) {
-            /* Recurse */
-            create_tasks(s, cp, cps);
-          }
-        }
-        /* Otherwise, at least one of the cells is not split, build a direct
-         * interaction. */
-      } else {
-
-        /* Set the data. */
-        data[0] = ci;
-        data[1] = cj;
+      /* Set the data. */
+      data[0] = ci;
+      data[1] = cj;
+      
+      /* Create the task. */
+      tid = qsched_addtask(s, task_type_pair, task_flag_none, data,
+			   sizeof(struct cell *) * 2, ci->count * cj->count);
+      
+      /* Add the resources. */
+      qsched_addlock(s, tid, ci->res);
+      qsched_addlock(s, tid, cj->res);
 
-        /* Create the task. */
-        tid = qsched_addtask(s, task_type_pair, task_flag_none, data,
-                             sizeof(struct cell *) * 2, ci->count * cj->count);
+    }
 
-        /* Add the resources. */
-        qsched_addlock(s, tid, ci->res);
-        qsched_addlock(s, tid, cj->res);
-      }
-    } /* Cells are direct neighbours */
   }   /* Otherwise it's a pair */
 }