From 21332566d04f5ca49f21b8a7f07a82ba7250a144 Mon Sep 17 00:00:00 2001
From: macrocosme <vohl@astron.nl>
Date: Mon, 7 Dec 2020 20:08:16 +0100
Subject: [PATCH 1/6] Functional omp version with timing output for test

---
 include/rficlean.h  |  10 +--
 src/cleanit.c       | 199 +++++++++++++++++++++++---------------------
 src/rficlean_data.c | 113 ++++++++++++++++++++-----
 3 files changed, 200 insertions(+), 122 deletions(-)

diff --git a/include/rficlean.h b/include/rficlean.h
index ef5dace..194dcf8 100644
--- a/include/rficlean.h
+++ b/include/rficlean.h
@@ -8,16 +8,16 @@ long byte_offset;
 int headerless,obits,iflip,nred, iwhite,nints,nsub,maxhist,numhist,nsubchans,tnorm,fnorm,igetstat,nharm,pcl,nsect,zerodm,bl_start,nblocks;
 char inpfile[128], outfile[128], gminfofile[128],gmhdrfile[128],psfile[100];
 FILE *input, *output, *gminfo, *gmhdr;
-double tempra,tempdec, psrf, *last_mspec, *wrms, *wmean, *ai;
+double tempra,tempdec, psrf, *last_mspec, **wrms, **wmean, **ai;
 float forcefthresh,fthresh,rthresh,sthresh,last_tvar,last_fvar,chanfrac,sampfrac,clipthresh;
 float *fftstat, *chanstat, *predist, *xpredist, *postdist, *xpostdist, *finaldist, *xfinaldist ;
 double *tfvar, *tfmean,meanvar,rmsvar ;
-double *chandata, *mspec, *rspec, *vspec, *wspec, *wt;
+double **chandata, *mspec, *rspec, *vspec, *wspec, *wt;
 long int *coff;
 
-fftw_complex *in;
-fftw_complex *out;
-fftw_plan fplan, bplan;
+fftw_complex **in;
+fftw_complex **out;
+fftw_plan *fplan, *bplan;
 
 
 int strings_equal (char *string1, char *string2);
diff --git a/src/cleanit.c b/src/cleanit.c
index c898aa7..4a302f4 100644
--- a/src/cleanit.c
+++ b/src/cleanit.c
@@ -1,16 +1,17 @@
 /*
  * cleanit  --- the actual cleaning/processing of data.
- * 
- * 
+ *
+ *
  * Yogesh Maan <maan@astron.nl>   2018.
- * 
- * 
+ *
+ *
  */
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #include "rficlean.h"
 #include <fftw3.h>
+#include <omp.h>
 
 //--------------------------------------------------------------
 int all_samef(double *arr, long int np)
@@ -133,7 +134,7 @@ void fftclean(fftw_complex *in, fftw_complex *out, long int npts, int ioff)
 
   nfft = npts ;
   nfft2 = nfft/2;
-  fftw_execute ( fplan );
+  fftw_execute ( fplan[omp_get_thread_num()] );
   anfft = (double)nfft;
   for (i=0;i<nfft;i++) {
       in[i][0] = out[i][0];
@@ -145,40 +146,40 @@ void fftclean(fftw_complex *in, fftw_complex *out, long int npts, int ioff)
     for (ii=0;ii<nred;ii++){
       kk = redplan[ii];
       for (j=0;j<kk;j++){
-        ai[j] = out[k][1];  // imaginary part
+        ai[omp_get_thread_num()][j] = out[k][1];  // imaginary part
         k++;
       }
-      robust_meanrms(ai,kk);
-      wmean[ii] = ai[kk];
-      wrms[ii] = ai[kk+1];
+      robust_meanrms(ai[omp_get_thread_num()],kk);
+      wmean[omp_get_thread_num()][ii] = ai[omp_get_thread_num()][kk];
+      wrms[omp_get_thread_num()][ii] = ai[omp_get_thread_num()][kk+1];
     }
 
   }    // end of block for whitened estimates of mean,rms
   else {
     k = nfft-1;
-    for (i=0;i<k;i++) ai[i]=out[i+1][1];  // exclude DC
-    ai[k] = ai[k+1] = 0.0;
-    robust_meanrms(ai,k);
+    for (i=0;i<k;i++) ai[omp_get_thread_num()][i]=out[i+1][1];  // exclude DC
+    ai[omp_get_thread_num()][k] = ai[omp_get_thread_num()][k+1] = 0.0;
+    robust_meanrms(ai[omp_get_thread_num()],k);
     for (ii=0;ii<nred;ii++){
-      wmean[ii] = ai[k];
-      wrms[ii] = ai[k+1];
+      wmean[omp_get_thread_num()][ii] = ai[omp_get_thread_num()][k];
+      wrms[omp_get_thread_num()][ii] = ai[omp_get_thread_num()][k+1];
     }
   }  // end of no-whitening block
 
   // now the actual cleaning
   i1 = 1;  // exclude dc
   for (ii=0;ii<nred;ii++){
-    ath = fthresh*wrms[ii] + wmean[ii];
+    ath = fthresh*wrms[omp_get_thread_num()][ii] + wmean[omp_get_thread_num()][ii];
     bth = ath;
     i2 = i1 + redplan[ii];
     for (i=i1;i<i2;i++) {
       if( (fabs(out[i][0])>ath || fabs(out[i][1])>bth) && psrifs[i]>0 ){
-        fftstat[ioff+i] = fftstat[ioff+i] + 1.0;
+        fftstat[ioff+i] += 1.0;
         for (j=i,ih=1; ih<=nharm && j<nfft2; ih++,j=i*ih){
-          in[j][0] = wmean[ii];
-          in[j][1] = wmean[ii];
-          in[nfft-j][0] = wmean[ii];
-          in[nfft-j][1] = wmean[ii];
+          in[j][0] = wmean[omp_get_thread_num()][ii];
+          in[j][1] = wmean[omp_get_thread_num()][ii];
+          in[nfft-j][0] = wmean[omp_get_thread_num()][ii];
+          in[nfft-j][1] = wmean[omp_get_thread_num()][ii];
         }
      }
     }
@@ -197,13 +198,13 @@ void fftclean(fftw_complex *in, fftw_complex *out, long int npts, int ioff)
       in[nfft-j-1][0] = out[nfft-j-1][0];
       in[nfft-j-1][1] = out[nfft-j-1][1];
       in[j-1][0] = out[j-1][0];
-      in[j-1][1] = out[j-1][1]; 
+      in[j-1][1] = out[j-1][1];
       in[nfft-j+1][0] = out[nfft-j+1][0];
-      in[nfft-j+1][1] = out[nfft-j+1][1]; 
+      in[nfft-j+1][1] = out[nfft-j+1][1];
     }
   }  */
 
-  fftw_execute ( bplan );
+  fftw_execute ( bplan[omp_get_thread_num()] );
   for (i=0;i<nfft;i++) {
     in[i][0] = out[i][0]/anfft;
   }
@@ -219,19 +220,19 @@ void tsclip(double *tdata, long int npts, float thresh)
   double anpt,ath;
   double m1,r1;
 
