Параллельные алгоритмы

Циклы

Под циклом будем понимать цикл типа while, рассматривая цикл типа for как частный случай. Будем также полагать, что с каждым циклом связана предшествующая ему некоторая инициализирующая часть Init, содержащая группу операторов. Инициализация необходима для обеспечения корректной работы цикла, после ее завершения должно выполняться предусловие цикла и стать истинным инвариант цикла.

С каждым циклом связывается одна или несколько переменных, называемых параметрами цикла, изменяющих свое значение при каждом выполнении тела цикла. В цикле for изменение параметров цикла осуществляется в заголовке цикла. В цикле while изменение параметров цикла выполняется явно в теле цикла.

Для простоты будем полагать, что тело цикла может содержать только операторы присваивания, условные операторы, составные операторы и операторы цикла. Мы полагаем, что цикл не содержит операторов вызова процедур и вызова функций с побочным эффектом.

Обозначим через V - множество переменных цикла. Сюда входят все переменные, встречающиеся в теле цикла, и параметры цикла, заданные в заголовке цикла for. Множество V представим как два непересекающихся подмножества:

$$V=V_v \bigcup V_{uv} $$

В множество $$V_v$$ входят переменные, которые могут изменять свое значение в ходе выполнения тела цикла. В это множество входят, например, параметры цикла. В множество $$V_{uv}$$ входят переменные, которые сохраняют постоянное значение на всех итерациях в ходе выполнения цикла.

Аналогично обозначим через E - множество выражений, встречающихся в цикле. Множество E представим как два непересекающихся подмножества:

$$E=E_v \bigcup E_{uv}$$

В множество $$E_{uv}$$ входят выражения, которые содержат только константы и переменные из $$V_{uv}$$. Независимо от итерации, на которой вычисляется значение этих выражений, результат вычислений будет один и тот же.

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

Если программа работает недопустимо медленно, то возникает вопрос, как ускорить ее выполнение? Основным приемом является выбор эффективного алгоритма, дающего решение исходной задачи. Классическим примером может служить задача сортировки массивов. Если необходимо сортировать большое число массивов с малым числом элементов, то вполне допустимы простые методы сортировки со сложностью $$O(n^2)$$. Когда n велико, такие методы становятся неэффективными, и следует применять методы со сложностью $$O(n \cdot log n)$$. Если элементы сортируемого массива принадлежат небольшому числу классов (двум - четырем классам), то следует применять методы, имеющие сложность O(n). Примером может служить массив персон, который нужно отсортировать по полу, разделив мужчин и женщин. Другим примером является сортировка новых учеников школы Хогвартс из романа о Гарри Поттере, где сортирующая шляпа делила учеников в зависимости от их свойств на четыре класса. Такую сортировку можно выполнить за линейное время.

Выбор наилучшего алгоритма определяется не только задачей, но и архитектурой компьютера, на котором эта задача будет решаться. Алгоритм должен наилучшим образом использовать ресурсы компьютера. Алгоритм, наилучший для компьютера с одним процессором, может быть далеко не лучшим на компьютере с множеством процессоров, поскольку не будет использовать возможности распараллеливания вычислений. Распараллеливание вычислений во многом определяется возможностью распараллеливания циклов, что и будет основной темой рассмотрения в данной главе. Но прежде давайте рассмотрим более простую тему, также связанную с циклами, и имеющую ту же цель - ускорить выполнение циклов нашей программы.

Чистка цикла

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

Рассмотрим вначале чистку выражений. Пусть expr - некоторое выражение из Е, встречающееся в теле цикла, и subexpr - некоторое его подвыражение. Если подвыражение subexpr принадлежит множеству $$E_{uv}$$, то выражение expr можно упростить. Для этого достаточно определить в Init локальную переменную loc, тип которой совпадает с типом subexpr, задать оператор присваивания loc = subexpr, затем заменить в выражении expr подвыражение subexpr переменной loc. Эту замену можно осуществить во всех вхождениях subexpr, встречающихся в выражениях из Е. Предполагается, конечно, что такая замена корректна, - subexpr не содержит вызовов функций с побочным эффектом и вхождение в expr таково, что замена subexpr на (loc) не меняет значения выражения expr.

Заметьте, что замена подвыражений локальными переменными является весьма полезным приемом, даже в тех случаях, когда вычисление подвыражения невозможно вынести из цикла. Упрощение выражений облегчает отладку, повышает надежность программы, а при наличии нескольких вхождений позволяет ускорить вычисление выражения.

Рассмотрим теперь чистку операторов. Оператор stat может быть вынесен из тела цикла и включен в конец инициализирующей части Init при условии, что все переменные этого оператора принадлежат множеству $$V_{uv}$$ и все выражения принадлежат множеству $$E_{uv}$$. Заметьте, наше определение множества $$V_{uv}$$ не исключает включения в него переменных, входящих в левые части операторов присваивания и, следовательно, получающих значения в теле цикла. Но на такие переменные накладываются дополнительные условия. Во-первых, получаемые ими значения должны быть результатом вычисления выражений из множества $$E_{uv}$$. Во-вторых, эти переменные используются в выражениях только после получения ими значений в результате присваивания.

Хороший оптимизирующий компилятор может выполнять чистку цикла. Мой анализ показал, что этого нельзя сказать как о компиляторе C#, включенном в состав Visual Studio 2010, так и о JIT компиляторе, входящем в состав Framework .Net 4.0.

Вот простой пример, демонстрирующий отсутствие автоматической чистки цикла компиляторами в C# программах:

static void Main(string[] args)
        {
            double a1 = 0.5, a2 = 0.1, a3 = -1.5;
            double x = 2.0, y = 0;
            int n  = 100000;
            DateTime start, finish;
            start = DateTime.Now;
            for (int i = 0; i < n; i++)
            {
                y = (1 + 2*(a1* Math.Pow(x, 3) +
                        a2 * Math.Pow(x, 2) + a3 * Math.Pow(x, 3)) - 
                    5 * (a1* Math.Pow(x, 3) +
                        a2 * Math.Pow(x, 2) + a3 * Math.Pow(x, 3)) +
                    3 * (a1* Math.Pow(x, 3) +
                        a2 * Math.Pow(x, 2) + a3 * Math.Pow(x, 3))) /
                    (Math.Sin(a1* Math.Pow(x, 3) + 
                        a2 * Math.Pow(x, 2) + a3 * Math.Pow(x, 3)) *
                    Math.Sin(a1* Math.Pow(x, 3) +
                        a2 * Math.Pow(x, 2) + a3 * Math.Pow(x, 3)) + 
                    Math.Cos(a1* Math.Pow(x, 3) +
                        a2 * Math.Pow(x, 2) + a3 * Math.Pow(x, 3)) *
                    Math.Cos(a1* Math.Pow(x, 3) +
                        a2 * Math.Pow(x, 2) + a3 * Math.Pow(x, 3)));
            }
            finish = DateTime.Now;
            Console.WriteLine(" y = " + y);
            Console.WriteLine("Время вычислений в тиках = " +
                (finish.Ticks - start.Ticks));
        }

