From 2f8f8c582c00f78c6fe1634c315f507ebf4b5619 Mon Sep 17 00:00:00 2001
From: Yohai Meiron <sahmes@gmail.com>
Date: Fri, 1 May 2020 00:27:41 -0400
Subject: [PATCH] Attempting to get fast active search working

---
 io.h         |   2 +-
 phigrape.cpp | 155 +++++++++++++++++++++++++++++++--------------------
 2 files changed, 97 insertions(+), 60 deletions(-)

diff --git a/io.h b/io.h
index 1baa152..770a0f5 100644
--- a/io.h
+++ b/io.h
@@ -5,7 +5,7 @@
 bool is_hdf5(std::string file_name);
 // This function is implemented independently of the HDF5 library
 
-void ascii_read(const std::string file_name, int *step_num, int *N, double *t, double *m, double3 *x, double3 *v);
+void ascii_read(const std::string file_name, int& step_num, int& N, double& t, double **m, double3 **x, double3 **v);
 
 void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, int precision=10);
 
diff --git a/phigrape.cpp b/phigrape.cpp
index 69f6d93..8d176d8 100644
--- a/phigrape.cpp
+++ b/phigrape.cpp
@@ -82,7 +82,9 @@ Last redaction : 2019.04.16 12:55
 
 #ifdef ETICS
 #include "grapite.h"
+//#define ACT_DEF_GRAPITE
 #endif
+const bool act_def_grapite = true;
 
 Config *config;
 