-  for (i=0;i<npts;i++) ai[i]=tdata[i];
-  robust_meanrms(ai,npts);
-  m1 = ai[npts];
-  r1 = ai[npts+1];
+  for (i=0;i<npts;i++) ai[omp_get_thread_num()][i]=tdata[i];
+  robust_meanrms(ai[omp_get_thread_num()],npts);
+  m1 = ai[omp_get_thread_num()][npts];
+  r1 = ai[omp_get_thread_num()][npts+1];
   if(r1<=0.0){
     r1 = 0.0;
     m1 = 0.0;
   }
   ath = m1+thresh*r1;
   for (i=0;i<npts;i++) {
-    if( fabs(ai[i])>ath) ai[i] = (ai[i]/fabs(ai[i]))*ath;
+    if( fabs(ai[omp_get_thread_num()][i])>ath) ai[omp_get_thread_num()][i] = (ai[omp_get_thread_num()][i]/fabs(ai[omp_get_thread_num()][i]))*ath;
   }
-  for (i=0;i<npts;i++) tdata[i]=ai[i];
+  for (i=0;i<npts;i++) tdata[i]=ai[omp_get_thread_num()][i];
   tdata[npts] = m1;
   tdata[npts+1] = r1;
 }
@@ -245,19 +246,19 @@ void tsclean(double *tdata, long int npts, float thresh)
   double anpt,ath;
   double m1,r1;
 
-  for (i=0;i<npts;i++) ai[i]=tdata[i];
-  robust_meanrms(ai,npts);
-  m1 = ai[npts];
-  r1 = ai[npts+1];
+  for (i=0;i<npts;i++) ai[omp_get_thread_num()][i]=tdata[i];
+  robust_meanrms(ai[omp_get_thread_num()],npts);
+  m1 = ai[omp_get_thread_num()][npts];
+  r1 = ai[omp_get_thread_num()][npts+1];
   if(r1<=0.0){
     r1 = 0.0;
     m1 = 0.0;
   }
   ath = thresh*r1;
   for (i=0;i<npts;i++) {
-    if( fabs((ai[i]-m1))>ath) ai[i] = m1;
+    if( fabs((ai[omp_get_thread_num()][i]-m1))>ath) ai[omp_get_thread_num()][i] = m1;
   }
-  for (i=0;i<npts;i++) tdata[i]=ai[i];
+  for (i=0;i<npts;i++) tdata[i]=ai[omp_get_thread_num()][i];
   tdata[npts] = m1;
   tdata[npts+1] = r1;
 }
@@ -271,17 +272,17 @@ void spfind(double *cdata, long int npts, float thresh, double *wt)
   double m1,r1,m2;
 
   n = npts-1;