Нетрудно видеть, что для рассматриваемого цикла все переменные, кроме параметра цикла, принадлежат множеству $$V_{uv}$$, а все выражения принадлежат множеству $$E_{uv}$$. В цикле многократно встречается подвыражение, принадлежащее $$E_{uv}$$, так что его вычисление можно вынести из тела цикла. Более того, сам оператор присваивания можно также вынести из тела цикла в инициализирующую часть, после чего тело цикла не будет содержать операторов, так что хороший оптимизирующий компилятор может удалить оператор цикла. Продолжая оптимизацию, можно заметить, что оператор присваивания можно заменить эквивалентным оператором

y = 1;

Однако ничего подобного не происходит. Если проанализировать IL код, построенный компилятором C# для версии Release с включенным флажком оптимизации кода, то можно видеть, что построенный код содержит несколько сотен ячеек и многократно выполняет одни и те же действия. Вот некоторый фрагмент этого кода:

IL_014c:  ldloc.0
  IL_014d:  ldloc.3
  IL_014e:  ldc.r8     3.
  IL_0157:  call       float64 [mscorlib]System.Math::Pow(float64,
                                                          float64)
  IL_015c:  mul
  IL_015d:  ldloc.1
  IL_015e:  ldloc.3
  IL_015f:  ldc.r8     2.
  IL_0168:  call       float64 [mscorlib]System.Math::Pow(float64,
                                                          float64)
  IL_016d:  mul
  IL_016e:  add
  IL_016f:  ldloc.2
  IL_0170:  ldloc.3
  IL_0171:  ldc.r8     3.
  IL_017a:  call       float64 [mscorlib]System.Math::Pow(float64,
                                                          float64)
  IL_017f:  mul
  IL_0180:  add
  IL_0181:  call       float64 [mscorlib]System.Math::Sin(float64)
  IL_0186:  mul
  IL_0187:  ldloc.0
  IL_0188:  ldloc.3
  IL_0189:  ldc.r8     3.
  IL_0192:  call       float64 [mscorlib]System.Math::Pow(float64,
                                                          float64)
  IL_0197:  mul
  IL_0198:  ldloc.1
  IL_0199:  ldloc.3
  IL_019a:  ldc.r8     2.
  IL_01a3:  call       float64 [mscorlib]System.Math::Pow(float64,
                                                          float64)
  IL_01a8:  mul
  IL_01a9:  add
  IL_01aa:  ldloc.2
  IL_01ab:  ldloc.3
  IL_01ac:  ldc.r8     3.
  IL_01b5:  call       float64 [mscorlib]System.Math::Pow(float64,
                                                          float64)

К сожалению, оптимизацию этого кода не выполняет и JIT компилятор, работающий на заключительном этапе. Приведу время, затраченное на выполнение программы:

Время вычислений как функция от длины цикла

n	10 000	100 000	1 000 000
T (тиках)	410 023	4 340 249	38 422 197

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

Распараллеливание цикла

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

Приведенное выше определение распараллеливания неконструктивно. Условия, при выполнении которых цикл может быть распараллелен, достаточно сложны. Мы ограничимся рассмотрением примеров, для каждого из которых будем пояснять, почему тот или иной цикл может или не может быть распараллелен.

Суммирование

Начнем с простейшей задачи - вычисление суммы элементов массива:

$$S=\sum_{i=0}^{n-1}x_i$$

Классический алгоритм выглядит так:

S = 0;
for(int i = 0; i < n; i++)
  S = S + x[i];

Алгоритм последовательный и не допускает распараллеливания, поскольку на каждом шаге цикла используется значение S, вычисленное на предыдущем шаге. Заметьте, цикл не допускает распараллеливания, если в теле цикла есть оператор присваивания, у которого одна и та же скалярная переменная встречается как в левой, так и в правой части оператора (S в нашем примере).

Поскольку в алгоритме используется только один цикл типа for с шагом, равным единица, то алгоритм имеет линейную временную сложность:

$$T\_Suc(n) = O(n)$$

Когда суммирование должен выполнять только один процессор, то это оптимальный по времени алгоритм. Но суммирование можно вести самыми разными способами, если к вычислению суммы привлечь несколько процессоров. В главе 1 мы уже анализировали "пирамидальный" алгоритм суммирования, допускающий распараллеливание. Вот возможная запись такого алгоритма на С#:

int[] y = new int[n];
Array.Copy(x, y, n);
int m = n;
while (m != 1)
{
    for (int i = 0, j = m - 1; i < j; i++, j--)
        y[i] = y[i] + y[j];
    m = (m + 1) / 2;
}
S = y[0];

В этом более сложном варианте алгоритма появились два цикла. Внутренний цикл for может быть распараллелен. Нетрудно видеть, что множества переменных, используемых на каждой итерации этого цикла, взаимно не пересекаются. Выполнение этого условия достаточно для распараллеливания. Внешний цикл while распараллелить нельзя, поскольку итерации "склеиваются" общей переменной m.

Для компьютера с одним процессором этот алгоритм будет хуже классического алгоритма по ряду причин:

Ему нужна дополнительная память, поскольку промежуточные результаты необходимо запоминать. Для этого введен массив y той же размерности, что и исходный массив x.

Алгоритм усложнен, - вместо одного цикла классического алгоритма, вводятся два цикла. Алгоритм менее понятен и отладка его может вызывать затруднение. Несомненно, что программисту потребуется больше времени на разработку и отладку такого алгоритма.

Алгоритм при последовательном выполнении хотя и имеет ту же временную сложность - O(n), но константа у него больше, чем у классического алгоритма, поскольку во внутреннем цикле используются три переменные с индексами, вместо одной, как в классическом алгоритме.

В чем же достоинство этого алгоритма? Одно несомненное достоинство у него есть - он допускает распараллеливание. Если запускать его на идеальном метакомпьютере с неограниченным числом процессоров, то тогда, используя n/2 процессоров, все вычисления внутреннего цикла можно выполнять параллельно.

Сложность внутреннего цикла при распараллеливании будет равна O(1), а общая сложность определяется внешним циклом и равна O(log n). В главе 1, где этот алгоритм уже анализировался, помимо ускорения в рассмотрение вводилась и такая характеристика алгоритма как эффективность. У пирамидального алгоритма эффективность низкая, поскольку для достижения максимального ускорения требуется n/2 процессоров - число, пропорциональное размерности массива. Метакомпьютеры и даже суперкомпьютеры с большим числом процессоров не всегда под рукой, но и при их наличии приходится заботиться об эффективности.

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

Разбиение массива на группы можно выполнить разными способами. Разумными представляются две стратегии. Первая состоит в том, чтобы исходный массив нарезать на р отрезков и поручить каждому процессору суммирование элементов соответствующего отрезка. Вторая стратегия основана на том, чтобы каждый процессор суммировал элементы, отстоящие на расстоянии р. Если соответствующим образом выбрать начальный элемент для каждого процессора, то это позволит вести параллельное суммирование.

Вот вариант записи сегментного алгоритма, соответствующего первой стратегии:

int count_p = Environment.ProcessorCount;          
          int[] y = new int[count_p];           
            int m = n / count_p + 1;
            int start = 0, finish = 0;
           for(int k =0; k < count_p; k++)
            {
               start = k * m; 
               finish =  (k+1)*m < n ? (k+1)*m : n;
               for (int i = start; i < finish; i++)
                    y[k] += x[i];               
            }
            S = y[0];
            for (int i = 1; i < count_p; i++)
                S += y[i];