@@ -248,42 +250,44 @@ public:
     double get_minimum_time(const double t[], const double dt[])
     {
         double min_t_loc, min_t;
-#ifdef ACT_DEF_GRAPITE
-        min_t_loc = grapite_get_minimum_time();
-#else
-        min_t_loc = t[myRank*n_loc]+dt[myRank*n_loc];
-        for (int j=myRank*n_loc+1; j<(myRank+1)*n_loc; j++) {
-            double tmp = t[j] + dt[j];
-            if (tmp < min_t_loc) min_t_loc = tmp;
+        if (act_def_grapite) {
+            min_t_loc = grapite_get_minimum_time();
+            printf("gggggggggggg  min_t_loc=%.10e\n", min_t_loc);
+        } else {
+            min_t_loc = t[myRank*n_loc]+dt[myRank*n_loc];
+            for (int j=myRank*n_loc+1; j<(myRank+1)*n_loc; j++) {
+                double tmp = t[j] + dt[j];
+                if (tmp < min_t_loc) min_t_loc = tmp;
+            }
         }
-#endif
         /* Reduce the "global" min_t from min_t_loc "local" on all processors) */
         MPI_Allreduce(&min_t_loc, &min_t, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
         return min_t;
     }
     void get_active_indices(const double min_t, const double t[], const double dt[], int ind_act[], int *n_act)
+#warning refrence not pointer
     {
-#ifdef ACT_DEF_GRAPITE
-        int n_act_loc;
-        grapite_active_search(min_t, ind_act_loc, &n_act_loc);
-        if (myRank > 0)
-            for (int i=0; i<n_act_loc; i++)
-                ind_act_loc[i] += myRank*n_loc;
-        int n_act_arr[256], displs[256]; // Assuming maximum of 256 processes... seems safe.
-        MPI_Allgather(&n_act_loc, 1, MPI_INT, n_act_arr, 1, MPI_INT, MPI_COMM_WORLD);
-        *n_act = n_act_arr[0];
-        for (int i=1; i<n_proc; i++)
-            *n_act += n_act_arr[i];
-        displs[0] = 0;
-        for (int i=1; i<n_proc; i++)
-            displs[i]=displs[i-1]+n_act_arr[i-1];
-        MPI_Allgatherv(ind_act_loc, n_act_loc, MPI_INT, ind_act, n_act_arr, displs, MPI_INT, MPI_COMM_WORLD);
-#else
-        *n_act = 0;
-        for (int i=0; i<N; i++) {
-            if (t[i]+dt[i] == min_t) ind_act[(*n_act)++] = i;
-        } /* i */
-#endif // ACT_DEF_GRAPITE
+        if (act_def_grapite) {
+            int n_act_loc;
+            grapite_active_search(min_t, ind_act_loc, &n_act_loc);
+            if (myRank > 0)
+                for (int i=0; i<n_act_loc; i++)
+                    ind_act_loc[i] += myRank*n_loc;
+            int n_act_arr[256], displs[256]; // Assuming maximum of 256 processes... seems safe.
+            MPI_Allgather(&n_act_loc, 1, MPI_INT, n_act_arr, 1, MPI_INT, MPI_COMM_WORLD);
+            *n_act = n_act_arr[0];
+            for (int i=1; i<n_proc; i++)
+                *n_act += n_act_arr[i];
+            displs[0] = 0;
+            for (int i=1; i<n_proc; i++)
+                displs[i]=displs[i-1]+n_act_arr[i-1];
+            MPI_Allgatherv(ind_act_loc, n_act_loc, MPI_INT, ind_act, n_act_arr, displs, MPI_INT, MPI_COMM_WORLD);
+        } else {
+            *n_act = 0;
+            for (int i=0; i<N; i++) {
+                if (t[i]+dt[i] == min_t) ind_act[(*n_act)++] = i;
+            } /* i */
+        }
     }
 private:
     int myRank, n_proc, n_loc, N;
@@ -670,6 +674,7 @@ int main(int argc, char *argv[])
 #endif
 
         double min_t = active_search.get_minimum_time(t, dt);
+        printf("zzzzzzzzzzzzzzzzz %.10e\n", min_t);
 
 #ifdef TIMING
         get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
@@ -683,39 +688,71 @@ int main(int argc, char *argv[])
         active_search.get_active_indices(min_t, t, dt, ind_act, &n_act);
 
 // TODO deal with it below
-#ifdef ACT_DEF_GRAPITE
-#error please fix here
-#endif
-#if defined(ACT_DEF_GRAPITE) && (defined(ADD_BH1) || defined(ADD_BH2))
-#ifdef ADD_BH1
-#define ACT_DEF_GRAPITE_NUMBH 1
-#else
-#define ACT_DEF_GRAPITE_NUMBH 2
-#endif
-        int acf_def_grapite_bh_count = 0;
-        for (i=0; i<n_act; i++) {
-            if (ind_act[i]==0) {
-                i_bh1 = i;
-                acf_def_grapite_bh_count++;
-        }
-#ifdef ADD_BH2
-            else if (ind_act[i]==1) {
-                i_bh2 = i;
-                acf_def_grapite_bh_count++;
+// #ifdef ACT_DEF_GRAPITE
+// #error please fix here
+// #endif
+
+        static int printouts = 0;
+        int n_bh = config->live_smbh_count;
+        if (n_bh>0) {
+            if (act_def_grapite) {
+                int act_def_grapite_bh_count = 0;
+                int i_bh[n_bh];
+                for (int i=0; i<n_act; i++) {
+                    if (ind_act[i]<n_bh) {
+                        i_bh[ind_act[i]] = i;
+                        if (++act_def_grapite_bh_count == n_bh) break;
+                    }
+                }
+                i_bh1 = i_bh[0];
+                if (n_bh == 2) i_bh2 = i_bh[1];
+            } else {
+                i_bh1 == 0;
+                if (n_bh == 2) i_bh2 = 1;
             }
-#endif
-            if (acf_def_grapite_bh_count==ACT_DEF_GRAPITE_NUMBH) break;
         }
-        if (i==n_act) {
-            fprintf(stderr, "ERROR: black holes were not found in the active particle list");
-            exit(1);
+        printf("previously got i_bh1=%d and i_bh2=%d\n", i_bh1, i_bh2);
+        for (int i=0; i<n_act; i++) {
+            if (ind_act[i]==0) i_bh1=i;
+            if (ind_act[i]==1) i_bh2=i;
         }
-#else
-        if (config->live_smbh_count > 0) {
-            i_bh1 = 0;
-            i_bh2 = 1;
-        }
-#endif
+        printf("now finding    i_bh1=%d and i_bh2=%d\n", i_bh1, i_bh2);
+        if (++printouts >= 10) return 0;
+
+
+
+// #if defined(ACT_DEF_GRAPITE) && (defined(ADD_BH1) || defined(ADD_BH2))
+
+// #ifdef ADD_BH1
+// #define ACT_DEF_GRAPITE_NUMBH 1
+// #else
+// #define ACT_DEF_GRAPITE_NUMBH 2
+// #endif
+
+//         int act_def_grapite_bh_count = 0;
+//         for (i=0; i<n_act; i++) {
+//             if (ind_act[i]==0) {
+//                 i_bh1 = i;
+//                 act_def_grapite_bh_count++;
+//         }
+// #ifdef ADD_BH2
+//             else if (ind_act[i]==1) {
+//                 i_bh2 = i;
+//                 act_def_grapite_bh_count++;
+//             }
+// #endif
+//             if (act_def_grapite_bh_count==config->live_smbh_count) break;
+//         }
+//         if (i==n_act) {
+//             fprintf(stderr, "ERROR: black holes were not found in the active particle list");
+//             return -1;
+//         }
+// #else
+//         if (config->live_smbh_count > 0) {
+//             i_bh1 = 0;
+//             i_bh2 = 1;
+//         }
+// #endif
 
 #ifdef TIMING
         get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);