-  for (i=0;i<n;i++) ai[i]=fabs(cdata[i+1]-cdata[i]) ;
-  ai[n] = 0.0;
-  ai[n+1] = 0.0;
-  robust_meanrms(ai,n);
-  m1 = ai[n];
-  r1 = ai[n+1];
+  for (i=0;i<n;i++) ai[omp_get_thread_num()][i]=fabs(cdata[i+1]-cdata[i]) ;
+  ai[omp_get_thread_num()][n] = 0.0;
+  ai[omp_get_thread_num()][n+1] = 0.0;
+  robust_meanrms(ai[omp_get_thread_num()],n);
+  m1 = ai[omp_get_thread_num()][n];
+  r1 = ai[omp_get_thread_num()][n+1];
   if(r1<=0.0) r1 = 10000.0;
 
   ath = thresh*r1;
   for (i=0;i<n;i++) {
-    if( fabs(ai[i]-m1)>ath){
+    if( fabs(ai[omp_get_thread_num()][i]-m1)>ath){
       wt[i] = -1.0;
       wt[i+1] = -1.0;
     }
@@ -298,17 +299,17 @@ void tsfind(double *ttdata, long int npts, float thresh, double *wt)
   double m1,r1,m2;
 
 
-  for (i=0;i<npts;i++) ai[i]=ttdata[i] ;
-  ai[npts] = 0.0;
-  ai[npts+1] = 0.0;
-  robust_meanrms(ai,npts);
-  m1 = ai[npts];
-  r1 = ai[npts+1];
+  for (i=0;i<npts;i++) ai[omp_get_thread_num()][i]=ttdata[i] ;
+  ai[omp_get_thread_num()][npts] = 0.0;
+  ai[omp_get_thread_num()][npts+1] = 0.0;
+  robust_meanrms(ai[omp_get_thread_num()],npts);
+  m1 = ai[omp_get_thread_num()][npts];
+  r1 = ai[omp_get_thread_num()][npts+1];
   if(r1<=0.0) r1 = 10000.0;
 
   ath = thresh*r1;
   for (i=0;i<npts;i++) {
-    if( fabs((ai[i]-m1))>ath){
+    if( fabs((ai[omp_get_thread_num()][i]-m1))>ath){
       wt[i] = -1.0;
     }
   }
@@ -358,31 +359,35 @@ void cleanit(float *data, int nchans, long int nadd)
 // get some pre-cleaning statistics
 
   for (channum=0; channum<nchans; channum++) coff[channum]=nadd*channum ;
+
+  #pragma omp parallel for private(channum,inc,isame,ii,jj) schedule(dynamic)
   for (channum=0; channum<nchans; channum++) {
     inc = coff[channum];
   /* Select the correct channel */
-    for (ii = 0, jj = channum; ii < nadd; ii++, jj += nchans) chandata[ii] = (double) data[jj];
-    for (ii = 0; ii<nadd; ii++)
-    { in[ii][0] = chandata[ii];
-      in[ii][1] = 0.0; }
-    if( pcl == 0 ){
+    for (ii = 0, jj = channum; ii < nadd; ii++, jj += nchans) chandata[omp_get_thread_num()][ii] = (double) data[jj];
+    for (ii = 0; ii<nadd; ii++) {
+      in[omp_get_thread_num()][ii][0] = chandata[omp_get_thread_num()][ii];
+      in[omp_get_thread_num()][ii][1] = 0.0;
+    }
+    if( pcl == 0 ) {
       isame = 0;
     } else {
-      isame = all_samef(chandata,nadd);
+      isame = all_samef(chandata[omp_get_thread_num()],nadd);
     }
-    if( isame == 0){
-      fftclean(in,out,nadd,inc); // clean some periodic RFIs
-      for (ii = 0; ii<nadd; ii++) chandata[ii] = in[ii][0];
-      tsclip(chandata,nadd,sthresh); // clean some spiky RFIs
-      mspec[channum]=chandata[nadd];
-      rspec[channum]=chandata[nadd+1];
+    if( isame == 0) {
+      fftclean(in[omp_get_thread_num()],out[omp_get_thread_num()],nadd,inc); // clean some periodic RFIs
+      for (ii = 0; ii<nadd; ii++) chandata[omp_get_thread_num()][ii] = in[omp_get_thread_num()][ii][0];
+      tsclip(chandata[omp_get_thread_num()],nadd,sthresh); // clean some spiky RFIs
+      mspec[channum]=chandata[omp_get_thread_num()][nadd];
+      rspec[channum]=chandata[omp_get_thread_num()][nadd+1];
     }
     else {
-      mspec[channum]=chandata[0];
+      mspec[channum]=chandata[omp_get_thread_num()][0];
       rspec[channum]=0.0;
     }
     vspec[channum] = rspec[channum]*rspec[channum];
-    for (ii = 0, jj = channum; ii < nadd; ii++, jj += nchans) data[jj]= (float) chandata[ii];
+    for (ii = 0, jj = channum; ii < nadd; ii++, jj += nchans)
+      data[jj] = (float) chandata[omp_get_thread_num()][ii];
   }
 // get some post- fft and timeseries cleaning statistics
 
@@ -392,19 +397,19 @@ void cleanit(float *data, int nchans, long int nadd)
   spfind(vspec,nchans,rthresh,wspec);
   tsfind(mspec,nchans,rthresh,wspec);
   tsfind(vspec,nchans,rthresh,wspec);
-  if(pcl>0){
+  if(pcl>0) {
      kk = (int)(0.9*nchans);
      for (i=kk;i<nchans;i++) wspec[i]=0.0;
      for (i=0;i<kk;i++) wspec[i]=0.0;
   }
-  for (i=0;i<nchans; i++){
+  for (i=0;i<nchans; i++) {
     if(wspec[i] > 0.0) {
       lg = mspec[i];
       lgr = rspec[i];
       break;
     }
   }
-  for (i=0;i<nchans;i++){
+  for (i=0;i<nchans;i++) {
     if(wspec[i]<0.0) {
       mspec[i] = lg;
       rspec[i] = lgr;
@@ -417,8 +422,8 @@ void cleanit(float *data, int nchans, long int nadd)
 
   // log the information which channels got rejected
   for (i=0;i<nchans;i++) if(wspec[i]<0.0) chanstat[i]=chanstat[i]+1.0;
-  
 
+  #pragma omp parallel for private(channum,jj,ii)
   for (channum=0; channum<nchans; channum++) {
     if(wspec[channum] < 0.0) {
        for (ii = 0, jj = channum; ii < nadd; ii++, jj += nchans) data[jj] = mspec[channum];
@@ -426,34 +431,40 @@ void cleanit(float *data, int nchans, long int nadd)
   }
 // see if the whole block needs to be rejected
   an = 0.0;
-  for (ii=0; ii<nchans; ii++){ if(wspec[ii]<0.0) an=an+1.0; }
-  if((an/nchans)>chanfrac){
+  for (ii=0; ii<nchans; ii++) if(wspec[ii]<0.0) an=an+1.0;
+  if((an/nchans)>chanfrac) {
+    #pragma omp parallel for private(channum,ii,jj) collapse(2)
     for (channum=0; channum<nchans; channum++) {
-         for (ii = 0, jj = channum; ii < nadd; ii++, jj += nchans) data[jj] = last_mspec[channum]; }
+      for (ii = 0; ii < nadd; ii++) {
+          jj = channum + (ii * nchans);
+          data[jj] = last_mspec[channum];
+      }
+    }
   } else {
-    for (ii=0; ii<nchans; ii++){ if(wspec[channum] > 0.0) last_mspec[ii] = mspec[ii];}
+    for (ii=0; ii<nchans; ii++) if(wspec[channum] > 0.0) last_mspec[ii] = mspec[ii];
   }
-  
+
 
 
 // Try clipping some channels in individual samples
+  #pragma omp parallel for private(t,c,nxc) schedule(dynamic)
   for (t=0; t<nadd; t++){
       nxc = t*nchans;
-      for (c=0; c<nchans; c++) chandata[c] = data[nxc+c];
-      spclip(chandata,nchans,clipthresh);
-      for (c=0; c<nchans; c++) data[nxc+c] = chandata[c];
+      for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][c] = data[nxc+c];
+      spclip(chandata[omp_get_thread_num()],nchans,clipthresh);
+      for (c=0; c<nchans; c++) data[nxc+c] = chandata[omp_get_thread_num()][c];
   }
 // Now some timeseries cleaning
   an = (double)nadd;
   for (t=0; t<nadd; t++){
       nxc = t*nchans;
-      chandata[t]=0.0;
-      for (c=0; c<nchans; c++) chandata[t] = chandata[t]+data[nxc+c];
-      chandata[t]=chandata[t]/an;
+      chandata[omp_get_thread_num()][t]=0.0;
+      for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][t] = chandata[omp_get_thread_num()][t]+data[nxc+c];
+      chandata[omp_get_thread_num()][t]=chandata[omp_get_thread_num()][t]/an;
   }
   for (i=0; i<nadd; i++) wt[i]=+1.0;
-  tsfind(chandata,nadd,sthresh,wt);
-  spfind(chandata,nadd,sthresh,wt);
+  tsfind(chandata[omp_get_thread_num()],nadd,sthresh,wt);
+  spfind(chandata[omp_get_thread_num()],nadd,sthresh,wt);
   for (t=0; t<nadd; t++){
       nxc = t*nchans;
       if(wt[t] < 0.0){
@@ -470,8 +481,8 @@ void cleanit(float *data, int nchans, long int nadd)
       isame = 0;
       for(jj=0; jj<nvar; jj++){
         ii = nxc + jj*kk;
-        for (c=0; c<kk; c++) chandata[c] = data[ii+c];
-        if(all_samef(chandata,kk)>0) isame = isame + 1;
+        for (c=0; c<kk; c++) chandata[omp_get_thread_num()][c] = data[ii+c];
+        if(all_samef(chandata[omp_get_thread_num()],kk)>0) isame = isame + 1;
       }
       if(isame >= (int)(0.8*nvar)){
         for (c=0; c<nchans; c++) data[c+nxc] = 0.0;
@@ -489,10 +500,10 @@ void cleanit(float *data, int nchans, long int nadd)
            rspec[t] = 0.0;
           }
           else{
-           for (c=0; c<nchans; c++) chandata[c] = data[nxc+c];
-           robust_meanrms(chandata,nchans);
-           //simple_meanrms(chandata,nchans);
-           rspec[t] = chandata[nchans+1]*chandata[nchans+1];
+           for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][c] = data[nxc+c];
+           robust_meanrms(chandata[omp_get_thread_num()],nchans);
+           //simple_meanrms(chandata[omp_get_thread_num()],nchans);
+           rspec[t] = chandata[omp_get_thread_num()][nchans+1]*chandata[omp_get_thread_num()][nchans+1];
           }
       }
       t = nadd;
@@ -502,21 +513,21 @@ void cleanit(float *data, int nchans, long int nadd)
     }
     nvar = nadd/nsect;
     for (jj=0; jj<nsect; jj++){
-      for (c=0; c<nchans; c++) chandata[c] = 0.0;
+      for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][c] = 0.0;
       for (t=0; t<nvar; t++){
         nxc = (jj*nvar+t)*nchans;
-        for (c=0; c<nchans; c++) chandata[c] = chandata[c]+data[nxc+c];
+        for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][c] = chandata[omp_get_thread_num()][c]+data[nxc+c];
       }
-      for (c=0; c<nchans; c++) chandata[c] = chandata[c]/nvar;
-      //robust_meanrms(chandata,nchans);
-      simple_meanrms(chandata,nchans);
-      an = chandata[nchans+1]*chandata[nchans+1];
+      for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][c] = chandata[omp_get_thread_num()][c]/nvar;
+      //robust_meanrms(chandata[omp_get_thread_num()],nchans);
+      simple_meanrms(chandata[omp_get_thread_num()],nchans);
+      an = chandata[omp_get_thread_num()][nchans+1]*chandata[omp_get_thread_num()][nchans+1];
       if(an<=0.0) an=10000000.0;
       for (t=0; t<nvar; t++){
         nxc = (jj*nvar+t)*nchans;
         for (c=0; c<nchans; c++) data[nxc+c] = data[nxc+c]*last_tvar/an;
       }
-    } 
+    }
   }
 
   // Try weighting channels by variance
@@ -532,7 +543,7 @@ void cleanit(float *data, int nchans, long int nadd)
       else {
          for (ii = 0, jj = channum; ii < nadd; ii++, jj += nchans) data[jj] = 0.0;
       }
-    }   
+    }
   }
 
 