Чуть проще шаговый алгоритм, соответствующий второй стратегии:

int[] y = new int[count_p];           
            for (int k = 0; k < count_p; k++)
            {               
                for (int i = k; i < n; i += count_p)
                    y[k] += x[i];
            }
            S = y[0];
            for (int i = 1; i < count_p; i++)
                S += y[i];

Оба эти варианта допускают параллельное выполнение тела внешнего цикла. В этом случае сложность алгоритма определяется сложностью операторов, составляющих тело этого цикла, фактически, сложностью внутреннего цикла - O(n/count_p). Ускорение для обоих вариантов равно count_p - числу процессоров, а эффективность равна единице. Таковы теоретические оценки. В последующих главах посмотрим, что можно получить на практике.

Задача суммирования крайне важна, поскольку встречается в самых разных приложениях. Она легко обобщается на случай, когда суммируются не элементы массива, а функции, зависящие от параметра i:

$$S=\sum^{n-1}_{i=0}f(i)$$

Все, что сказано о суммировании, касается и других подобных задач - нахождение произведения элементов массива, максимального или минимального элемента и других задач этого класса.

Суммирование рядов

Нахождение суммы конечного ряда принципиально не отличается от нахождения суммы элементов массива. Поговорим о том, как вычислять сумму бесконечного сходящегося ряда:

$$S=\sum^{\mathcal1}_{i=1}a(i)$$

Необходимым условием сходимости ряда является стремление $$a_i$$ к нулю, когда i стремится к бесконечности. В программировании это условие, в отличие от математики, является и достаточным условием сходимости. В математике это не так, - примером является расходящийся гармонический ряд с общим членом ряда 1/i. В мире компьютеров все дискретно, нет иррациональности, нет бесконечности, вычисления не всегда точны и могут иметь некоторую погрешность. При нахождении на компьютере суммы гармонического ряда, начиная с некоторого i*, значение общего члена станет равным нулю (так называемому машинному нулю) из-за ограниченности разрядной сетки, отводимой для хранения числа.

В программировании не ставится задача вычисления точного значения S в формуле (3.4), - достаточно вычислить это значение с некоторой точностью. В сравнении с конечными суммами дополнительная сложность в построении алгоритма состоит в том, что заранее неизвестно, сколь много членов ряда необходимо вычислить, чтобы найти сумму ряда с нужной точностью. Классический алгоритм основан на том, что суммирование прекращается, как только очередной член суммы $$a_i$$ становится по модулю меньше заданной точности $$\varepsilon$$. При этом предполагается, что выполняется условие сходимости ряда, так что все не учитываемые члены ряда будут по модулю заведомо меньше $$\varepsilon$$.

Вот пример записи такого алгоритма:

double eps = 1E-15;
           double i = 1;
           double a = 1;
           double S = 0;
           while (Math.Abs(a) > eps)
           {
               //Вычисление общего члена ряда 
               a = 1.0 / ((4 * i + 1) * (4 * i - 1));   // пример 
               S += a;
               i++;
           }

С программистской точки зрения алгоритмы 3.1 и 3.5 во многом схожи. Разница в том, что в первом случае используется цикл for, во-втором - цикл while. Из-за этого различия труднее оценить временную сложность алгоритма, поскольку нет такого естественного параметра как n в алгоритме 3.1. Число суммирований зависит как от формулы, задающей общий член ряда, так и от выбранной точности вычислений - $$\varepsilon$$.

Другой подход к вычислению суммы сходящегося ряда состоит в том, чтобы вместо точности $$\varepsilon$$ задавать N - число суммируемых элементов, сводя исходную задачу к задаче вычисления суммы конечного ряда. Иногда алгоритм усложняют, вводя дополнительный цикл while, в котором на каждом шаге цикла N увеличивается, например, вдвое. Цикл заканчивается, когда два последующих вычисленных значений суммы отличаются на величину, меньшую заданной точности $$\varepsilon$$. Оценить временную сложность такого алгоритма затруднительно.

Вернемся к рассмотрению алгоритма 3.5. По уже указанным причинам он не допускает распараллеливания. Но, как и для вычисления конечных сумм, нетрудно построить допускающую распараллеливание версию, ориентированную на p процессоров. Для конечных сумм мы рассматривали два варианта - сегментный и шаговый алгоритм.

Сегментный алгоритм требует нарезки на сегменты примерно равной длины, для этого нужно знать N - число суммируемых элементов. Поэтому этот вариант алгоритма следует применять тогда, когда заранее выбирается N.

Модификацию шагового алгоритма, допускающего распараллеливание, выполнить несложно. Для каждого процессора нужно задать начальное значение, после чего процессор будет вести суммирование, пока общий член не станет меньше заданной точности $$\varepsilon$$.

double i = 1;
            double a = 0;
            double[] y = new double[count_p];
            for (int k = 0; k < count_p; k++)
            {
                i = k+1;
                a = 1.0 / ((4 * i + 1) * (4 * i - 1)); //начальное значение               
                while (Math.Abs(a) > eps)
                {                    
                    y[k] += a;
                    i += count_p ;
                    a = 1.0 / ((4 * i + 1) * (4 * i - 1));
                }
            }
            S = y[0];
            for (int k = 1; k < count_p; k++)
                S += y[k];

Вычисление тригонометрических функций

Сходящиеся ряды появляются во многих приложениях. В частности вычисление многих функций основано на использовании их разложения в сходящийся ряд. Рассмотрим задачу вычисления значения непрерывной дифференцируемой функции, используя ее разложение в ряд Тэйлора.

С программистской точки зрения задача сводится к вычислению бесконечной суммы сходящегося ряда:

$$f(x)=\sum^{\mathcal1}_{i=0}a_i$$

При построении эффективного последовательного алгоритма, как правило, удается построить рекуррентную формулу, существенно снижающую трудоемкость вычислений

$$a_{k+1}=g(a_k)$$

Учитывая сходимость $$a_k$$, бесконечная сумма заменяется конечной суммой, когда суммирование заканчивается при условии, что $$a_k$$ по модулю становится меньше заданной точности вычислений $$\varepsilon$$. В другом варианте задаются достаточно большим значением N и вычисления прекращаются при i равном N.

Как распараллелить вычисление суммы? Понятно, что при наличии параллельно работающих P процессоров, каждый из них может вычислять свою часть суммы:

$$\sum a_i=\sum a_{i1}+\sum a_{i2}+…+\sum a_{ip}$$

Если при суммировании задавать N и выбирать его кратным P (N = k * P), то каждый из процессоров может вычислять k членов суммы. Например, первый процессор будет вычислять сумму первых k членов, второй - суммирует следующую группу из k членов и так далее. Этот алгоритм распараллеливания мы называем сегментным алгоритмом. Другой способ распараллеливания вычислений, называемый шаговым алгоритмом, состоит в том, что процессоры суммируют члены ряда, отстоящие друг от друга на расстоянии P.

