title: OpenMP C++ Data race with zheevr
body: <p>I'm trying to parallelize a for loop in c++ with openMP. I have many matrices (of class Matrix) that need to be exponentiated with zheevr. The implementation gives a Data Race. </p>
  
  <p>The parallelized for loop is</p>
  
  <pre class="lang-cpp prettyprint-override"><code>/* in main */
  #pragma omp parallel for shared(Aexp_omp, A, Z, W) private(j) firstprivate(idt)
  for (j=0; j&lt; nMat; ++j) {
      Aexp_omp[j] = EigenMethod(A[j],-idt,&amp;(Z[j]),W[j]);
  }
  </code></pre>
  
  <p>Here, <code> Aexp_omp</code>, <code>A</code>, <code>Z</code> are all arrays of a class <code>Matrix</code>.  The exact line seems to be the second call to zheevr in <code>EigenMethod</code>. The code executes properly when the #pragma omp critical is not commented out (but putting pragma omp critical here seems to defeat the purpose of parellelizing this code).</p>
  
  <p>I do not understand where the data race is comming from since, if I remember properly, all variables decleared in a loop are private to the thread?</p>
  
  <p>Any insights as to why this is happening? Have I wrongly decleared the <code>pragma omp parallel for</code> part of the code?</p>
  
  <p>Many thanks for your help.</p>
  
  <h2>Edit: Full Code</h2>
  
  <p>To provide all the information, here is the full code (It's a little long, my apologies). It compiles on my machines with the makefile </p>
  
  <pre class="lang-cpp prettyprint-override"><code>omp_matrix: omp_matrix.cpp
          g++ -fopenmp -lgsl -lgslcblas -llapack -lblas -lm -o omp_matrix omp_matrix.cpp
  </code></pre>
  
  <p>It gives random behaviour when running due to, I suspect the data race. The code has to files, the one with main (and some extra namesspace) and the Matrix class.</p>
  
  <pre class="lang-cpp prettyprint-override"><code>/* omp_matrix.cpp */
  #include &lt;iostream&gt;
  #include &lt;omp.h&gt;
  #include &lt;cstdio&gt;
  #include &lt;ctime&gt;
  #include &lt;string&gt;
  #include &lt;cmath&gt;
  #include &lt;vector&gt;
  #include &lt;complex&gt;
  #include &lt;limits&gt;
  #include &lt;cstdlib&gt;
  
  const char Trans[]  = {'N','T','C'};
  const char UpLo[]   = {'U','L'};
  const char Jobz[]   = {'V','N'};
  const char Range[]  = {'A','V','I'};
  #ifdef __cplusplus
  extern "C" {
  #endif
  
  void zgemm_(const char* TransA, const char* TransB, const size_t* M, const size_t* N, const size_t* K, const std::complex&lt;double&gt;* alpha, const std::complex&lt;double&gt;* A, const size_t* lda, const std::complex&lt;double&gt;* B, const size_t* ldb, const std::complex&lt;double&gt;* beta, std::complex&lt;double&gt;* C, size_t* ldc);
  int zheevr_(const char* jobz, const char* range, const char* uplo, const size_t* n, std::complex&lt;double&gt;* a, size_t* lda, double* vl, double* vu, size_t* il, size_t* iu, double* abstol, size_t* m, double* w, std::complex&lt;double&gt;* z, size_t* ldz, int* isuppz, std::complex&lt;double&gt;* work, int* lwork, double* rwork, int* lrwork, int* iwork, int* liwork, int* info);
  
  #ifdef __cplusplus
  }
  #endif
  
      #include "Matrix.hpp"
  
  namespace MOs {
      //Identity Matrix
      template &lt;class T&gt; 
      inline void Identity(Matrix&lt;T&gt;&amp; I){
      size_t dim = I.GetRows();
      if (dim != I.GetColumns())
              exit(1);
      for(size_t i=0; i&lt;dim; i++)
          for(size_t j=0; j&lt;dim; j++)
                  I(i,j) = i == j ? T(1) : T(0);
      return; 
      }
  
      // The H.O. Annilation operator
      template &lt;class T&gt; 
      inline void Destroy(Matrix&lt;T&gt;&amp; a){
              size_t dim = a.GetRows();
              if (dim != a.GetColumns())
                      exit(1);
              for(size_t i=0; i&lt;dim; i++)
                      for(size_t j=0; j&lt;dim; j++)
                              a(i,j) = j == i+1 ? std::sqrt(T(j)) : T(0); 
              return; 
      }
  
      //Take the Hermitian conjugate of a complex Matrix
      inline Matrix&lt;std::complex&lt;double&gt; &gt; Dagger(const Matrix&lt;std::complex&lt;double&gt; &gt;&amp; A){
              size_t cols = A.GetColumns(), rows= A.GetRows();
              Matrix&lt;std::complex&lt;double&gt; &gt; temp(cols,rows);
              for (size_t i=0; i &lt; rows; i++)
                      for(size_t j=0; j &lt; cols; j++)
                              temp(j,i) = conj(A(i,j));
              return temp;
      }
  
      template &lt;class T&gt;
      inline Matrix&lt;T&gt; TensorProduct(const Matrix&lt;T&gt;&amp; A, const Matrix&lt;T&gt;&amp; B){
              size_t rows1 = A.GetRows(), rows2 = B.GetRows(), cols1 = A.GetColumns(), cols2 = B.GetColumns();
              size_t rows_new = rows1*rows2, cols_new = cols1*cols2, n, m;
              Matrix&lt;T&gt; temp(rows_new,cols_new);
              for(size_t i=0; i&lt;rows1; i++)
                      for(size_t j=0; j&lt;cols1; j++)
                              for(size_t p=0; p&lt;rows2; p++)
                                      for(size_t q=0; q&lt;cols2; q++){
                                              n = i*rows2+p;    
                                              m = j*cols2+q; //  0 (0 + 1)  + 1*dimb=2 + (0 + 1 )  (j*dimb+q)         
                                              temp(n,m) = A(i,j)*B(p,q);
                              }
                      return temp;
              }
  }
  
  Matrix&lt;std::complex&lt;double&gt; &gt; EigenMethod(const Matrix&lt;std::complex&lt;double&gt; &gt;&amp; Ain,  std::complex&lt;double&gt; x, Matrix&lt;std::complex&lt;double&gt; &gt;* Z=NULL, double *W=NULL );
  
  using namespace std;
  
  int main(int argc, char const *argv[]) {
      // dim is the dimension of the system
      size_t  dimQ    = 2;
      size_t  dim = dimQ*dimQ*dimQ;
  
      // Define the matrices used to build the matrix to be exponentiated
      Matrix&lt;complex&lt;double&gt; &gt; o(dimQ,dimQ), od(dimQ,dimQ), nq(dimQ,dimQ);
      MOs::Destroy(o);
      od = MOs::Dagger(o); nq=od*o;
      Matrix&lt;complex&lt;double&gt; &gt; IdentQ(dimQ, dimQ), Hd(dim, dim);
      MOs::Identity(IdentQ);
      Hd = 0.170*MOs::TensorProduct(MOs::TensorProduct(od,IdentQ),o) + 0.124* MOs::TensorProduct(MOs::TensorProduct(IdentQ,od),o);
  
      complex&lt;double&gt; idt = complex&lt;double&gt;(0.0,0.2);                 // time unit multiplied by i
  
      size_t  nMat    = 2;
      double* c   = new double[nMat];
      c[0] = 0.134;
      c[1] = -0.326;
  
      // Decalre matrices and allocate memory for them
      Matrix&lt;complex&lt;double&gt; &gt;*   A    = new Matrix&lt;complex&lt;double&gt; &gt;[nMat];
      Matrix&lt;complex&lt;double&gt; &gt;*   Aexp     = new Matrix&lt;complex&lt;double&gt; &gt;[nMat];
      Matrix&lt;complex&lt;double&gt; &gt;*   Aexp_omp = new Matrix&lt;complex&lt;double&gt; &gt;[nMat];
      for (size_t k = 0; k &lt; nMat; ++k)
              A[k] = c[k]*Hd;
  
      // METHOD 1: Serial. Exponentiate one matrix after another
      double start_s = omp_get_wtime();
      for (size_t k = 0; k &lt; nMat; ++k)
              Aexp[k] = EigenMethod(A[k],-idt);
      double end_s = omp_get_wtime();
  
      cout &lt;&lt; "Aexp[0] = \n" &lt;&lt; Aexp[0] &lt;&lt; endl;
      cout &lt;&lt; "Aexp[1] = \n" &lt;&lt; Aexp[1] &lt;&lt; endl;
  
      // METHOD 2: Parallel. Exponentiate all matrices at the same time
      // THIS DOES NOT WORK. Data race condition in EigenMetod?
      Matrix&lt;complex&lt;double&gt; &gt;* Z = new Matrix&lt;complex&lt;double&gt; &gt;[nMat];
      double** W = new double*[nMat];
      for (size_t k = 0; k &lt; nMat; ++k) {
              Z[k].initialize(dim,dim);
              W[k] = new double[dim];
      }
  
      size_t j;
      double start_p1 = omp_get_wtime( );
      #pragma omp parallel for shared(Aexp_omp,A,Z,W) private(j) firstprivate(idt)
      for (j=0; j&lt; nMat; ++j) {
              Aexp_omp[j] = EigenMethod(A[j],-idt,&amp;(Z[j]),W[j]);
      }
      double end_p1 = omp_get_wtime( );
  
      cout &lt;&lt; "Done" &lt;&lt; endl;
  
      cout &lt;&lt; "Aexp_omp[0] = \n" &lt;&lt; Aexp_omp[0] &lt;&lt; endl;
      cout &lt;&lt; "Aexp_omp[1] = \n" &lt;&lt; Aexp_omp[1] &lt;&lt; endl;
  
      cout &lt;&lt; "Serial time: " &lt;&lt; end_s - start_s &lt;&lt; ", parallel time method 2: " &lt;&lt; end_p1-start_p1 &lt;&lt; endl;
  
      delete [] c;
      delete [] A;
      delete [] Aexp;
      delete [] Aexp_omp;
  
  } /* end of int main */
  
  Matrix&lt;complex&lt;double&gt; &gt; EigenMethod(const Matrix&lt;complex&lt;double&gt; &gt;&amp; A, complex&lt;double&gt; x, Matrix&lt;complex&lt;double&gt; &gt;* Z, double *W){ 
  
      bool returnZ(Z!=NULL), returnW(W!=NULL);                            // flags to see if the function received valid Z and W adresses
  
      //if (MOs::IsHermitian(A)==0) UFs::MyError("Routine ExpM::EigenMethod: The input Matrix is not Hermitian.");
      size_t N = A.GetRows(), LDA=A.GetLD(), IL, IU, M, LDZ=N;
      double VL, VU, ABSTOL=numeric_limits&lt;double&gt;::min(); 
      complex&lt;double&gt;*    WORK;
      double*         RWORK;
      int*            IWORK;
  
      //finding the optimal work sizes
      int LWORK=-1, LRWORK=-1, LIWORK=-1, ok=0;
      WORK    = new std::complex&lt;double&gt;[1];
      RWORK   = new double[1];
      IWORK   = new int[1];
  
      // ZHEEV computes all eigenvalues and, optionally, eigenvectors of a complex Hermitian Matrix A.
      // Jobz[0] = V i.e. compute eigenvalues and eigen-vectors, Range[0] = A i.e. all eigen values will be found, UpLo[1] = L i.e. lower triangle of matrix is stored
      zheevr_ (&amp;Jobz[0],&amp;Range[0],&amp;UpLo[1],&amp;N,NULL,&amp;LDA,&amp;VL,&amp;VU,&amp;IL,&amp;IU,&amp;ABSTOL,&amp;M,NULL,NULL,&amp;LDZ,NULL,WORK,&amp;LWORK,RWORK,&amp;LRWORK,IWORK,&amp;LIWORK,&amp;ok);
  
      //Setting the optimal workspace
      LWORK   = (int)real(WORK[0]);
      LRWORK  = (int)RWORK[0];
      LIWORK  = IWORK[0];
  
      delete [] WORK; 
      delete [] RWORK;
      delete [] IWORK;
  
      //Running the algorithim 
      if (!returnZ) Z = new Matrix&lt;std::complex&lt;double&gt; &gt;(LDZ,LDZ);
      if (!returnW) W = new double[LDA];
  
      Matrix&lt;std::complex&lt;double&gt; &gt; temp(A);
  
      int *ISUPPZ;
      ISUPPZ  = new int[2*N]; 
      WORK    = new std::complex&lt;double&gt;[LWORK];
      RWORK   = new double[LRWORK];
      IWORK   = new int[LIWORK];
      //#pragma omp critical
      //{
      // THIS LINE IS DOING SOMETHING WRONG
      zheevr_ (&amp;Jobz[0],&amp;Range[0],&amp;UpLo[1],&amp;N,temp.GetMat(),&amp;LDA,&amp;VL,&amp;VU,&amp;IL,&amp;IU,&amp;ABSTOL,&amp;M,W,Z-&gt;GetMat(),&amp;LDZ,ISUPPZ,WORK,&amp;LWORK,RWORK,&amp;LRWORK,IWORK,&amp;LIWORK,&amp;ok);   
      //}
  
      if (ok!=0) {
          cerr &lt;&lt; "Routine EigenMethod: An error occured in geting the diagonalization of A" &lt;&lt; endl;
          exit(1);
      }
  
      //Getting the Diag*Dagger(Z), this is where the exponentiation is done using the eigenvalues
      for (size_t i=0; i &lt; N; i++)
          for(size_t j=0; j &lt; N; j++)
              temp(i,j) = exp(x*W[i])*conj((*Z)(j,i));
  
      temp = (*Z)*temp;   
  
      if(!returnW) delete [] W; 
      if(!returnZ) delete Z; 
  
      delete [] ISUPPZ;
      delete [] WORK;
      delete [] RWORK;
      delete [] IWORK;
  
      return temp;
  }
  </code></pre>
  
  <p>Next is class Matrix</p>
  
  <pre class="lang-cpp prettyprint-override"><code>#ifndef Matrix_h
  #define Matrix_h
  
  enum OutputStyle {Column, List, Array};
  
  template &lt;class T&gt;
  class Matrix {  
  //friend functions get to use the private varibles of the class as well as have different classes as inputs
  template &lt;class S&gt;
  friend std::ostream&amp; operator &lt;&lt; (std::ostream&amp; output, const Matrix&lt;S&gt;&amp; A){
      for (size_t i=0; i&lt;A.rows_; i++){
          for (size_t j=0; j&lt;A.cols_; j++){
              output &lt;&lt; A.mat_[j*A.rows_ +i] &lt;&lt; "\t";
          }
          output &lt;&lt; std::endl;
      }
      return output;
  }
  
  // overloads beta*A
  template &lt;class S1, class S2&gt;
  friend Matrix&lt;S1&gt; operator*(const S2&amp; beta, const Matrix&lt;S1&gt;&amp; A){
      size_t rows= A.rows_, cols = A.cols_;
      Matrix&lt;S1&gt; temp(rows,cols);
      for (size_t i=0; i &lt; rows; i++)
          for(size_t j=0; j &lt; cols; j++)
              temp(i,j) = beta*A(i,j);
      return temp;
  }
  
  
  friend Matrix&lt;std::complex&lt;double&gt; &gt; operator*(const Matrix&lt;std::complex&lt;double&gt; &gt;&amp; A, const Matrix&lt;std::complex&lt;double&gt; &gt;&amp; B) {
      static Matrix&lt;std::complex&lt;double&gt; &gt; C;
      C.initialize(A.rows_,B.cols_);
      std::complex&lt;double&gt; alpha =1.0, beta =0.0;
      zgemm_(&amp;Trans[0], &amp;Trans[0], &amp;A.rows_, &amp;B.cols_, &amp;A.cols_, &amp;alpha, A.mat_, &amp;A.LD_, B.mat_, &amp;B.LD_, &amp;beta, C.mat_, &amp;C.LD_);
      return C;
  }
  
  
  public:
  //Construtors
  Matrix() : rows_(0), cols_(0), size_(0), LD_(0), outputstyle_(Array){ mat_=0;}
  Matrix(size_t rows, size_t cols) : rows_(rows), cols_(cols), size_(rows*cols), LD_(rows), outputstyle_(Array), mat_(new T [size_]){}
  
  Matrix(const Matrix&lt;T&gt;&amp; rhs) : rows_(rhs.rows_), cols_(rhs.cols_), size_(rhs.size_), LD_(rows_), outputstyle_(rhs.outputstyle_), mat_(new T [size_]) {
  for(size_t p = 0; p &lt; size_; p++)
      mat_[p] = rhs.mat_[p];
  }
  
  //Initialize an empty Matrix() to Matrix(size_t  rows, size_t cols)
  inline void initialize(size_t  rows, size_t cols) {
      if (rows_ != rows || cols_ != cols ) {
          if (mat_ != 0) delete [] (mat_);
          rows_   = rows;
          cols_   = cols;
          size_   = rows_*cols_;
          LD_ = rows;
          mat_    = new T[size_];
      }
  }
  
  //Destructor
  virtual ~Matrix()   {if (mat_ != 0) delete[] (mat_);}
  
  //Assignment operator
  inline Matrix&lt;T&gt;&amp; operator=(const Matrix&lt;T&gt;&amp; rhs){
      if (rows_ != rhs.rows_ || cols_ != rhs.cols_ ) {
          if (mat_ != 0) delete [] (mat_);
          rows_=rhs.rows_;
          cols_=rhs.cols_;
          size_=rows_*cols_;
          LD_=rhs.LD_;
          mat_= new T[size_];
      }
      for (size_t p=0; p &lt; size_; p++)
          mat_[p] = T(rhs.mat_[p]);
      return *this;
  }
  
  inline T&amp; operator()(size_t i, size_t j){
      if (i &gt;= rows_ || j &gt;= cols_) exit(1);
      return mat_[j*rows_+i];
  }
  
  inline T operator()(size_t i, size_t j) const{
      if (i &gt;= rows_ || j &gt;= cols_) exit(1);
      return mat_[j*rows_+i];
  }
  
  
  //overloading functions. 
  inline Matrix&lt;T&gt; operator+(const Matrix&lt;T&gt;&amp; A){
      if(rows_ != A.rows_ || cols_ != A.cols_)
          exit(1);
      Matrix&lt;T&gt; temp(rows_,cols_);
      for (unsigned int p=0; p &lt; size_; p++)      
          temp.mat_[p] = mat_[p] + A.mat_[p];
      return temp;
  }
  
  //Member Functions
  inline size_t GetColumns() const            { return cols_; }
  inline size_t GetRows() const               { return rows_; }
  inline size_t GetLD() const             { return LD_;   }
  inline size_t size() const              { return size_; }
  T* GetMat() const                   { return mat_;  }
  
  protected:
  size_t rows_, cols_, size_, LD_;
  //rows_ and cols_ are the rows and columns of the Matrix
  //size_ = rows*colums dimensions of the vector representation   
  //LD is the leading dimeonsion and for Column major order is in general eqaul to rows
  enum OutputStyle outputstyle_;
  
  T * mat_;
  
  };
  #endif /* Matrix_h */
  </code></pre>