diff --git a/src/rficlean_data.c b/src/rficlean_data.c
index 4c6aac5..23950e8 100644
--- a/src/rficlean_data.c
+++ b/src/rficlean_data.c
@@ -1,8 +1,8 @@
 /*
- * 
- * rficlean_data --- handle reading/writing of data blocks, and 
+ *
+ * rficlean_data --- handle reading/writing of data blocks, and
  *                   intitialization etc.
- * 
+ *
  * Yogesh Maan <maan@astron.nl>   2018
  *
  */
@@ -12,9 +12,11 @@
 #include <math.h>
 #include "rficlean.h"
 #include <fftw3.h>
+#include <omp.h>
+#include <time.h>
 
-void rficlean_data(FILE *input, FILE *output) 
-{ 
+void rficlean_data(FILE *input, FILE *output)
+{
   char string[80],plotdevice[100];
   float *fblock,min,max,hpower,realtime,fsaved[2], *ts0dm, atemp;
   unsigned short *sblock;
@@ -22,6 +24,13 @@ void rficlean_data(FILE *input, FILE *output)
   int nsaved=0,opened=0, wpout=10;
   long int ns,nsblk,nout,iter,i,j,k, jt1, jt2, iblock;
   long int itemp, isum, nsize, istart, ii, jj, kk, n0;
+  time_t s, start, cpu_time;
+  clock_t t_global, t_function, t_function_sum;
+  int i_function;
+  i_function = 0;
+
+  s = time(NULL);
+  t_global = clock();
 
   last_tvar = -1.0;
   last_fvar = -1.0;
@@ -41,7 +50,7 @@ void rficlean_data(FILE *input, FILE *output)
 
   printf ("\n Preparing all buffers...\n");
 /***/
-    chandata = (double *) malloc(nsize*sizeof(double));
+    // chandata = (double *) malloc(nsize*sizeof(double));
     vspec = (double *) malloc(nsize*sizeof(double));
     rspec = (double *) malloc(nsize*sizeof(double));
   coff = (long int *) malloc((nchans)*sizeof(long int));
@@ -50,17 +59,34 @@ void rficlean_data(FILE *input, FILE *output)
   wspec = (double *) malloc(nsize*sizeof(double));
 /***/
 
-  in = fftw_malloc (sizeof(fftw_complex)*(naddt+5));
-  out = fftw_malloc (sizeof(fftw_complex)*(naddt+5));
-  fplan = fftw_plan_dft_1d( naddt, in, out,FFTW_FORWARD, FFTW_ESTIMATE );
-  bplan = fftw_plan_dft_1d( naddt, in, out,FFTW_BACKWARD, FFTW_ESTIMATE );
-  ai = (double *) malloc(sizeof(double)*(nsize));
+  // in = fftw_malloc (sizeof(fftw_complex)*(naddt+5));
+  // out = fftw_malloc (sizeof(fftw_complex)*(naddt+5));
+  // fplan = fftw_plan_dft_1d( naddt, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
+  // bplan = fftw_plan_dft_1d( naddt, in, out, FFTW_BACKWARD, FFTW_ESTIMATE );
+  // ai = (double *) malloc(sizeof(double)*(nsize));
   psrifs = (long int *) malloc(sizeof(long int)*(naddt));
 
+  #ifdef _OPENMP
+  if (fftw_init_threads()) fftw_plan_with_nthreads(omp_get_max_threads());
+  #endif
 
+  chandata = (double **) malloc(omp_get_max_threads() * sizeof(double *));
+  in = (fftw_complex**) malloc (omp_get_max_threads() * sizeof(fftw_complex*));
+  out = (fftw_complex**) malloc (omp_get_max_threads() * sizeof(fftw_complex*));
+  fplan = (fftw_plan*) malloc (omp_get_max_threads() * sizeof(fftw_plan));
+  bplan = (fftw_plan*) malloc (omp_get_max_threads() * sizeof(fftw_plan));
+  ai = (double **) malloc(omp_get_max_threads() * sizeof(double *));
+  for (i=0; i<omp_get_max_threads(); i++) {
+    chandata[i] = (double *) malloc(nsize * sizeof(double));
+    ai[i] = (double *) malloc(nsize * sizeof(double));
+    in[i] = fftw_malloc ((naddt+5) * sizeof(fftw_complex));
+    out[i] = fftw_malloc ((naddt+5) * sizeof(fftw_complex));
+    fplan[i] = fftw_plan_dft_1d(naddt, in[i], out[i], FFTW_FORWARD, FFTW_ESTIMATE);
+    bplan[i] = fftw_plan_dft_1d(naddt, in[i], out[i], FFTW_BACKWARD, FFTW_ESTIMATE);
+  }
 
   printf (" Preparing the median-filtering plan ... ");
-  // plan the median-filter window sizes for spectral whitening 
+  // plan the median-filter window sizes for spectral whitening
   //itemp = 16; // starting length of median-filter
   itemp = 8; // starting length of median-filter
   iter = 0;
@@ -76,8 +102,14 @@ void rficlean_data(FILE *input, FILE *output)
   iter++;
   redplan[iter] = naddt/2-isum-1;
   nred = iter+1;
-  wmean = (double *) malloc((nred)*sizeof(double));
-  wrms = (double *) malloc((nred)*sizeof(double));
+  // wmean = (double *) malloc((nred)*sizeof(double));
+  // wrms = (double *) malloc((nred)*sizeof(double));
+  wmean = (double **) malloc(omp_get_max_threads() * sizeof(double *));
+  wrms = (double **) malloc(omp_get_max_threads() * sizeof(double *));
+  for (i=0; i<omp_get_max_threads(); i++) {
+    wmean[i] = (double *) malloc(nred * sizeof(double));
+    wrms[i] = (double *) malloc(nred * sizeof(double));
+  }
   printf (" Done! \n");
 
   // identify the fft-bin corresponding to pulsar's fundamental to be protected
@@ -158,7 +190,12 @@ void rficlean_data(FILE *input, FILE *output)
       jj = naddt*nchans;
       for (j=ns;j<jj;j++) fblock[j] = 0.0;
     }
+
+    if (i_function == 0) t_function_sum = t_function = clock();
+    else t_function = clock();
     cleanit(fblock,nchans,naddt);
+    t_function_sum += ((double) clock() - t_function);
+    i_function += 1;
     //-------------------------------------------------
     // compute post-cleaning 0-DM tseries
     for (j=0;j<n0;j++) {
@@ -183,7 +220,7 @@ void rficlean_data(FILE *input, FILE *output)
     //-------------------------------------------------
 
     istart = istart+n0;
-    
+
 
     if (!opened) {
       opened=1;
@@ -250,13 +287,13 @@ void rficlean_data(FILE *input, FILE *output)
   fclose(input);
   fclose(output);
 
-  fftw_destroy_plan ( fplan );
-  fftw_destroy_plan ( bplan );
-  fftw_free ( in );
-  fftw_free ( out );
-  free (ai);
-  free (wmean);
-  free (wrms);
+  // fftw_destroy_plan ( fplan );
+  // fftw_destroy_plan ( bplan );
+  // fftw_free ( in );
+  // fftw_free ( out );
+  // free (ai);
+  // free (wmean);
+  // free (wrms);
   free (tfvar);
   free (tfmean);
   free (fblock);
@@ -279,7 +316,37 @@ void rficlean_data(FILE *input, FILE *output)
   free(wspec);
   free(wt);
   free(coff);
-  free(chandata);
+
+  for (i=0; i<omp_get_max_threads(); i++) {
+    fftw_free (in[i]);
+    fftw_free (out[i]);
+    fftw_destroy_plan ( fplan[i] );
+    fftw_destroy_plan ( bplan[i] );
+    free (ai[i]);
+    free (wmean[i]);
+    free (wrms[i]);
+    free(chandata[i]);
+  }
+  free (fplan);
+  free (bplan);
+  free (in);
+  free (out);
+  free (ai);
+  free (wmean);
+  free (wrms);
+  free (chandata);
+
+  FILE *file_timing = fopen("timing.csv", "a+");
+  fprintf(file_timing, "%d,%d,%d,%.3f,%.3f\n",
+    naddt,
+    omp_get_max_threads(),
+    time(NULL) - s,
+    ((double) clock() - t_global) / (CLOCKS_PER_SEC) / (double) omp_get_max_threads(),
+    ((double) t_function_sum / (i_function * CLOCKS_PER_SEC)) / (double) omp_get_max_threads()
+  );
+  // fprintf(file_timing, "%s,%d,%d\n", time_section, -1, time(NULL) - s);
+  fclose(file_timing);
 
   printf (" Freed the buffers. All done!\n\n ");
+  printf ("%d seconds\n\n", time(NULL) - s);
 }

From ad61447c20ecb71a86251ce915ae03e78ba99ff6 Mon Sep 17 00:00:00 2001
From: Dany Vohl <macrocosme@users.noreply.github.com>
Date: Mon, 7 Dec 2020 20:14:25 +0100
Subject: [PATCH 2/6] Add threads and omp

---
 Makefile | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Makefile b/Makefile
index b11637e..e619edb 100644
--- a/Makefile
+++ b/Makefile
@@ -6,8 +6,8 @@ BINDIR=bin
 MYBIN=/home/maan/pulsar_softwares/bin/
 
 CC=gcc
-CFLAGS=-I$(IDIR) -Wno-unused-result -O3 -march=native
-LIBS=-lm -lfftw3 -lcpgplot
+CFLAGS=-I$(IDIR) -Wno-unused-result -O3 -march=native -fopenmp
+LIBS=-lm -lfftw3_threads -lfftw3 -lcpgplot
 
 _DEPS = header.h  rficlean.h  version.h
 DEPS = $(patsubst %,$(IDIR)/%,$(_DEPS))

From ee6dcda7645df48d8f49fc0575b3580215d53150 Mon Sep 17 00:00:00 2001
From: macrocosme <vohl@astron.nl>
Date: Mon, 11 Jan 2021 15:34:04 +0100
Subject: [PATCH 3/6] Adding ifdefs _OPENMP

---
 src/cleanit.c       | 186 ++++++++++++++++++++++++++------------------
 src/rficlean_data.c |  34 ++++----
 2 files changed, 131 insertions(+), 89 deletions(-)

diff --git a/src/cleanit.c b/src/cleanit.c
index 4a302f4..839efef 100644
--- a/src/cleanit.c
+++ b/src/cleanit.c
@@ -134,7 +134,15 @@ void fftclean(fftw_complex *in, fftw_complex *out, long int npts, int ioff)
 
   nfft = npts ;
   nfft2 = nfft/2;
-  fftw_execute ( fplan[omp_get_thread_num()] );
+
+  #ifdef _OPENMP
+  int thread_num = omp_get_thread_num();
+  #else
+  int thread_num=0;
+  #endif
+
+  fftw_execute ( fplan[thread_num] );
+
   anfft = (double)nfft;
   for (i=0;i<nfft;i++) {
       in[i][0] = out[i][0];
@@ -146,40 +154,40 @@ void fftclean(fftw_complex *in, fftw_complex *out, long int npts, int ioff)
     for (ii=0;ii<nred;ii++){
       kk = redplan[ii];
       for (j=0;j<kk;j++){
-        ai[omp_get_thread_num()][j] = out[k][1];  // imaginary part
+        ai[thread_num][j] = out[k][1];  // imaginary part
         k++;
       }
-      robust_meanrms(ai[omp_get_thread_num()],kk);
-      wmean[omp_get_thread_num()][ii] = ai[omp_get_thread_num()][kk];
-      wrms[omp_get_thread_num()][ii] = ai[omp_get_thread_num()][kk+1];
+      robust_meanrms(ai[thread_num],kk);
+      wmean[thread_num][ii] = ai[thread_num][kk];
+      wrms[thread_num][ii] = ai[thread_num][kk+1];
     }
 
   }    // end of block for whitened estimates of mean,rms
   else {
     k = nfft-1;
-    for (i=0;i<k;i++) ai[omp_get_thread_num()][i]=out[i+1][1];  // exclude DC
-    ai[omp_get_thread_num()][k] = ai[omp_get_thread_num()][k+1] = 0.0;
-    robust_meanrms(ai[omp_get_thread_num()],k);
+    for (i=0;i<k;i++) ai[thread_num][i]=out[i+1][1];  // exclude DC
+    ai[thread_num][k] = ai[thread_num][k+1] = 0.0;
+    robust_meanrms(ai[thread_num],k);
     for (ii=0;ii<nred;ii++){
-      wmean[omp_get_thread_num()][ii] = ai[omp_get_thread_num()][k];
-      wrms[omp_get_thread_num()][ii] = ai[omp_get_thread_num()][k+1];
+      wmean[thread_num][ii] = ai[thread_num][k];
+      wrms[thread_num][ii] = ai[thread_num][k+1];
     }
   }  // end of no-whitening block
 
   // now the actual cleaning
   i1 = 1;  // exclude dc
   for (ii=0;ii<nred;ii++){
-    ath = fthresh*wrms[omp_get_thread_num()][ii] + wmean[omp_get_thread_num()][ii];
+    ath = fthresh*wrms[thread_num][ii] + wmean[thread_num][ii];
     bth = ath;
     i2 = i1 + redplan[ii];
     for (i=i1;i<i2;i++) {
       if( (fabs(out[i][0])>ath || fabs(out[i][1])>bth) && psrifs[i]>0 ){
         fftstat[ioff+i] += 1.0;
         for (j=i,ih=1; ih<=nharm && j<nfft2; ih++,j=i*ih){
-          in[j][0] = wmean[omp_get_thread_num()][ii];
-          in[j][1] = wmean[omp_get_thread_num()][ii];
-          in[nfft-j][0] = wmean[omp_get_thread_num()][ii];
-          in[nfft-j][1] = wmean[omp_get_thread_num()][ii];
+          in[j][0] = wmean[thread_num][ii];
+          in[j][1] = wmean[thread_num][ii];
+          in[nfft-j][0] = wmean[thread_num][ii];
+          in[nfft-j][1] = wmean[thread_num][ii];
         }
      }
     }
@@ -204,7 +212,7 @@ void fftclean(fftw_complex *in, fftw_complex *out, long int npts, int ioff)
     }
   }  */
 
-  fftw_execute ( bplan[omp_get_thread_num()] );
+  fftw_execute ( bplan[thread_num] );
   for (i=0;i<nfft;i++) {
     in[i][0] = out[i][0]/anfft;
   }
@@ -220,19 +228,25 @@ void tsclip(double *tdata, long int npts, float thresh)
   double anpt,ath;
   double m1,r1;
 
-  for (i=0;i<npts;i++) ai[omp_get_thread_num()][i]=tdata[i];
-  robust_meanrms(ai[omp_get_thread_num()],npts);
-  m1 = ai[omp_get_thread_num()][npts];
-  r1 = ai[omp_get_thread_num()][npts+1];
+  #ifdef _OPENMP
+  int thread_num = omp_get_thread_num();
+  #else
+  int thread_num = 0;
+  #endif
+
+  for (i=0;i<npts;i++) ai[thread_num][i]=tdata[i];
+  robust_meanrms(ai[thread_num],npts);
+  m1 = ai[thread_num][npts];
+  r1 = ai[thread_num][npts+1];
   if(r1<=0.0){
     r1 = 0.0;
     m1 = 0.0;
   }
   ath = m1+thresh*r1;
   for (i=0;i<npts;i++) {
-    if( fabs(ai[omp_get_thread_num()][i])>ath) ai[omp_get_thread_num()][i] = (ai[omp_get_thread_num()][i]/fabs(ai[omp_get_thread_num()][i]))*ath;
+    if( fabs(ai[thread_num][i])>ath) ai[thread_num][i] = (ai[thread_num][i]/fabs(ai[thread_num][i]))*ath;
   }
-  for (i=0;i<npts;i++) tdata[i]=ai[omp_get_thread_num()][i];
+  for (i=0;i<npts;i++) tdata[i]=ai[thread_num][i];
   tdata[npts] = m1;
   tdata[npts+1] = r1;
 }