Последовательный алгоритм зачастую имеет несомненные достоинства, состоящие в том, что во многих практически значимых задачах удается построить простую рекуррентную формулу для вычисления $$a_k$$ и простую формулу для вычисления начального члена суммы. Поскольку вычисление $$a_k$$ может требовать сложных вычислений, то применение простых рекуррентных соотношений позволяет существенно увеличить эффективность последовательного алгоритма. При применении сегментного алгоритма распараллеливания удается сохранить рекуррентную формулу, применяемую в последовательном алгоритме. Однако для каждого процессора необходимо вычислить начальное значение и это существенно снижает эффект, получаемый за счет распараллеливания. Для шагового алгоритма начальные значения вычисляются не столь сложно, но рекуррентная формула становится намного сложнее. Это типичная картина - распараллеливание требует жертв - усложнения алгоритма. В результате может оказаться, что привлечение дополнительных процессоров может приводить не к снижению времени вычислений, а к его росту. В подобных задач может существовать оптимальное число процессоров p*, после достижения которого время вычислений начнет возрастать.

Приведем пример, демонстрирующий указанные проблемы. В качестве функции f(x) выберем функцию ArcSin(x), для которой справедливо следующее разложение в ряд Тэйлора:

$$ArcSin(x)=x+\frac{x^3}{2 \cdot 3}+\frac{1 \cdot 3 \cdot x^5}{2 \cdot 4 \cdot 5}+…+\frac{1 \cdot 3 \cdot 5 \cdot …(2\cdot n-1)x^{2n+1}}{2 \cdot 4 \cdot 6 \cdot …(2 \cdot n)\cdot (2\cdot n+1)}+…$$

Точная формулировка задачи состоит в следующем. Дано вещественное число x такое, что $$|x| \le 1$$. Требуется найти с заданной точностью значение функции ArcSin(x), вычисляемое как сумма бесконечного сходящегося ряда (3.8). Нетрудно получить следующие соотношения:

$$a_0=x$$ $$a_i=\frac{(2i-1)!!}{(2i)!!\cdot(2i+1)}x^{2i+1}$$ $$a_{i+1}=a_i\cdot \frac{k}{(k+1)\cdot(k+2)}x^2,$$

где k=2i+1

Последовательный алгоритм построен на шаблоне, заданном алгоритмом 3.5. Отличие состоит в том, что добавляется аргумент x и текущий член вычисляется в соответствии с рекуррентной формулой, заданной соотношением (3.9):

double x2 = x * x;
            double i = 0;
            double a = x;   //начальное значение 
            double k = 0;
            S = 0;
            while (Math.Abs(a) > eps)
            {
                S += a;                
                k = 2 * i + 1;
                a *= k * k * x2 / ((k + 1) * (k + 2));
                i++;
            }

Шаговый алгоритм можно строить на основе шаблона, заданного алгоритмом 3.6. Нужно лишь корректно задать начальные значения и применить для вычисления общего члена рекуррентную формулу, связывающую i-й и i+p -й члены ряда. Для функции Arcsin(x) эта формула имеет вид:

