1. Enabled OpenMP

2. Enabled AVX - earlier commit had disabled AVX

PairHMM_JNI/Makefile

@@ -1,4 +1,4 @@
-#OMPCFLAGS=-fopenmp
+OMPCFLAGS=-fopenmp
 #OMPLDFLAGS=-lgomp
 
 #CFLAGS=-O2 -std=c++11 -W -Wall -march=corei7-avx -Wa,-q            -pedantic $(OMPCFLAGS) -Wno-unknown-pragmas

PairHMM_JNI/pairhmm-1-base.cc

@@ -19,6 +19,12 @@ LoadTimeInitializer g_load_time_initializer;
 #define BATCH_SIZE  10000
 #define RUN_HYBRID
 
+double getCurrClk() {
+  struct timeval tv ;
+  gettimeofday(&tv, NULL);
+  return (double)tv.tv_sec + (double)tv.tv_usec / 1000000.0;
+}
+
 int main(int argc, char** argv)
 {
   if(argc < 2)
@@ -29,6 +35,10 @@ int main(int argc, char** argv)
   bool use_old_read_testcase = false;
   if(argc >= 3 && string(argv[2]) == "1")
     use_old_read_testcase = true;
+  unsigned chunk_size = 100;
+  if(argc >= 4)
+    chunk_size = strtol(argv[3],0,10);
+
 
   initialize_function_pointers(); 
 
@@ -45,12 +55,24 @@ int main(int argc, char** argv)
     assert(ifptr.is_open());
   }
 
+  vector<testcase> tc_vector;
+  tc_vector.clear();
   testcase tc;
   while(1)
   {
     int break_value = use_old_read_testcase ? read_testcase(&tc, fptr) : read_mod_testcase(ifptr,&tc,true);
     if(break_value < 0)
       break;
+    tc_vector.push_back(tc);
+  }
+  vector<double> results_vec;
+  results_vec.clear();
+  results_vec.resize(tc_vector.size());
+  double start_time = getCurrClk();
+#pragma omp parallel for schedule(dynamic,chunk_size)  num_threads(12)
+  for(unsigned i=0;i<tc_vector.size();++i)
+  {
+    testcase& tc = tc_vector[i];
     float result_avxf = g_compute_full_prob_float(&tc, 0);
     double result = 0;
     if (result_avxf < MIN_ACCEPTED) {
@@ -60,12 +82,23 @@ int main(int argc, char** argv)
     else
       result = (double)(log10f(result_avxf) - log10f(ldexpf(1.f, 120.f)));
 
+    results_vec[i] = result;
+  }
+  cout << "Time taken "<<getCurrClk()-start_time << "\n";
+  for(unsigned i=0;i<tc_vector.size();++i)
+  {
+    testcase& tc = tc_vector[i];
     double baseline_result = compute_full_prob<double>(&tc);
     baseline_result = log10(baseline_result) - log10(ldexp(1.0, 1020.0)); 
-    cout << std::scientific << baseline_result << " "<<result<<"\n";
-    delete tc.rs;
-    delete tc.hap;
+    double abs_error = fabs(baseline_result-results_vec[i]);
+    double rel_error = (baseline_result != 0) ? fabs(abs_error/baseline_result) : 0;
+    if(abs_error > 1e-5 && rel_error > 1e-5)
+      cout << std::scientific << baseline_result << " "<<results_vec[i]<<"\n";
+    delete tc_vector[i].rs;
+    delete tc_vector[i].hap;
   }
+  results_vec.clear();
+  tc_vector.clear();
   if(use_old_read_testcase)
     fclose(fptr);
   else

PairHMM_JNI/utils.cc

@@ -35,8 +35,7 @@ bool is_sse42_supported()
 
 void initialize_function_pointers()
 {
-  //if(is_avx_supported())
-  if(false)
+  if(is_avx_supported())
   {
     cout << "Using AVX accelerated implementation of PairHMM\n";
     g_compute_full_prob_float = compute_full_prob_avxs<float>;