@@ -246,19 +260,25 @@ void tsclean(double *tdata, long int npts, float thresh)
   double anpt,ath;
   double m1,r1;
 
-  for (i=0;i<npts;i++) ai[omp_get_thread_num()][i]=tdata[i];
-  robust_meanrms(ai[omp_get_thread_num()],npts);
-  m1 = ai[omp_get_thread_num()][npts];
-  r1 = ai[omp_get_thread_num()][npts+1];
+  #ifdef _OPENMP
+  int thread_num = omp_get_thread_num();
+  #else
+  int thread_num=0;
+  #endif
+
+  for (i=0;i<npts;i++) ai[thread_num][i]=tdata[i];
+  robust_meanrms(ai[thread_num],npts);
+  m1 = ai[thread_num][npts];
+  r1 = ai[thread_num][npts+1];
   if(r1<=0.0){
     r1 = 0.0;
     m1 = 0.0;
   }
   ath = thresh*r1;
   for (i=0;i<npts;i++) {
-    if( fabs((ai[omp_get_thread_num()][i]-m1))>ath) ai[omp_get_thread_num()][i] = m1;
+    if( fabs((ai[thread_num][i]-m1))>ath) ai[thread_num][i] = m1;
   }
-  for (i=0;i<npts;i++) tdata[i]=ai[omp_get_thread_num()][i];
+  for (i=0;i<npts;i++) tdata[i]=ai[thread_num][i];
   tdata[npts] = m1;
   tdata[npts+1] = r1;
 }