$$a_{i+p}=a_i \cdot \frac{2i+1}{2i+2p+1}\cdot x^{2p}\cdot \prod^{p}_{k=1}(1-\frac{1}{2 \cdot (i+k)}$$

Вот код соответствующего алгоритма:

double x2 = x * x;
            double i = 0;            
            double k = 0;
            double[] y = new double[count_p];
            double a = 0;
            double x2p = x2;

            //Вычисление начальных значений 
            y[0] = x;
            for (int j = 1; j < count_p; j++)
            {
               k = 2 * i + 1;
               y[j] = y[j - 1] * k * k * x2 / ((k + 1) * (k + 2));                
               i++;
               x2p *= x2;
            }

            for(int j = 0; j < count_p; j++)
            {
                a = y[j];
                i = 0;
                double pr = 1;
                while (Math.Abs(a) > eps)
                {
                    pr = 1;                   
                    for(int r = 1; r <= count_p; r++)
                        pr *= (1- 1/(2*(i + r)));
                    a *= pr * (2*i + 1)/((2*(i +count_p) +1)) * x2p;                    
                    y[j] += a;
                    i += count_p;
                }
            }
            S = y[0];
            for (int j = 1; j < count_p; j++)
                S += y[j];

Конечно формула (3.10) существенно сложнее, чем формула (3.9) для последовательного алгоритма. Это плата за распараллеливание, которая может оказаться чрезмерной в ряде случаев. Результаты численных экспериментов для данной конкретной функции будут приведены в последующих главах.

Для сегментного алгоритма рекуррентное соотношение дается формулой (3.9), как и для последовательного алгоритма. Но вычисление начального значения для соответствующего сегмента потребует серьезных вычислительных затрат в сравнении с последовательным алгоритмом. Поэтому для задач, подобных вычислению функции ArcSin(x), использование p процессоров не даст выигрыша во времени в p раз. Не буду приводить код сегментного алгоритма, полагая, что он достаточно понятен.

Вычисление определенного интеграла

Задача вычисления определенного интеграла

$$I=\int_a^b f(x)$$

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

$$I=\sum_{i=0}^{N-1} f(x_i)\cdot dx$$

Геометрически значение интеграла задает площадь фигуры, образованной графиком подынтегральной функции f(x). Простейший численный метод - метод прямоугольников - вычисляет значение интеграла, как сумму площадей N прямоугольников, у которых основание равно dx, а высота равна значению подынтегральной функции в точке $$x_i$$. При выбранном значении N значение dx и координата $$x_i$$ вычисляются по следующим формулам:

$$dx=\frac{b-a}{N}$$ $$x_i=a+i \cdot dx$$

Если выбрать N достаточно большим, то формула (3.12) дает хорошую аппроксимацию значения интеграла. Рис. 3.1 иллюстрирует сущность метода прямоугольников.

(рис 3.1) Метод прямоугольников численного вычисления интеграла

На рисунке N равно двум и приближенное значение интеграла равно сумме площадей двух выделенных прямоугольников. Конечно же, при таком малом N трудно ожидать хорошей аппроксимации. Значение N следует существенно увеличить. Но как велико должно быть N? Это зависит как от величины интервала интегрирования, так и от поведения функции. Для осциллирующих функций N может быть очень велико.

Численный метод интегрирования предполагает итерирование по N. Это означает построение цикла, на каждой итерации которого значение N увеличивается (обычно в два раза). Если значения вычисленных сумм на соседних итерациях по модулю отличаются на величину меньшую заданной точности, то итерирование прекращается.

Метод прямоугольников хорош тем, что нетрудно написать реализацию, в которой на каждом шаге цикла по N значения функции рассчитываются только в новых точках. Сумма, вычисленная на предыдущей итерации, используется для вычисления нового значения.

Последовательный алгоритм в этом случае имеет временную сложность, заданную соотношением:

$$T(I(f))=O(N\cdot T(f))$$

Здесь N - это то конечное значение, при котором достигается требуемая точность. Цикл итерирования по N не вносит дополнительной сложности вычисления.

Давайте построим реализацию этого алгоритма. Первым делом зададим класс, описывающий подынтегральные функции:

public delegate double Integral_function(double x);

Все функции этого класса принимают один аргумент типа double и возвращают значение такого же типа.

Построим теперь класс NewIntegral, содержащий различные методы вычисления интеграла. Вот как выглядит общая часть этого класса, содержащая описание полей класса, методов свойств и конструктора объектов:

/// <summary>
    /// Вычисление определенного интеграла 
    /// разными методами
    /// </summary>
    public class NewIntegral
    {
        double a, b;    //пределы интегрирования
        Integral_function f;    //подынтегральная функция
        int p;          //число сегментов разбиения 
        double eps;     //точность вычисления
        double result;  //результат вычисления
        object locker = new object();
        double[] results;
        public double Result
        {
            get { return result; }
        }
        /// <summary>
        /// конструктор 
        /// </summary>
        /// <param name="a">нижний предел интегрирования</param>
        /// <param name="b">верхний предел интегрирования</param>
        /// <param name="p">число сегментов разбиения интервала
                                     интегрирования</param>
        /// <param name="f">подынтегральная функция</param>
        /// <param name="eps">точность вычисления интеграла</param>
        public NewIntegral(double a, double b, int p,
            Integral_function f, double eps)
        {
            this.a = a;
            this.b = b;
            this.p = p;
            this.f = f;
            this.eps = eps;
            results = new double[p];
        }

Добавим теперь в наш класс метод, реализующий последовательное вычисление интеграла по рассмотренной выше схеме прямоугольников:

/// <summary>
        /// Последовательный алгоритм
        /// </summary>
        /// <param name="a">начало отрезка интегрирования</param>
        /// <param name="b">конец отрезка интегрирования</param>
        void DefiniteIntegral(double a, double b, out double result)
        {
            int n = 2;
            double dx = (b - a) / 2;
            double S0 = 0, S = 0;
            double x = 0;
            bool success = false;
            for (int i = 0; i < n; i++)
            {
                x = a + i * dx;
                S0 += f(x);
            }
            S0 *= dx;
            while (!success)
            {
                n = 2 * n;
                dx = dx / 2;
                S = 0;
                for (int i = 1; i < n; i += 2)
                {
                    x = a + i * dx;
                    S += f(x);
                }
                S = S * dx + S0 / 2;
                if (Math.Abs(S - S0) > eps)
                    S0 = S;
                else
                    success = true;
            }
            result = S;
        }

Простой и эффективный алгоритм 3.9 полностью реализует описанную выше идею. Вначале вычисляется сумма $$S_0$$ при n, равном 2. Затем в цикле по while осуществляется итерирование по n. На каждой итерации используется ранее посчитанная сумма, к которой добавляются слагаемые, построенные для новых точек разбиения отрезка интегрирования.

Заметьте, подынтегральная функция f задана как поле класса.

Алгоритм 3.9 прост и эффективен для последовательного выполнения, но требует модификации для случая параллельного выполнения. К счастью, он допускает естественное распараллеливание. Более того, распараллеливание можно вести на двух уровнях. Во-первых, интервал интегрирования можно разбить на р отрезков и независимо вычислять интеграл на соответствующем отрезке. Далее останется только суммировать полученные значения. При этом на каждом отрезке можно использовать одну и ту же последовательную версию. Поскольку объем требуемой работы на каждом отрезке уменьшается, то при параллельном выполнении можно ожидать уменьшения общего объема работы.

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

Приведу метод, в котором распараллеливание ведется за счет разбиения интервала интегрирования на p отрезков:

/// <summary>
        /// Последовательный алгоритм
        /// Вычисление с разбиением интервала интегрирования
        /// на p сегментов
        /// </summary>
        public void SequenceIntegralWithSegments()
        {           
            double dx = (b - a) / p;           
            double start = 0, finish = 0;
            for (int i = 0; i < p; i++)
            {
                start = a + i * dx;
                finish = start + dx;
                DefiniteIntegral(start,
                   finish, out results[i]);
            }
            result = 0;
            for (int i = 0; i < p; i++)
            {
                result += results[i];
            }
        }

В цикле по числу отрезков вызывается последовательный алгоритм, вычисляющий значение интеграла на соответствующем отрезке. Этот цикл допускает распараллеливание. При наличии нескольких процессоров все они могут параллельно вычислять интеграл на своем отрезке интегрирования.

Заметьте, этот вариант алгоритма оказывается предпочтительнее и в случае проведения вычислений одним процессором. Причина в том, что на разных отрезках интегрирования функция может вести себя по-разному. Поскольку для каждого отрезка эффективно подбирается свое число N - число разбиений интервала интегрирования в методе прямоугольников, то общее время решения задачи может быть уменьшено.

Вычисление числа ПИ

Для параллельных вычислений классической тестовой задачей является задача вычисления числа ПИ. Обсудим и мы некоторые алгоритмы решения этой задачи, тем более что вся необходимая подготовка для ее решения уже выполнена.

Число ПИ и тригонометрические функции

Число ПИ естественным образом появляется при вычислении тригонометрических функций. Поскольку синус тридцати градусов равен одной второй, то вычислить ПИ можно по формуле:

$$\pi=6 \cdot ArcSin(0,5)$$

Ранее мы рассмотрели алгоритм вычисления функции ArcSin, способы распараллеливания этого алгоритма, и трудности, возникающие при распараллеливании.

Число ПИ и определенный интеграл

Вспомним, чему равна производная функции ArcTg(x)

$$\frac{d(arctg(x))}{dx}=\frac{1}{1+x^2}= > arctg(y) = \int_0^y \frac{1}{1+x^2} dx$$

Отсюда непосредственно следует:

$$\pi=4 \int_0^1\frac{1}{1+x^2} $$

Опять-таки, интеграл мы уже умеем считать, и знаем, как можно распараллелить алгоритм вычисления.

Существует и множество других разложений, позволяющих вычислить число ПИ с некоторой точностью. Одно из них в качестве примера использовалось нами при рассмотрении алгоритмов суммирования рядов (алгоритмы 3.5 и 3.6).

Алгоритмы линейной алгебры

Многие алгоритмы линейной алгебры легко распараллеливаются достаточно естественным образом. Все суперкомпьютеры при анализе их эффективности проходят тест Linpack, представляющий пакет различных алгоритмов линейной алгебры. Давайте рассмотрим некоторые классические алгоритмы линейной алгебры и обсудим возможности их распараллеливания.

Умножение матриц

Пусть матрица C размерности [m, n] является произведением двух прямоугольных матриц A и B, имеющих соответственно размерности [m, q] и [q, n]. Каждый элемент матрицы C является скалярным произведением двух векторов, заданных соответствующей строкой матрицы A и столбца матрицы B:

$$c_{i,j}=\sum_{k=1}^q a_{i,k} \cdot b_{k,j}, \text{ где } i=1..n; j=1..m$$

Очевидно, что временная сложность последовательного алгоритма равна:

$$T(Sec\_Algorithm) = O(m\cdot n\cdot q)$$

Что может дать распараллеливание? В идеальном случае, когда у нас в распоряжении суперкомпьютер или, еще лучше, метакомпьютер с неограниченным числом процессоров, все элементы матрицы C можно считать параллельно, что в m*n раз сокращает время вычислений. Если же для вычисления суммы (3.18) применять пирамидальный алгоритм, то сложность умножения матриц становится равной:

$$T(Ideal\_Parallel) = O(log_2 q)$$

Ускорение сверхзамечательное, но и его можно улучшить! Если использовать векторные процессоры, то вычисление скалярного произведения (3.18) становится элементарной операцией и тогда:

$$T(Super\_Ideal\_Parallel) = O(1)$$

Все прекрасно. Одна беда - эффективность крайне низкая, поскольку потребуется m*n*q/2 векторных процессоров. Конечно же, если задача важна, а именно такие задачи решаются на суперкомпьютерах, то ускорение важнее эффективности.

Если же в нашем распоряжении есть компьютер с конечным и относительно небольшим числом процессоров - p, существенно меньшим m, то время умножения матрицы можно сократить в p раз. Распараллеливание проще всего выполнить естественным образом, разбив матрицу A на полосы примерно одинаковой ширины. Если m не кратно p, то ширина полос может отличаться на единицу. Если например p равно 8, а n - 100, то четыре полосы могут иметь по 12 строк, а другие четыре полосы - по 13 строк. Каждый из 8 процессоров будет выполнять умножение своей полосы на матрицу B, получая соответствующую полосу матрицы C. В результате их совместной и параллельной работы задача умножения матриц будет решена. Алгоритм распараллеливания настолько прост, что я не буду приводить его код. По сути, он немногим отличается от последовательного алгоритма.

Алгоритм умножения слабо заполненных матриц

Давайте рассмотрим более интересную задачу умножения слабо заполненных матриц, у которых большинство элементов равны нулю. Это позволит нам рассмотреть алгоритмы со сложно организованными данными, использующими не только массивы, но и списки и структуры. Когда на практике приходится иметь дело с очень большими матрицами, то как правило они являются слабо заполненными. Иногда они могут иметь специальную блочную структуру с нулевыми блоками. Но будем рассматривать общий случай, не учитывающий возможную структуру строения матрицы.

Слабо заполненную матрицу будем хранить в виде массива из n элементов. В зависимости от того, как мы хотим хранить матрицу - по строкам или по столбцам, параметр n будет задавать число строк или число столбцов матрицы, и каждый элемент массива задает соответственно строку или столбец слабо заполненной матрицы. Каждый элемент этого массива представляет собой список структур. Каждая структура хранит два числа - значение ненулевого элемента и его индекс в строке или в столбце.

Если перемножаемые матрицы A и B заполнены на 10%, то такое представление матриц позволяет примерно в 5 раз сократить требуемую для их хранения память. Существенно сократится и время умножения матриц, поскольку все действия будут выполняться только над ненулевыми элементами.

Будем предполагать, что нули равномерно распределены в слабо заполненной матрице. Пусть p - вероятность того, что элемент такой матрицы не равен нулю, и эта вероятность мала. Что можно сказать о степени заполнения матрицы C = A * B? Нетрудно видеть, что для больших значений n произведение слабо заполненных матриц этим свойством обладать не будет. Вероятность того, что элемент матрицы C не равен нулю, можно посчитать по формуле:

$$p_c=1-(1-p^2)^n$$

При n = 100 и p = 0,1 вероятность $$p_с$$ = 0,634, но при n = 1000 она практически приближается к единице и равна 0,9999.

Ситуация меняется, если вероятность p является убывающей функцией от n. Пусть, например, p = 1/n, тогда

$$p_c=1-(1-\frac {1}{n^2})^n \approx \frac{1}{n}$$

Учитывая соотношения (3.20) и (3.21), можно сформулировать следующее утверждение:

Теорема:

Если при умножении слабо заполненных матриц A и B вероятность ненулевого элемента p постоянна, то с ростом размера матриц n вероятность заполнения результирующей матрицы C стремится к единице.

Если при умножении слабо заполненных матриц A и B вероятность ненулевого элемента p зависит от n и равна 1/n, то результирующая матрица C также является слабо заполненной с той же вероятностью заполнения, как и исходные матрицы $$p_c = p$$.

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

Начнем с первого случая и построим алгоритм для случая умножения слабо заполненных матриц, полагая, что их произведение является плотно заполненной матрицей. Первый сомножитель - матрицу A - будем хранить по строкам, а второй сомножитель - матрицу B - по столбцам. Каждую из этих матриц будем представлять в виде одномерного массива, элементы которого являются списками с ненулевыми элементами соответствующих строк и столбцов. Результирующая матрица C будет представлена традиционным способом в виде двумерного массива. Вот пример возможной программы на C#:

/// <summary>
        /// Умножение слабо заполненных матриц
        /// [n * m] * [m * q] = [n * q]
        /// </summary>
        /// <param name="A"></param>
        /// <param name="B"></param>
        /// <param name="C"></param>
        public void MultMatr(ArrayList[] A, ArrayList[] B, double[,] C)
        {           
            ArrayList listA, listB;
            Item itemA, itemB;
            int na = 0, nb = 0;
            int ia = 0, ib = 0;
            //Цикл по числу строк матрицы А
            for(int i =0; i < A.Length; i++)
            {
                listA = A[i];   //i-я строка А как список ненулевых элементов
                na = listA.Count;
                // Цикл по числу столбцов матрицы В
                for (int j = 0; j < B.Length; j++)
                {
                    listB = B[j];   //j-й столбец
                    nb = listB.Count;
                    ia = 0; ib = 0;
                    C[i, j] = 0;
                    while (ia < na  ib < nb)
                    {
                        itemA = (Item)listA[ia];
                        itemB = (Item)listB[ib];
                        if (itemA.index == itemB.index)
                        {
                            C[i, j] += itemA.val * itemB.val;
                            ia++;                           
                            ib++;                            
                        }
                        else
                            if (itemA.index < itemB.index)
                                ia++;                 
                            else
                                ib++; 
                    }
                }
            }
        }

Циклы for в этой программе допускают распараллеливание, позволяя все элементы результирующей матрицы C считать параллельно и одновременно. Понятно, что внутренний цикл while требует последовательного выполнения. Поэтому, если бы в нашем распоряжении были бы n*m свободных процессоров, то время вычисления матрицы C определялось бы временем работы внутреннего цикла и зависело бы как от размера матриц n, так и от вероятности заполнения, так что временную сложность можно было бы представить как O(p * n).

Следует отметить, что эта программа является оптимальной и для случая ее выполнения одним процессором.

Если ориентироваться на многоядерные компьютеры с фиксированным числом процессоров, то разумно рассмотреть алгоритм, предполагающий управляемое программистом распараллеливание. В этом случае естественный способ распараллеливания состоит в разбиении исходной матрицы A на полосы. Каждый процессор для умножения своей полосы может использовать выше приведенный алгоритм умножения прямоугольных матриц.

Вот пример возможной записи такой программы на C#:

/// <summary>
 /// Умножение слабо заполненных матриц
 /// [n * m] * [m * q ] = [n * q]
 /// Умножение идет p полосами 
 /// Предполагается, что каждую полосу умножает один процессор
 /// </summary>
 /// <param name="A"></param>
 /// <param name="B"></param>
 /// <param name="C"></param>
 /// <param name="p">число полос</param>
 public void MultMatrBar(ArrayList[] A, ArrayList[] B, double[,] C, int p)
        {          
            int n, m, q;
            n = A.Length;
            m = n / p;
            ArrayList[] barA = new ArrayList[m];
            q = B.Length;
            double[,] barC = new double[m, q];

            //Цикл по числу полос кроме последней
            for (int i = 0; i < p - 1; i++)
            {
                Array.Copy(A, i * m, barA, 0, m);
                MultMatr(barA, B, barC);
                Array.Copy(barC, 0, C, i * q, m * q); 
            }
            // ширина последней полосы
            // может отличаться 
            n = n - m * (p - 1);
            barA = new ArrayList[n];
            barC = new double[n, q];
            Array.Copy(A,(p - 1) * m, barA, 0, n);
            MultMatr(barA, B, barC);
            Array.Copy(barC, 0, C, (p - 1) * q, n * q); 
        }

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

Сортировка массивов

Рассмотрим возможности распараллеливания алгоритмов сортировки массивов.

Пузырьковая сортировка

Начнем с простейшего метода сортировки - пузырьковой сортировки. Идея последовательного алгоритма проста и элегантна. Массив можно рассматривать как некий вертикальный сосуд, заполненный элементами - пузырьками. Для массива из n элементов выполняется n-1 проход. Каждый проход начинается снизу - со дна сосуда, производя последовательный обмен соседних элементов, если нижний элемент "легче" верхнего. В результате прохода "легкие" элементы всплывают. На первом проходе самый легкий элемент окажется вверху сосуда. На i-м проходе на свое место всплывет i-й легкий элемент. Число сравнений на каждом проходе уменьшается на единицу. Временная сложность алгоритма - O(n2). Рис. 3.2 иллюстрирует алгоритм пузырьковой сортировки:

(рис 3.2) Алгоритм пузырьковой сортировки

Приведу текст записи классического варианта алгоритма на языке C#:

/// <summary>
        /// Классический вариант пузырьковой сортировки
        /// Со сложностью O(n * n)
        /// </summary>
        /// <param name="mas">сортируемый массив</param>
        public void BubbleSortClassic(double[] mas)
        {
            int n = mas.Length;            
            double temp = 0;
            for (int k = 0;  k < n - 1; k++)
            { //цикл по числу проходов               
                for (int i = n - 1; i > k; i--)
                {  // цикл всплытия
                    if (mas[i] < mas[i - 1])
                    {//swap
                        temp = mas[i];
                        mas[i] = mas[i - 1];
                        mas[i - 1] = temp;
                    }
                }
            }
        }

Если применить этот вариант алгоритма к уже отсортированному массиву, то он работает неэффективно, поскольку будет выполнять бесполезные проходы, не выполняя никаких обменов, поскольку всплывать некому. Первая возможная оптимизация состоит в том, чтобы на каждом проходе фиксировать существование обменов - факт всплытия элементов. Если на проходе не было ни одного обмена, то массив уже отсортирован и дальнейшие проходы бесполезны. Приведу текст такого варианта:

/// <summary>
        /// Пузырьковая сортировка
        /// Учитывает возможную отсортированность массива
        /// Проходы прекращаются, если отсутствуют обмены 
        /// на предыдущем проходе.
        /// Булевская переменная change следит за обменами
        /// </summary>
        /// <param name="mas">сортируемый массив</param>
        public void BubbleSort(double[] mas)
        {
            int n = mas.Length;
            bool change = true;
            double temp = 0;
            int k = 0;
            for (k = 0; change  k < n-1; k++)
            {
                change = false;
                for(int i = n-1; i > k; i--)
                {
                    if (mas[i] < mas[i - 1])
                    {//swap
                        temp = mas[i];
                        mas[i] = mas[i - 1];
                        mas[i - 1] = temp;
                        change = true;
                    }
                }
            }
        }

Усложнение алгоритма незначительное, но и эффект на хорошо перемешанных массивах минимален. Этот вариант стоит применять, когда есть предположения о частичной упорядоченности массива. Более интересна другая оптимизация. Если на очередном проходе запоминать индекс первого обмена, то на следующем проходе не требуется начинать проверку с самого нижнего элемента. Начальную точку определяет сохраненный индекс. Эта оптимизация хотя и не изменяет порядка сложности алгоритма, но зачастую работает быстрее, чем классический вариант. Вот код этого варианта алгоритма:

/// <summary>
        /// Вариант пузырьковой сортировки 
        /// с запоминанием индекса первого обмена
        /// Показывает лучшие результаты. 
        /// </summary>
        /// <param name="mas"></param>
        public void BubbleSortIndex(double[] mas)
        {
            int n = mas.Length;
            bool change = true;
            int index = n - 1;
            int i0 = 0;
            double temp = 0;
            for (int k = 0; change  k < n - 1; k++)
            {
                change = false;
                i0 = index < n - 1 ? index : n - 1;
                for (int i = i0; i > k; i--)
                {
                    if (mas[i] < mas[i - 1])
                    {//swap
                        temp = mas[i];
                        mas[i] = mas[i - 1];
                        mas[i - 1] = temp;
                        if(!change) index = i + 1;
                        change = true;                        
                    }
                }
            }

Во всех рассмотренных вариантах алгоритм последователен по своей сути, и никакой из циклов принципиально не распараллеливается, поскольку сравниваются два соседних элемента, и следующее сравнение зависит от результата предыдущего сравнения.

Нетрудно построить параллельную версию пузырьковой сортировки, ориентированную на выполнение сортировки p процессорами. Сортируемый массив можно разделить на p частей, каждая из которых независимо сортируется, используя последовательный вариант пузырьковой сортировки. Эта работа может выполняться параллельно. Затем происходит слияние отсортированных частей, что требует последовательного выполнения. Возникают два вопроса - как делить массив на части и как выполнять слияние отсортированных фрагментов?

При рассмотрении задачи суммирования элементов массива мы уже говорили, что разумными способами деления массива на p частей являются сегментное деление и шаговое. В данном случае целесообразно воспользоваться шаговым алгоритмом, поскольку он обеспечивает более быстрое всплывание легких элементов массива.

Слияние фрагментов массива можно выполнять, используя дополнительный массив, в который в нужном порядке будут сливаться элементы отсортированных фрагментов, как это делается в классической процедуре сортировки слиянием, когда сливаются два массива. В этом случае сложность слияния определяется как O(p * n). Нужно поставить на место n элементов, и на каждом шаге элемент нужно выбрать из p кандидатов. В целом сложность параллельного варианта пузырьковой сортировки задается соотношением $$O(n^2/p^2 +n\cdot p)$$.

Можно выполнять слияние на месте, но тогда сложность в худшем случае будет O(n), поскольку на каждом шаге придется восстанавливать отсортированность того фрагмента массива, который содержал минимальный элемент. После обмена место минимального элемента займет другой элемент, который может нарушить отсортированность фрагмента. Поскольку элемент, который попадает в верхушку отсортированного фрагмента, является минимальным элементом другого фрагмента, то для восстановления отсортированности обычно нужно выполнять небольшое число обменов. Учитывая, что вероятность худшего варианта мала, такая версия алгоритма может представлять практический интерес в ситуации, когда приходится экономить память.

Приведу текст реализации шагового варианта алгоритма, в котором слияние использует дополнительный массив:

/// <summary>
        /// Параллельный вариант пузырьковой сортировки         
        /// </summary>
        /// <param name="mas">сортируемый массив</param>
        /// <param name="p">число процессоров</param>
        public void BubbleParallel(double[] mas, int p)
        {
            int n = mas.Length;
            if (p > n) p = n;
            bool[] change = new bool[p];
            int[] index = new int[p];
            for (int i = 0; i < p; i++)
                index[i] = n - i - 1;
            int i0 = 0, m = 0;
            double temp = 0;
            //Цикл по числу процессоров
            for (int j = 0; j < p; j++)               
               {// процессор j сортирует подпоследовательность элементов, 
                //начинающихся индексом n -j -1 и отстоящих на расстоянии р
                   //цикл по числу проходов m
                   m = index[j] / p;
                   for (int k = 0; k < m; k++)
                   {
                       //цикл всплытия легкого элемента на k-м проходе
                       i0 = index[j] < n - j - 1 ? index[j] : n - j - 1;
                       change[j] = false;
                       for (int i = i0; i - p >= k * p; i = i - p)
                       {
                           if (mas[i] < mas[i - p])
                           {//swap
                               temp = mas[i];
                               mas[i] = mas[i - p];
                               mas[i - p] = temp;

                               if (!change[j]) index[j] = i + p;
                               change[j] = true;
                           }
                       }
                   }
               }
            //Слияние отсортированных последовательностей
            // Merge(mas, p);
            Merge1(mas, p);
        }

Процедура слияния p отсортированных фрагментов массива, где элементы каждого фрагмента отстоят на расстоянии p, имеет следующий вид:

/// <summary>
        /// Слияние упорядоченных последовательностей
        /// Последовательности представляют отрезки с шагом p
        /// Используется дополнительная память
        /// </summary>
        /// <param name="mas">сортируемый массив</param>
        /// <param name="p">число процессоров</param>
        public void Merge1(double[] mas, int p)
        {
            int n = mas.Length;
            int m = n / p;
            int index_min = 0;
            double min = 0;
            int i = 0;
            double[] tmas = new double[n];
            int[] start = new int[p], finish = new int[p];
            for (i = 1; i <= p; i++)
            {
                finish[p-i] = n - i;
                start[p - i] = finish[p - i] % p;               
            }
            for (int k = 0; k < n; k++)
            {//пересылка k-ого элемента
                //поиск кандидата
                i = 0;
                while (start[i] > finish[i]) i++;
                index_min = i; min = mas[start[i]];

                for (int j = i + 1; j < p; j++)
                {
                    //цикл по кандидатам
                    if (start[j] <= finish[j])
                    {
                        if (mas[start[j]] < min)
                        {
                            min = mas[start[j]];
                            index_min = j;
                        }
                    }
                }
                //pass                        
                tmas[k] = mas[start[index_min]];
                start[index_min] += p;
            }
            for (i = 0; i < n; i++)
                mas[i] = tmas[i];
        }

Эксперименты показывают, что этот вариант сортировки работает значительно эффективнее классического варианта пузырьковой сортировки уже на массивах длины 100. И это происходит даже в том случае, когда параллельные вычисления не выполняются и обработка идет последовательно. Объясняется это тем, что из-за пошагового характера алгоритма данный вариант сортировки приближается к версии сортировки Шелла.

/// <summary>
        /// Чет-нечет сортировка
        /// Вариация пузырьковой сортировки
        /// </summary>
        public void OddEvenSortSeq()
        {
            int n = mas.Length;
            int m = n/2;
            double temp = 0;
            for (int k = 0; k < n; k++)
            { //цикл по числу проходов               
                if (k % 2 == 0)                
                    for (int j = n - 1; j > 0; j -= 2)
                    {
                        if (mas[j] < mas[j - 1])
                        {
                            temp = mas[j];
                            mas[j] = mas[j - 1];
                            mas[j - 1] = temp;
                        }
                    }                
                else
                    for (int j = n - 2; j > 0; j -= 2)
                    {
                        if (mas[j] < mas[j - 1])
                        {
                            temp = mas[j];
                            mas[j] = mas[j - 1];
                            mas[j - 1] = temp;
                        }
                    }
            }
        }

/// <summary>
        /// Быстрая сортировка Хоара
        /// Вызов рекурсивной версии
        /// </summary>
        /// <param name="mas"></param>
        public void QuickSort(double[] mas)
        {
            QSort(mas, 0, mas.Length - 1);
        }

        /// <summary>
        /// Рекурсивная версия быстрой сортировки
        /// </summary>
        /// <param name="mas">сортируемый массив</param>
        /// <param name="start">индекс начала сортируемой части</param>
      /// <param name="finish">индекс конца сортируемой части массива </param>
        void QSort(double[] mas, int start, int finish)
        {
           
            if (finish - start > 0)
            {
                double cand = 0, temp = 0;
                int l = start, r = finish;
                cand = mas[(r + l) / 2];
                while (l <= r)
                {
                    while (mas[l] < cand) l++;
                    while (mas[r] > cand) r--;
                    if (l <= r)
                    {
                        temp = mas[l];
                        mas[l] = mas[r];
                        mas[r] = temp;
                        l++; r--;
                    }
                }
                QSort(mas, start, r);
                QSort(mas, l, finish);
            }
        }
        /// <summary>
        /// Быстрая сортировка Хоара
        /// Параллельная версия
        /// </summary>
        /// <param name="mas"></param>
        /// <param name="p">число процессоров</param>
        public void QuickSortParallel(double[] mas, int p)
        {
            int n = mas.Length;
            int m = n / p;
            for (int i = 0; i < p; i++)
            {
                int start = i * m;
                int finish = i != p - 1? start + m - 1 : n - 1;
                QSort(mas, start, finish);                
            }//Слияние             
          //MergeQ(mas, p);
            MergeQ1(mas, p);
        }

        /// <summary>
        /// Слияние упорядоченных последовательностей
        /// Последовательности представляют подряд идущие отрезки
        /// Используется дополнительная память
        /// </summary>
        /// <param name="mas">сортируемый массив</param>
        /// <param name="p">число процессоров</param>
        public void MergeQ1(double[] mas, int p)
        {
            int n = mas.Length;
            int m = n / p;
            int index_min = 0;            
            double min = 0;            
            int i = 0;
            double[] tmas = new double[n];
            int[] start = new int[p], finish = new int[p];
            for (i = 0; i < p; i++)
            {
                start[i] = i * m;
                finish[i] = i != p - 1 ? start[i] + m - 1 : n - 1;
            }
            for (int k = 0; k < n; k++)
            {//пересылка k-ого элемента
                //поиск кандидата
                i = 0;
                while (start[i] > finish[i]) i++;
                index_min = i; min = mas[start[i]]; 
                    
                for (int j = i+1; j < p; j++)
                {
                    //цикл по кандидатам
                    if (start[j] <= finish[j])
                    {
                        if (mas[start[j]] < min)
                        {
                            min = mas[start[j]];
                            index_min = j;                            
                        }
                    }
                }                    
                    //pass                        
                        tmas[k] = mas[start[index_min]];
                        start[index_min]++;
            }
            for (i = 0; i < n; i++)
                mas[i] = tmas[i];
        }
    }
}