diff --git a/makefile b/makefile
index 31a29da..4619878 100644
--- a/makefile
+++ b/makefile
@@ -5,6 +5,7 @@ itpp=${pkgs}/itpp/4.3.1/lib/ -litpp
 # INC_DIR=weilei_lib
 # CXX = g++ -Wall -g -m64 -std=c++0x ${opt} ${itpp}
 CXX = g++ -Wall -g -m64 -std=c++11 -pthread ${opt} ${itpp}
+CXX = g++ -Wall -g -m64 -std=c++11 -fopenmp ${opt} ${itpp}
 CMD= -lm -lgmpxx -lgmp
 
 INC_DIR=~/working/weilei_lib
@@ -39,6 +40,8 @@ my_bp3.out:my_bp3.c $(files)
 	$(command)
 my_bp4.out:my_bp4.cpp $(files)
 	$(command)
+my_bp5.out:my_bp5.cpp $(files)
+	$(command)
 #bp_decoding3.out:bp_decoding3.c mm_read.c mm_read.h mmio.c mmio.h mm_write.c mm_write.h lib.cpp lib.h my_lib.h makefile
 bp_decoding3.out:bp_decoding3.c $(files)
 	$(command)
diff --git a/my_bp5.cpp b/my_bp5.cpp
new file mode 100644
index 0000000..e04379c
--- /dev/null
+++ b/my_bp5.cpp
@@ -0,0 +1,303 @@
+//Weilei  Apr 6, 2020
+// copied from my_bp1.c
+//#include <future>
+#include <string>
+#include <iostream>
+#include <itpp/itbase.h>
+#include <itpp/itcomm.h>
+#include <sstream>
+#include <stdio.h>
+#include "my_lib.h"
+#include <math.h>
+#include "bp_decoder.h"
+using namespace std;
+using namespace itpp;
+
+
+
+// Read the code from files and do BP decoding
+//input:source file for stabilzier matrix; error propability p ;
+
+int decode( BP_Decoder, GF2mat G, GF2mat H, double p, mat * data, int col_index, int row_index, int cycles, int feedback, double time_out, int num_data_points);
+  
+int main(int argc, char **argv){
+  Parser parser;
+  parser.init(argc,argv);
+  //p.set_silentmode(true);
+
+  string filename_G, filename_H;
+  //parser.get(filename_G,"filename_G");
+  Real_Timer timer;   double remained_time; //remained time for each size
+  //vector<future<int>> pool;
+  int cores=32; parser.get(cores, "cores");
+  //vector<future<int>>::size_type pool_size = cores+2; //max number of threads, 15 with decreasing size for best performace
+  //std::chrono::milliseconds span (100);
+  //std::chrono::milliseconds final_thread_time (10000);//10 secs before prelimilary result print
+  string filename_data;
+  filename_data="gnuplot/result/my-bp5-test.gnudat"; 
+  parser.get(filename_data,"filename_data");
+  
+  int feedback=5; parser.get(feedback,"feedback");
+  //  int cycles = 1000000;//70 sec for 2,000,000
+  //  double time_out=200;//time out in seconds for each data points (p and size)
+  int num_data_points = 100;//data entry for each data points
+  parser.get(num_data_points, "num_data_points");
+  //change parameter p, code size
+  double p;
+  //parser.get(p,"p");
+  //  int sizes[]= {13,11,9,7,5};
+  int sizes[]= {13,9,5};
+  int size_of_sizes = sizeof(sizes)/sizeof(*sizes);
+  string stabilizer_folder="data/toric/stabilizer";
+  double ip_begin=-0.7;
+  double ip_end=-1.8;    
+  int data_rows = (ip_begin-ip_end)/0.1;
+  mat data(data_rows,5 * size_of_sizes );   //return result in a mat, 5 columns for each size. format defines in header  
+  data.zeros();
+
+  int exit_iteration=9; parser.get(exit_iteration,"exit_iteration");
+  int schedule_mode=4; parser.get(schedule_mode, "schedule_mode");
+
+  //split tasks into smaller chunks
+  int chunk_num_data_points=10;//number of data points in each chunk
+  int chunk_size=num_data_points/chunk_num_data_points; //number of chunks for each p and size
+  int chunk_num_for_each_size = chunk_size*((ip_begin-ip_end)/0.1);
+
+  //timeout should be 5 times longer in hpcc
+  double chunk_time_out = 100.0;//time_out/chunk_size;
+  parser.get(chunk_time_out,"chunk_time_out");
+  int chunk_cycles=chunk_num_data_points*1000;//1000 for prob 1/1000
+
+  mat chunk_data(data_rows*chunk_size,5*5);
+  chunk_data.zeros();
+
+
+  int col_index=-5;
+  for ( int size : sizes){
+    //    cout<<"size = "<<size<<endl;
+    col_index +=5;//5 col for each size
+    filename_G = stabilizer_folder + "/toric_S_x_size_" + to_string(size) + ".mm";
+    filename_H = stabilizer_folder + "/toric_S_z_size_" + to_string(size) + ".mm";
+    //filename_result = error_folder + "/toric_S_size_" + to_string(size) + ".mm_rate";
+    GF2mat G =MM_to_GF2mat(filename_G);      
+    GF2mat H =MM_to_GF2mat(filename_H);
+    //set up for BP_Decoder
+    BP_Decoder bp_decoder;
+    bp_decoder.init(H);
+    bp_decoder.set_decode_mode_str("min sum");    
+    bp_decoder.set_exit_iteration(exit_iteration);
+    bp_decoder.set_debug_mode(false);
+    bp_decoder.set_schedule_mode(schedule_mode);
+    bp_decoder.print_info();
+      
+    int row_index=-1;
+    timer.tic();
+    for ( double ip=ip_begin;ip > ip_end;ip-=0.1){
+      //cout<<"ip = "<<ip<<endl;
+      //      row_index ++;
+      //atof(argv[4]);
+      p=pow(10,ip);
+
+      for ( int ic=0; ic<chunk_size; ic++){
+	row_index ++;
+
+	//manage thread within limit of pool size
+	/*
+  while ( pool.size() >= pool_size ){
+	  //wait until some of them finish
+	  //break;
+	  for(vector <future<int>> :: iterator it = pool.begin(); it != pool.end(); ++it){
+	    //printf("%d", *it);
+	    if ( it->wait_for(span) == future_status::ready){
+	      //	    cout<<"."<<endl;
+	      pool.erase(it);	    
+	      //cout<<"remove thread. pool_size="<<pool.size()<<endl;
+	      it--;
+	      //break;
+	    }
+	  }
+	}
+   */ 
+	// future<int> fut = async(launch::async, decode,G, H, p, filename_result_p);
+	//pool.push_back(move(fut));
+	//pool.push_back(async(launch::async, decode, bp_decoder, G, H, p, & data,col_index, row_index, cycles, feedback, time_out, num_data_points) );
+	//pool.push_back(async(launch::async, decode, bp_decoder, G, H, p, & chunk_data,col_index, row_index, chunk_cycles, feedback, chunk_time_out, chunk_num_data_points) );
+  decode( bp_decoder, G, H, p, & chunk_data,col_index, row_index, chunk_cycles, feedback, chunk_time_out, chunk_num_data_points); 
+
+	remained_time = timer.toc()/(row_index+1)*chunk_num_for_each_size;
+  std::cout
+     // cout<<"my_bp: add new thread. pool_size="<<pool.size()
+	  <<", size = "<<size
+	  <<", p = "<<p
+	  <<", row_index = "<<row_index<<"/"<<chunk_num_for_each_size
+	  <<", remained time for this p is "<< remained_time <<" sec"
+	  <<", col_index = "<<col_index
+	  <<", chunk_cycles = "<<chunk_cycles
+	  <<endl;
+      }
+    }
+  }
+  //finish adding all the threads and most of them are done
+  
+  //wait a bit for last threads
+  /*for(vector <future<int>> :: iterator it = pool.begin(); it != pool.end(); ++it){
+    it->wait_for(final_thread_time);
+  }*/
+
+  //print prelimilary result  
+  string header="# header\n# sizes: ";
+  for (int size : sizes){
+    header += to_string(size) + ", ";
+  }
+  header += ("\n# p, P_c converge rate, counts_total, counts_nonconverge, timer.toc()");
+  //  mat2gnudata(data,filename_data,header);
+  //cout<<"save prelimilary data to "<<filename_data.c_str()<<endl;
+
+  //wait all process to finish
+  /*for(vector <future<int>> :: iterator it = pool.begin(); it != pool.end(); ++it){
+    it->get();
+  }*/
+
+  //process chunk_data to data
+  int temp_index;
+  double value;
+  for ( int i =0; i<data.rows();i++){
+    for ( int j = 0; j< data.cols();j++){
+      value=0;
+      for ( int k =0;k<chunk_size; k++){
+	temp_index = i* chunk_size+k;
+	value += chunk_data(temp_index,j);
+      }
+      value = value/chunk_size;
+      data.set(i,j,value);
+    }
+  }
+
+  //print final result
+  mat2gnudata(data,filename_data,header);
+  cout<<"save final data to "<<filename_data.c_str()<<endl;
+  return 0;
+}
+
+//*********************** bp decoding function by itpp
+
+int decode( BP_Decoder bp_decoder, GF2mat G, GF2mat H, double p,  mat * data, int col_index, int row_index, int cycles, int feedback, double time_out, int num_data_points){
+  //parameter setup
+  //int decode_mode = 2;
+  //  int cycles=1000;//10000;//number of cycles: fro toric code, 10000 give reletively clear result
+  //  int exit_at_iteration=9;//parameter for bp decoding set_exit_condition()
+  //  int bound= (int) -4096 * log (p/(1-p));// -300; see note.pdf on how to choose bound
+  const int bound = 0;
+  const int max_repetition = feedback;//10 for best result. supposed to be zero in our set up, just for a check
+  
+  Real_Timer timer;
+  const int nvar = H.cols();
+  //int N = C.get_nvar(); // number of bits per codeword; N=2 x n x n is the size of the toric code.
+  
+  const bvec bitsin = zeros_b(nvar);//original zero vector
+ 
+  RNG_randomize();
+  timer.tic(); 
+
+  BSC bsc(p); // initialize BSC, error channel with probability pp
+  
+  
+  const double LLR=log( (1-p)/p);
+
+  int counts_converge=0;
+  int counts_nonconverge=0;
+  #pragma omp parallel for num_threads(8)
+  for (int i=0;i<cycles;i++){
+  if ( counts_nonconverge < num_data_points){
+    vec LLRin(nvar);
+    int ans=0; //the value return by bp_decoding(), which is the number of iterations and with a minus sign if it doesn't converge
+    //bvec rec_bits = bitsin;   // input vector with manually-input errors, for test
+
+    bvec rec_bits0 = bsc(bitsin);    // input vector with errors from bsc chanel
+    bvec rec_bits=rec_bits0;
+    QLLRvec llr_output;
+    bvec error=rec_bits;
+    bvec syndrome = (H*error);
+    //cout<<"syndrome = "<<syndrome<<endl;
+    LLRin.ones();
+    LLRin = LLRin*LLR;
+
+    vec LLRout(LLRin);//default, zero syndrome give zero error
+    ans = bp_decoder.decode(syndrome,LLRin, LLRout);
+
+    bvec bitsout = llr_output < bound; //30;//0;//output
+
+
+    //enhanced feedback
+    int current_iteration = ans;
+    for (int i_iteration = 0;i_iteration<max_repetition;i_iteration++){
+      if (ans >0){
+	break;
+      }
+      // usleep(100000);
+
+      LLRin = LLRout;
+      LLRout.zeros();
+      ans = bp_decoder.decode(syndrome,LLRin, LLRout);
+      //      ans = bp_decoder.bp_syndrome_llr(syndrome,LLRin, LLRout);
+      //      ans = bp_syndrome_llr(H,syndrome,LLRin, LLRout, exit_iteration, decode_mode);
+      bitsout = LLRout < bound;
+
+    }
+
+    bitsout=bitsout+rec_bits+rec_bits0;//cancel input vector
+    //bitsout=reduce_weight( bitsout, G);//remove trivial cycles
+    rec_bits=rec_bits0;//replace for compatibility for saving data
+
+    if(ans>=0){
+      if ( current_iteration < 0 ){
+	//cout<<"-------------------------------make it converge after iteration."<<endl;
+      }
+      counts_converge++;
+    }else{
+      counts_nonconverge++;
+      std::cout<<"counts_nonconverge = "<<counts_nonconverge<<std::endl;
+
+      /*
+      //time control
+      if (timer.toc()>time_out){
+	std::cout<<"timeout: "<<time_out<<std::endl;
+	break;
+      }
+      else if (counts_nonconverge == num_data_points/2){
+        //change number of cycles
+        cycles =min(cycles, i *2);
+        cycles=max(cycles, 50);
+        std::cout<<"update cycles to "<<cycles<<std::endl;
+        //std::cout<<i<<","<<counts_nonconverge<<", "
+        //	 << num_data_points<<std::endl;
+      }
+      */
+    }
+    }
+  }
+  //cout<<"counts_converge="<<counts_converge<<endl;
+  //double rate_converge=counts_converge*1.0/cycles;
+  //this statistical method create an overall bias
+  int counts_total = counts_converge+counts_nonconverge;
+  double rate_converge=counts_converge*1.0/counts_total;
+  cout<<"Run summary: nvar = "<<nvar<<" =  2 x "<< nvar/2;
+  cout<<", p = "<<p;
+  cout<<", Converge rate ="<<rate_converge
+      <<", cycles = "<<cycles
+      <<", counts_total = "<<counts_total
+      <<", counts_nonconverge = "<<counts_nonconverge<<endl;
+  //  save_result(p,rate_converge,filename_result_p);//no need to save this. can get it by counting the size of the matrix when doing gnuplot
+
+  //save result to mat data, using row_index and col_index
+  //data entries: p, rate
+  data->set(row_index,col_index,p);
+  data->set(row_index,col_index+1,rate_converge);
+  data->set(row_index,col_index+2,counts_total);
+  data->set(row_index,col_index+3,counts_nonconverge);
+  data->set(row_index,col_index+4,timer.toc());
+  
+  timer.toc_print();
+  return 0;
+}
+