@@ -271,18 +291,24 @@ void spfind(double *cdata, long int npts, float thresh, double *wt)
   double anpt,ath;
   double m1,r1,m2;
 
+  #ifdef _OPENMP
+  int thread_num = omp_get_thread_num();
+  #else
+  int thread_num=0;
+  #endif
+
   n = npts-1;
-  for (i=0;i<n;i++) ai[omp_get_thread_num()][i]=fabs(cdata[i+1]-cdata[i]) ;
-  ai[omp_get_thread_num()][n] = 0.0;
-  ai[omp_get_thread_num()][n+1] = 0.0;
-  robust_meanrms(ai[omp_get_thread_num()],n);
-  m1 = ai[omp_get_thread_num()][n];
-  r1 = ai[omp_get_thread_num()][n+1];
+  for (i=0;i<n;i++) ai[thread_num][i]=fabs(cdata[i+1]-cdata[i]) ;
+  ai[thread_num][n] = 0.0;
+  ai[thread_num][n+1] = 0.0;
+  robust_meanrms(ai[thread_num],n);
+  m1 = ai[thread_num][n];
+  r1 = ai[thread_num][n+1];
   if(r1<=0.0) r1 = 10000.0;
 
   ath = thresh*r1;
   for (i=0;i<n;i++) {
-    if( fabs(ai[omp_get_thread_num()][i]-m1)>ath){
+    if( fabs(ai[thread_num][i]-m1)>ath){
       wt[i] = -1.0;
       wt[i+1] = -1.0;
     }
@@ -298,18 +324,23 @@ void tsfind(double *ttdata, long int npts, float thresh, double *wt)
   double anpt,ath;
   double m1,r1,m2;
 
-
-  for (i=0;i<npts;i++) ai[omp_get_thread_num()][i]=ttdata[i] ;
-  ai[omp_get_thread_num()][npts] = 0.0;
-  ai[omp_get_thread_num()][npts+1] = 0.0;
-  robust_meanrms(ai[omp_get_thread_num()],npts);
-  m1 = ai[omp_get_thread_num()][npts];
-  r1 = ai[omp_get_thread_num()][npts+1];
+  #ifdef _OPENMP
+  int thread_num = omp_get_thread_num();
+  #else
+  int thread_num=0;
+  #endif
+
+  for (i=0;i<npts;i++) ai[thread_num][i]=ttdata[i] ;
+  ai[thread_num][npts] = 0.0;
+  ai[thread_num][npts+1] = 0.0;
+  robust_meanrms(ai[thread_num],npts);
+  m1 = ai[thread_num][npts];
+  r1 = ai[thread_num][npts+1];
   if(r1<=0.0) r1 = 10000.0;
 
   ath = thresh*r1;
   for (i=0;i<npts;i++) {
-    if( fabs((ai[omp_get_thread_num()][i]-m1))>ath){
+    if( fabs((ai[thread_num][i]-m1))>ath){
       wt[i] = -1.0;
     }
   }
@@ -355,6 +386,7 @@ void cleanit(float *data, int nchans, long int nadd)
   float thresh;
   double lg,lgr, an, atemp;
 
+  int thread_num=0;
 
 // get some pre-cleaning statistics
 
@@ -362,32 +394,35 @@ void cleanit(float *data, int nchans, long int nadd)
 
   #pragma omp parallel for private(channum,inc,isame,ii,jj) schedule(dynamic)
   for (channum=0; channum<nchans; channum++) {
+    #ifdef _OPENMP
+    thread_num = omp_get_thread_num();
+    #endif
     inc = coff[channum];
   /* Select the correct channel */
-    for (ii = 0, jj = channum; ii < nadd; ii++, jj += nchans) chandata[omp_get_thread_num()][ii] = (double) data[jj];
+    for (ii = 0, jj = channum; ii < nadd; ii++, jj += nchans) chandata[thread_num][ii] = (double) data[jj];
     for (ii = 0; ii<nadd; ii++) {
-      in[omp_get_thread_num()][ii][0] = chandata[omp_get_thread_num()][ii];
-      in[omp_get_thread_num()][ii][1] = 0.0;
+      in[thread_num][ii][0] = chandata[thread_num][ii];
+      in[thread_num][ii][1] = 0.0;
     }
     if( pcl == 0 ) {
       isame = 0;
     } else {
-      isame = all_samef(chandata[omp_get_thread_num()],nadd);
+      isame = all_samef(chandata[thread_num],nadd);
     }
     if( isame == 0) {
-      fftclean(in[omp_get_thread_num()],out[omp_get_thread_num()],nadd,inc); // clean some periodic RFIs
-      for (ii = 0; ii<nadd; ii++) chandata[omp_get_thread_num()][ii] = in[omp_get_thread_num()][ii][0];
-      tsclip(chandata[omp_get_thread_num()],nadd,sthresh); // clean some spiky RFIs
-      mspec[channum]=chandata[omp_get_thread_num()][nadd];
-      rspec[channum]=chandata[omp_get_thread_num()][nadd+1];
+      fftclean(in[thread_num],out[thread_num],nadd,inc); // clean some periodic RFIs
+      for (ii = 0; ii<nadd; ii++) chandata[thread_num][ii] = in[thread_num][ii][0];
+      tsclip(chandata[thread_num],nadd,sthresh); // clean some spiky RFIs
+      mspec[channum]=chandata[thread_num][nadd];
+      rspec[channum]=chandata[thread_num][nadd+1];
     }
     else {
-      mspec[channum]=chandata[omp_get_thread_num()][0];
+      mspec[channum]=chandata[thread_num][0];
       rspec[channum]=0.0;
     }
     vspec[channum] = rspec[channum]*rspec[channum];
     for (ii = 0, jj = channum; ii < nadd; ii++, jj += nchans)
-      data[jj] = (float) chandata[omp_get_thread_num()][ii];
+      data[jj] = (float) chandata[thread_num][ii];
   }
 // get some post- fft and timeseries cleaning statistics
 
@@ -449,22 +484,25 @@ void cleanit(float *data, int nchans, long int nadd)
 // Try clipping some channels in individual samples
   #pragma omp parallel for private(t,c,nxc) schedule(dynamic)
   for (t=0; t<nadd; t++){
+      #ifdef _OPENMP
+      thread_num = omp_get_thread_num();
+      #endif
       nxc = t*nchans;
-      for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][c] = data[nxc+c];
-      spclip(chandata[omp_get_thread_num()],nchans,clipthresh);
-      for (c=0; c<nchans; c++) data[nxc+c] = chandata[omp_get_thread_num()][c];
+      for (c=0; c<nchans; c++) chandata[thread_num][c] = data[nxc+c];
+      spclip(chandata[thread_num],nchans,clipthresh);
+      for (c=0; c<nchans; c++) data[nxc+c] = chandata[thread_num][c];
   }
 // Now some timeseries cleaning
   an = (double)nadd;
   for (t=0; t<nadd; t++){
       nxc = t*nchans;
-      chandata[omp_get_thread_num()][t]=0.0;
-      for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][t] = chandata[omp_get_thread_num()][t]+data[nxc+c];
-      chandata[omp_get_thread_num()][t]=chandata[omp_get_thread_num()][t]/an;
+      chandata[thread_num][t]=0.0;
+      for (c=0; c<nchans; c++) chandata[thread_num][t] = chandata[thread_num][t]+data[nxc+c];
+      chandata[thread_num][t]=chandata[thread_num][t]/an;
   }
   for (i=0; i<nadd; i++) wt[i]=+1.0;
-  tsfind(chandata[omp_get_thread_num()],nadd,sthresh,wt);
-  spfind(chandata[omp_get_thread_num()],nadd,sthresh,wt);
+  tsfind(chandata[thread_num],nadd,sthresh,wt);
+  spfind(chandata[thread_num],nadd,sthresh,wt);
   for (t=0; t<nadd; t++){
       nxc = t*nchans;
       if(wt[t] < 0.0){
@@ -481,8 +519,8 @@ void cleanit(float *data, int nchans, long int nadd)
       isame = 0;
       for(jj=0; jj<nvar; jj++){
         ii = nxc + jj*kk;
-        for (c=0; c<kk; c++) chandata[omp_get_thread_num()][c] = data[ii+c];
-        if(all_samef(chandata[omp_get_thread_num()],kk)>0) isame = isame + 1;
+        for (c=0; c<kk; c++) chandata[thread_num][c] = data[ii+c];
+        if(all_samef(chandata[thread_num],kk)>0) isame = isame + 1;
       }
       if(isame >= (int)(0.8*nvar)){
         for (c=0; c<nchans; c++) data[c+nxc] = 0.0;
@@ -500,10 +538,10 @@ void cleanit(float *data, int nchans, long int nadd)
            rspec[t] = 0.0;
           }
           else{
-           for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][c] = data[nxc+c];
-           robust_meanrms(chandata[omp_get_thread_num()],nchans);
-           //simple_meanrms(chandata[omp_get_thread_num()],nchans);
-           rspec[t] = chandata[omp_get_thread_num()][nchans+1]*chandata[omp_get_thread_num()][nchans+1];
+           for (c=0; c<nchans; c++) chandata[thread_num][c] = data[nxc+c];
+           robust_meanrms(chandata[thread_num],nchans);
+           //simple_meanrms(chandata[thread_num],nchans);
+           rspec[t] = chandata[thread_num][nchans+1]*chandata[thread_num][nchans+1];
           }
       }
       t = nadd;
@@ -513,15 +551,15 @@ void cleanit(float *data, int nchans, long int nadd)
     }
     nvar = nadd/nsect;
     for (jj=0; jj<nsect; jj++){
-      for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][c] = 0.0;
+      for (c=0; c<nchans; c++) chandata[thread_num][c] = 0.0;
       for (t=0; t<nvar; t++){
         nxc = (jj*nvar+t)*nchans;
-        for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][c] = chandata[omp_get_thread_num()][c]+data[nxc+c];
+        for (c=0; c<nchans; c++) chandata[thread_num][c] = chandata[thread_num][c]+data[nxc+c];
       }
