title: Problem occur when implementing a omp parallel nested loops
body: <p>I am learning omp parallel in C++. I want to assign values to a 4D array named &quot;F&quot; with dimensions (N, M, K, K). When I implement it without parallelization, the resulting matrix of F[N-1][M-1][0:K-1][0:K-1] is</p>
  <pre><code>0.0001 600 600
  600 0.0001 600
  600 600 0.0001
  </code></pre>
  <p>However, when it is implemented using omp parallel for, the resulting matrix is</p>
  <pre><code>0.0001 75 75
  75 0.0001 75
  75 75 0.0001
  </code></pre>
  <p>where the off dignoal elements are wrong. Here are codes:</p>
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;iostream&gt;
  #include &lt;omp.h&gt;
  
  using namespace std;
  
  double**** make4Darray(int a1, int a2, int a3, int a4) {
      double**** tmp;
      tmp = (double****)malloc(a1 * sizeof(double***));
      for (int n = 0; n &lt; a1; n++) {
          tmp[n] = (double***)malloc(a2 * sizeof(double**));
          for (int m = 0; m &lt; a2; m++) {
              tmp[n][m] = (double**)malloc(a3 * sizeof(double*));
              for (int k = 0; k &lt; a3; k++) {
                  tmp[n][m][k] = (double*)malloc(a4 * sizeof(double));
              }
          }
      }
      return tmp;
  }
  
  void delet4Darray(double**** tmp, int a1, int a2, int a3, int a4) {
      for (int n = 0; n &lt; a1; n++) {
          for (int m = 0; m &lt; a2; m++) {
              for (int k = 0; k &lt; a3; k++) {
                  free(tmp[n][m][k]);
              }
              free(tmp[n][m]);
          }
          free(tmp[n]);
      }
      free(tmp);
  }
  
  int main() {
      int M = 3000;
      int N = 2000;
      int K = 3;
      int nc = 8;
  
      double A = .0001;
      double B = 0;
      double C = 0;
      double D = 0;
      double**** F;
      F = make4Darray(N, M, K, K);
  
      for (int n = 0; n &lt; N; n++) {
          for (int m = 0; m &lt; M; m++) {
              B += A;
              C -= A;
              D += A;
              for (int k = 0; k &lt; K; k++) {
                  for (int k1 = 0; k1 &lt; k; k1++) {
                      F[n][m][k][k1] = B;
                      F[n][m][k1][k] = F[n][m][k][k1];
                  }
                  F[n][m][k][k] = A;
              }
          }
      }
  
      for (int k = 0; k &lt; K; k++) {
          for (int k1 = 0; k1 &lt; K; k1++)
              cout &lt;&lt; F[N-1][M - 1][k][k1] &lt;&lt; &quot; &quot;;
          cout &lt;&lt; endl;
      }
      cout &lt;&lt; endl;
      delet4Darray(F, N, M, K, K);
  
      // Set the number of cores for parallel
      omp_set_num_threads(nc);
      A = .0001;
      B = 0;
      C = 0;
      D = 0;
      F = make4Darray(N, M, K, K);
  
  #pragma omp parallel for reduction(+:B, D) reduction(-:C)
      for (int n = 0; n &lt; N; n++) {
          for (int m = 0; m &lt; M; m++) {
              B += A;
              C -= A;
              D += A;
              for (int k = 0; k &lt; K; k++) {
                  for (int k1 = 0; k1 &lt; k; k1++) {
                      F[n][m][k][k1] = B;
                      F[n][m][k1][k] = F[n][m][k][k1];
                  }
                  F[n][m][k][k] = A;
              }
          }
      }
  
      for (int k = 0; k &lt; K; k++) {
          for (int k1 = 0; k1 &lt; K; k1++)
              cout &lt;&lt; F[N-1][M - 1][k][k1] &lt;&lt; &quot; &quot;;
          cout &lt;&lt; endl;
      }
      cout &lt;&lt; endl;
  
      return 0;
  }
  </code></pre>
  <p>Thank you very much for your kind help.</p>