Почему декомпозиция LU с использованием Parallel.For не работает?

Question

Я пытаюсь решить разложение LU с помощью алгоритма Дулиттла - согласно этому документу. Без распараллеливания код работает нормально. Тем не менее, я хотел бы, чтобы этот код выполнялся параллельно - пытаясь создать параллельный внешний и два внутренних цикла. К сожалению, я все еще что-то упускаю. Если я сначала попытался запустить внешний цикл параллельно, я получил плохой результат.

Я пытался определить, где может быть столкновение. Я запер эти места потом, но все равно не получил нужный результат. Я добавил их в скопированный код в качестве комментариев. Что я делаю не так, мне нужно заблокировать другие места?

Какова правильная реализация внешнего цикла?

Как бы выглядели внутренние петли?

Заранее спасибо.

Реализация алгоритма (последовательная)

        //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 < n; i++)
        {
            for (int j = i; j < n; j++)
                upper[i, j] = arr[i, j];
            for (int j = i + 1; j < n; j++)
                lower[j, i] = arr[j, i];
            lower[i, i] = 1;
        }

        for(int i=0; i<n; i++)
        {
            for (int j = i; j < n; j++)
            {
                for (int k = 0; k < i; k++)
                { 
                    upper[i, j] = upper[i, j] - (lower[i, k] * upper[k, j]);
                }
            }

            for (int j = i + 1; j < n; j++)
            {
                for (int k = 0; k < i; k++)
                {
                    lower[j, i] = lower[j, i] - (lower[j, k] * upper[k, i]);
                }
                    lower[j, i] = lower[j, i] / upper[i, i];
            }
        }

Реализация алгоритма (параллельная)

        //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 < n; i++)
        {
            for (int j = i; j < n; j++)
                upper[i, j] = arr[i, j];
            for (int j = i + 1; j < n; j++)
                lower[j, i] = arr[j, i];
            lower[i, i] = 1;
        }

        //making outer loop parallel
        Parallel.For(0, n, (i, state) =>
        {
            //possibly make this loop parallel also
            for (int j = i; j < n; j++)
            {
                for (int k = 0; k < 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 < n; j++)
            {
                for (int k = 0; k < 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];
            }
        });

последовательный правильный результат

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

параллельный плохой результат

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 -->BAD     -2 -1  1

РЕДАКТИРОВАТЬ Я попытался заблокировать все подходы к полям с помощью одной блокировки. Я осознаю, что таким образом теряю всю параллелизацию. Однако я хотел добиться хотя бы правильного результата, к сожалению, безуспешно.

 static object mylock = new object();
            //making outer loop parallel
            Parallel.For(0, n, (i, state) =>
            {
                for (int j = i; j < n; j++)
                {
                    for (int k = 0; k < i; k++)
                    {
                        lock (mylock)
                        {
                            upper[i, j] = upper[i, j] - (lower[i, k] * upper[k, j]);
                        }
                    }
                }

                for (int j = i + 1; j < n; j++)
                {
                    for (int k = 0; k < 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];
                    }
                }
            });

Answer 1

Параллельные циклы записывают в один и тот же массив, верно?

upper[i, j] = upper[i, j] - (lower[i, k] * upper[k, j]);

Но не определено, когда какой цикл будет записывать в позицию в вашем массиве. Таким образом, два цикла НЕ будут записывать в один и тот же индекс, а будут читать из индекса, где другой цикл, возможно, уже записал. Вы не можете парализовать свой алгоритм таким образом.