-      for (c=0; c<nchans; c++) chandata[omp_get_thread_num()][c] = chandata[omp_get_thread_num()][c]/nvar;
-      //robust_meanrms(chandata[omp_get_thread_num()],nchans);
-      simple_meanrms(chandata[omp_get_thread_num()],nchans);
-      an = chandata[omp_get_thread_num()][nchans+1]*chandata[omp_get_thread_num()][nchans+1];
+      for (c=0; c<nchans; c++) chandata[thread_num][c] = chandata[thread_num][c]/nvar;
+      //robust_meanrms(chandata[thread_num],nchans);
+      simple_meanrms(chandata[thread_num],nchans);
+      an = chandata[thread_num][nchans+1]*chandata[thread_num][nchans+1];
       if(an<=0.0) an=10000000.0;
       for (t=0; t<nvar; t++){
         nxc = (jj*nvar+t)*nchans;
diff --git a/src/rficlean_data.c b/src/rficlean_data.c
index 23950e8..4676a5f 100644
--- a/src/rficlean_data.c
+++ b/src/rficlean_data.c
@@ -67,16 +67,19 @@ void rficlean_data(FILE *input, FILE *output)
   psrifs = (long int *) malloc(sizeof(long int)*(naddt));
 
   #ifdef _OPENMP
+  int max_threads = omp_get_max_threads();
   if (fftw_init_threads()) fftw_plan_with_nthreads(omp_get_max_threads());
+  #else
+  int max_threads = 1;
   #endif
 
-  chandata = (double **) malloc(omp_get_max_threads() * sizeof(double *));
-  in = (fftw_complex**) malloc (omp_get_max_threads() * sizeof(fftw_complex*));
-  out = (fftw_complex**) malloc (omp_get_max_threads() * sizeof(fftw_complex*));
-  fplan = (fftw_plan*) malloc (omp_get_max_threads() * sizeof(fftw_plan));
-  bplan = (fftw_plan*) malloc (omp_get_max_threads() * sizeof(fftw_plan));
-  ai = (double **) malloc(omp_get_max_threads() * sizeof(double *));
-  for (i=0; i<omp_get_max_threads(); i++) {
+  chandata = (double **) malloc(max_threads * sizeof(double *));
+  in = (fftw_complex**) malloc (max_threads * sizeof(fftw_complex*));
+  out = (fftw_complex**) malloc (max_threads * sizeof(fftw_complex*));
+  fplan = (fftw_plan*) malloc (max_threads * sizeof(fftw_plan));
+  bplan = (fftw_plan*) malloc (max_threads * sizeof(fftw_plan));
+  ai = (double **) malloc(max_threads * sizeof(double *));
+  for (i=0; i<max_threads; i++) {
     chandata[i] = (double *) malloc(nsize * sizeof(double));
     ai[i] = (double *) malloc(nsize * sizeof(double));
     in[i] = fftw_malloc ((naddt+5) * sizeof(fftw_complex));
@@ -104,9 +107,9 @@ void rficlean_data(FILE *input, FILE *output)
   nred = iter+1;
   // wmean = (double *) malloc((nred)*sizeof(double));
   // wrms = (double *) malloc((nred)*sizeof(double));
-  wmean = (double **) malloc(omp_get_max_threads() * sizeof(double *));
-  wrms = (double **) malloc(omp_get_max_threads() * sizeof(double *));
-  for (i=0; i<omp_get_max_threads(); i++) {
+  wmean = (double **) malloc(max_threads * sizeof(double *));
+  wrms = (double **) malloc(max_threads * sizeof(double *));
+  for (i=0; i<max_threads; i++) {
     wmean[i] = (double *) malloc(nred * sizeof(double));
     wrms[i] = (double *) malloc(nred * sizeof(double));
   }
@@ -317,7 +320,7 @@ void rficlean_data(FILE *input, FILE *output)
   free(wt);
   free(coff);
 
-  for (i=0; i<omp_get_max_threads(); i++) {
+  for (i=0; i<max_threads; i++) {
     fftw_free (in[i]);
     fftw_free (out[i]);
     fftw_destroy_plan ( fplan[i] );
@@ -337,12 +340,13 @@ void rficlean_data(FILE *input, FILE *output)
   free (chandata);
 
   FILE *file_timing = fopen("timing.csv", "a+");
-  fprintf(file_timing, "%d,%d,%d,%.3f,%.3f\n",
+  fprintf(file_timing, "%s,%d,%d,%d,%.3f,%.3f\n",
+    inpfile,
     naddt,
-    omp_get_max_threads(),
+    max_threads,
     time(NULL) - s,
-    ((double) clock() - t_global) / (CLOCKS_PER_SEC) / (double) omp_get_max_threads(),
-    ((double) t_function_sum / (i_function * CLOCKS_PER_SEC)) / (double) omp_get_max_threads()
+    ((double) clock() - t_global) / (CLOCKS_PER_SEC) / (double) max_threads,
+    ((double) t_function_sum / (i_function * CLOCKS_PER_SEC)) / (double) max_threads
   );
   // fprintf(file_timing, "%s,%d,%d\n", time_section, -1, time(NULL) - s);
   fclose(file_timing);

From 9ddd06b1dcd0791b16b92a235f6cb4f90bff9dda Mon Sep 17 00:00:00 2001
From: macrocosme <vohl@astron.nl>
Date: Mon, 11 Jan 2021 19:18:08 +0100
Subject: [PATCH 4/6] Adding ifdefs _OPENMP to includes

---
 src/cleanit.c       | 2 ++
 src/rficlean_data.c | 4 +++-
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/cleanit.c b/src/cleanit.c
index 839efef..0d2ae6d 100644
--- a/src/cleanit.c
+++ b/src/cleanit.c
@@ -11,7 +11,9 @@
 #include <math.h>
 #include "rficlean.h"
 #include <fftw3.h>
+#ifdef _OPENMP
 #include <omp.h>
+#endif
 
 //--------------------------------------------------------------
 int all_samef(double *arr, long int np)
diff --git a/src/rficlean_data.c b/src/rficlean_data.c
index 4676a5f..f7c5bfd 100644
--- a/src/rficlean_data.c
+++ b/src/rficlean_data.c
@@ -12,8 +12,10 @@
 #include <math.h>
 #include "rficlean.h"
 #include <fftw3.h>
-#include <omp.h>
 #include <time.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
 
 void rficlean_data(FILE *input, FILE *output)
 {

From 48b9f9235a586da029d421b6237fb9d592c49585 Mon Sep 17 00:00:00 2001
From: macrocosme <vohl@astron.nl>
Date: Wed, 17 Feb 2021 10:25:56 +0100
Subject: [PATCH 5/6] Fix double pointer, missing bracket, timing print

---
 include/rficlean.h  | 2 +-
 src/cleanit.c       | 1 +
 src/rficlean_data.c | 1 -
 3 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/include/rficlean.h b/include/rficlean.h
index e4916fe..d55ae05 100644
--- a/include/rficlean.h
+++ b/include/rficlean.h
@@ -9,7 +9,7 @@ long byte_offset;
 int headerless,obits,iflip,nred, iwhite,nints,nsub,maxhist,numhist,nsubchans,tnorm,fnorm,igetstat,nharm,pcl,nsect,zerodm,bl_start,nblocks;
 char inpfile[128], outfile[128], gminfofile[128],gmhdrfile[128],psfile[100];
 FILE *input, *output, *gminfo, *gmhdr;
-double tempra,tempdec, psrf, psrfdf, *last_mspec, *wrms, *wmean, *ai;
+double tempra,tempdec, psrf, psrfdf, *last_mspec, **wrms, **wmean, **ai;
 float forcefthresh,fthresh,rthresh,sthresh,last_tvar,last_fvar,chanfrac,sampfrac,clipthresh;
 float *fftstat, *chanstat, *predist, *xpredist, *postdist, *xpostdist, *finaldist, *xfinaldist ;
 double *tfvar, *tfmean,meanvar,rmsvar ;
diff --git a/src/cleanit.c b/src/cleanit.c
index 1260ee1..15cae67 100644
--- a/src/cleanit.c
+++ b/src/cleanit.c
@@ -525,6 +525,7 @@ void cleanit(float *data, int nchans, long int nadd)
         if(wt[t] < 0.0){
           for (c=0; c<nchans; c++) data[nxc+c] = mspec[c];
         }
+    }
   }
 
  // sometimes gpt results might need additional checks
diff --git a/src/rficlean_data.c b/src/rficlean_data.c
index 6ea3ea2..36eb608 100644
--- a/src/rficlean_data.c
+++ b/src/rficlean_data.c
@@ -350,5 +350,4 @@ void rficlean_data(FILE *input, FILE *output)
   free (chandata);
 
   printf (" Freed the buffers. All done!\n\n ");
-  printf ("%d seconds\n\n", time(NULL) - s);
 }

From 602f65e35fb7a68afe0dc13501e901902963149e Mon Sep 17 00:00:00 2001
From: macrocosme <vohl@astron.nl>
Date: Wed, 31 Mar 2021 19:22:36 +0200
Subject: [PATCH 6/6] Add missing thread_num as omp private

---
 src/cleanit.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/cleanit.c b/src/cleanit.c
index 15cae67..2318121 100644
--- a/src/cleanit.c
+++ b/src/cleanit.c
@@ -394,7 +394,7 @@ void cleanit(float *data, int nchans, long int nadd)
 
   for (channum=0; channum<nchans; channum++) coff[channum]=nadd*channum ;
 
-  #pragma omp parallel for private(channum,inc,isame,ii,jj) schedule(dynamic)
+  #pragma omp parallel for private(channum,inc,isame,ii,jj,thread_num) schedule(dynamic)
   for (channum=0; channum<nchans; channum++) {
     #ifdef _OPENMP
     thread_num = omp_get_thread_num();
@@ -497,7 +497,7 @@ void cleanit(float *data, int nchans, long int nadd)
 
 // Try clipping some channels in individual samples
   if(rfiSclip){
-    #pragma omp parallel for private(t,c,nxc) schedule(dynamic)
+    #pragma omp parallel for private(t,c,nxc,thread_num) schedule(dynamic)
     for (t=0; t<nadd; t++){
         #ifdef _OPENMP
         thread_num = omp_get_thread_num();