title: Segmentation fault on parallelization with OpenMP(Parallel for)
body: <p>When I run my code without "parallel for" it works perfectly . But when I use "Parallel for"
  , I receive Core dumped error . The problem occurs in calling operators .
  The operators works perfectly when I comment the parallel for syntax .</p>
  
  <p>This is the part of my code :</p>
  
  <blockquote>
  <pre><code>#include &lt;iostream&gt;
  #include &lt;vector&gt;
  #include &lt;cmath&gt;
  #include &lt;omp.h&gt;
  #include &lt;Eigen/Core&gt;
  #include &lt;Eigen/Dense&gt;
  
  #include "NEGF.h"
  
  using namespace std;
  using namespace Eigen;
  
  void NEGF::PHONON(){
  
      cout &lt;&lt; endl &lt;&lt; "*********************************" &lt;&lt; endl ;
      cout         &lt;&lt; "PHONON calculations is running ...   " &lt;&lt; endl ;
      cout         &lt;&lt; "*********************************" &lt;&lt; endl ;
  
      int Nhw0, ii, kk, iE,n_hw, iK,iKx, iteta, iq , XY, YX;
    double Kp ,Kt,Ktx , V_Value,gamma,dteta,Nw,KBT;
      Nhw0=floor(hw0/dE);
  
      vector&lt; complex&lt;double&gt; &gt; Gn_ems(nZ*nZ*nK);
      vector&lt; complex&lt;double&gt; &gt; Gn_abs(nZ*nZ*nK);
      vector&lt; complex&lt;double&gt; &gt; Vpop(nZ*nZ*nK);
      MatrixXcd M_Gamma(nZ,nZ) ;
      MatrixXcd M_Self(nZ,nZ) ;
      vector&lt;double&gt; teta(40);
  
      gamma=q*hw0/(0.5/epsAC -0.5/epsDC);
      KBT=Temp*0.026/300;
      Nw=1/(exp(hw0/(KBT))-1);
  
      for(ii=0;ii&lt;teta.size();ii++){
        teta[ii]= M_PI*ii/(teta.size()-1) ;
      }
  
      dteta=teta[1]-teta[0];
  #pragma omp parallel for
    for(iKx=0;iKx&lt;nK;iKx++){
      Ktx=KX[iKx];
      for(iE=0;iE&lt;nE;iE++){
          M_Gamma.setZero(nZ,nZ);
          M_Self.setZero(nZ,nZ);
  
          for(iteta=0;iteta&lt;teta.size();iteta++){
              for(iK=0;iK&lt;nK;iK++){
                Kt=KX[iK];
                if(iE+Nhw0 &lt; nE){
                    n_hw=iE+Nhw0 ;
                    XY=1;
                    operator()(Gn,Gn_ems,iK,n_hw,XY);
                }
                else{
                    XY=0;
                    operator()(Gn,Gn_ems,iK,iE,XY);
                }
  
                if(iE-Nhw0&gt;=0){
                    n_hw=iE-Nhw0 ;
                    YX=1;
                    operator()(Gn,Gn_abs,iK,n_hw,YX);
                }
                else{
                    YX=0;
                    operator()(Gn,Gn_abs,iK,iE,YX);
                }
                      Kp=sqrt(pow(Ktx,2)+pow(Kt,2)-2*Kt*Ktx*cos(teta[iteta])) ;
                      for(ii=0;ii&lt;nZ;ii++){
                          for(kk=0;kk&lt;nZ;kk++){
                              V_Value=exp(-sqrt(pow(Kp,2)+pow(qD,2))*abs(ii-kk)*dZ)*
                                      (1-pow(qD,2)*abs(ii-kk)*dZ/(2*sqrt(pow(Kp,2)+pow(qD,2)))-
                                              pow(qD,2)/(2*(pow(Kp,2)+pow(qD,2))))/(pow(Kp,2)+pow(qD,2));
                              operator()(Vpop,0,ii,kk,iK,V_Value);
                          }
                      }
  
                }
              for(iq=0;iq&lt;nK;iq++){
                for(ii=0;ii&lt;nZ;ii++){
                    for(kk=0;kk&lt;nZ;kk++){
                        M_Gamma(ii,kk)=M_Gamma(ii,kk)+ (2*M_PI*dK*KX[iq]*dteta*gamma*M_PI/pow(M_PI*2,2))
                                *operator()(Vpop,ii,kk,iq)*(operator()(Gn_abs,ii,kk,iq)+operator()(Gn_ems,ii,kk,iq));
                    }
                }
              }
              }
          M_Self=-Xi*(M_Gamma-M_Gamma.adjoint().eval())/2;
          operator()(GAMPOP,M_Gamma,iKx,iE);
          operator()(POPSELF,M_Self,iKx,iE);
  
          }
    }
  }
  </code></pre>
  </blockquote>