title: Why does LU decomposition using Parallel.For not work?
body: <p>I am trying to solve LU decomposition with the Doolittle Algorithm – according to this <a href="http://www.engr.colostate.edu/~thompson/hPage/CourseMat/Tutorials/CompMethods/doolittle.pdf" rel="nofollow noreferrer">document</a>. Without parallelization, code works fine. However, I would like to make this code run in parallel - trying to make a parallel outer and two inner loops. Unfortunately I am still missing something. If I tried to first make the outer loop run in parallel, I received a bad result.</p>
  
  <p>I tried to detect where there might be a collision. I locked those places afterwards, but I still did not receive the right result. I added them to the copied code as comments. What am I doing wrong, do I need lock out other places?
  <p><strong>What´s the right implementation of outer loop?</strong> 
  <p><strong>How would inner loops look like?</strong> 
  <p>Thank you in advance.</p>
  
  <p><br>
  <br>
  <em>Implementation of algorithm (sequential)</em></p>
  
  <pre><code>        //upper triangle
          var upper = new double[arr.GetLength(0), arr.GetLength(0)];
          //lower triangle
          var lower = new double[arr.GetLength(0), arr.GetLength(0)];
  
          //field initialization
          for (int i = 0; i &lt; n; i++)
          {
              for (int j = i; j &lt; n; j++)
                  upper[i, j] = arr[i, j];
              for (int j = i + 1; j &lt; n; j++)
                  lower[j, i] = arr[j, i];
              lower[i, i] = 1;
          }
  
          for(int i=0; i&lt;n; i++)
          {
              for (int j = i; j &lt; n; j++)
              {
                  for (int k = 0; k &lt; i; k++)
                  { 
                      upper[i, j] = upper[i, j] - (lower[i, k] * upper[k, j]);
                  }
              }
  
              for (int j = i + 1; j &lt; n; j++)
              {
                  for (int k = 0; k &lt; i; k++)
                  {
                      lower[j, i] = lower[j, i] - (lower[j, k] * upper[k, i]);
                  }
                      lower[j, i] = lower[j, i] / upper[i, i];
              }
          }
  </code></pre>
  
  <p><em>Implementation of algorithm (parallel)</em>
  </p>
  
  <pre><code>        //upper triangle
          var upper = new double[arr.GetLength(0), arr.GetLength(0)];
          //lower triangle
          var lower = new double[arr.GetLength(0), arr.GetLength(0)];
  
          //field initialization
          for (int i = 0; i &lt; n; i++)
          {
              for (int j = i; j &lt; n; j++)
                  upper[i, j] = arr[i, j];
              for (int j = i + 1; j &lt; n; j++)
                  lower[j, i] = arr[j, i];
              lower[i, i] = 1;
          }
  
          //making outer loop parallel
          Parallel.For(0, n, (i, state) =&gt;
          {
              //possibly make this loop parallel also
              for (int j = i; j &lt; n; j++)
              {
                  for (int k = 0; k &lt; i; k++)
                  {
                      //lower[i,k] is it potential problem?
                      /*
                       * I tried this solution
                       * double a;
                       * lock(lowerLock){
                       *   a = lower[i,k];
                       * }
                       * upper[i, j] = upper[i, j] - (a * upper[k, j])
                       */
                      upper[i, j] = upper[i, j] - (lower[i, k] * upper[k, j]);
                  }
              }
  
              //possibly make this loop parallel also
              for (int j = i + 1; j &lt; n; j++)
              {
                  for (int k = 0; k &lt; i; k++)
                  {
                      //upper [k,i] is it potential problem?
                      /*
                       * I tried this solution
                       * double b;
                       * lock(upperLock){
                       *   b = upper[k, i];
                       * }
                       * lower[j, i] = lower[j, i] - (lower[j, k] * b);
                       */
                      lower[j, i] = lower[j, i] - (lower[j, k] * upper[k, i]);
                  }
                      lower[j, i] = lower[j, i] / upper[i, i];
              }
          });
  </code></pre>
  
  <p><em>sequential right result</em></p>
  
  <pre><code>Concatenation  Upper triangle  Lower triangle
   2 -1 -2       2 -1 -2          1  0  0
  -2  4 -1       0  4 -1         -2  1  0
  -2 -1  3       0  0  3         -2 -1  1
  </code></pre>
  
  <p><em>parallel bad result</em></p>
  
  <pre><code>Concatenation  Upper triangle    Lower triangle
   2 -1 -2       2 -1 -2             1  0  0
  -2  4 -1       0  4 -1            -2  1  0
  -2 -1  3       0  0 10 --&gt;BAD     -2 -1  1
  </code></pre>
  
  <p><strong>EDIT</strong>
  I tried to lock all the approaches to the fields with one lock. I am aware of loosing all parallelization in this way. However, I wanted to achieve at least the right result, unfortunately without success.</p>
  
  <pre><code> static object mylock = new object();
              //making outer loop parallel
              Parallel.For(0, n, (i, state) =&gt;
              {
                  for (int j = i; j &lt; n; j++)
                  {
                      for (int k = 0; k &lt; i; k++)
                      {
                          lock (mylock)
                          {
                              upper[i, j] = upper[i, j] - (lower[i, k] * upper[k, j]);
                          }
                      }
                  }
  
                  for (int j = i + 1; j &lt; n; j++)
                  {
                      for (int k = 0; k &lt; i; k++)
                      {
                          lock (mylock)
                          {
                              lower[j, i] = lower[j, i] - (lower[j, k] * upper[k, i]);
                          }
                      }
                      lock (mylock)
                      {
                          lower[j, i] = lower[j, i] / upper[i, i];
                      }
                  }
              });
  </code></pre>