Итеративный параллелизм: умножение матриц
Итеративная последовательная программа использует для обработки данных и вычисления результатов циклы типа for и while. Итеративная параллельная программа содержит несколько итеративных процессов. Каждый процесс вычисляет результаты для подмножества данных, а затем эти результаты собираются вместе.
В качестве простого примера рассмотрим задачу из области научных вычислений. Предположим, даны матрицы а и Ь, у каждой по п строк и столбцов, и обе инициализированы. Цель — вычислить произведение матриц, поместив результат в матрицу с размером пхп. Для этого нужно вычислить п2 промежуточных произведений, по одному для каждой пары строк и столбцов.
Матрицы являются разделяемыми переменными, объявленными следующим образом. double a[n,n], b[n,n], c[n,n];
При условии, что п уже объявлено и инициализировано, это объявление резервирует память для трех массивов действительных чисел двойной точности. По умолчанию индексы строк и столбцов изменяются от 0 до п-1.
После инициализации массивов а и Ь можно вычислить произведение матриц по такой последовательной программе.
for [i = 0 to n-1] { for [j = 0 to n-1] {
28 Глава 1. Обзор области параллельных вычислений
# вычислить произведение а[i,*] и b[*,j]
c[i,j] = 0.0;
for [k = 0 to n-1]
c[i,j] = c[i,j] + a[i,k]*b[k,j]; } }
Внешние циклы (с индексами i и j) повторяются для каждой строки и столбца. Во внутреннем цикле (с индексом k) вычисляется промежуточное произведение строки i матрицы а и столбца j матрицы Ь; результат сохраняется в ячейке с [ i, j ]. Строка с символом # в начале является комментарием.
Умножение матриц — это пример приложения с массовым параллелизмом, поскольку программа содержит большое число операций, которые могут выполняться параллельно. Две операции могут выполняться параллельно, если они независимы. Предположим, что множество чтения операции содержит переменные, которые она читает, но не изменяет, а множество записи — переменные, которые она изменяет (и, возможно, читает).
Две операции являются независимыми, если их множества записи не пересекаются. Говоря неформально, процессы всегда могут безопасно читать переменные, которые не изменяются. Однако двум процессам в общем случае небезопасно выполнять запись в одну и ту же переменную или одному процессу читать переменную, которая записывается другим. (Эта тема рассматривается подробно в главе 2.)
При умножении матриц вычисления промежуточных произведений являются независимыми операциями. В частности, строки с 4 по 6 приведенной выше программы выполняют инициализацию и последующее вычисление элемента матрицы с. Внутренний цикл программы читает строку матрицы а и столбец матрицы Ь, а затем читает и записывает один элемент матрицы с. Следовательно, множество чтения для внутреннего произведения — это строка матрицы а и столбец матрицы Ь, а множество записи — элемент матрицы с.
Поскольку множества записи внутренних произведений не пересекаются, их можно выполнять параллельно. Возможны варианты, когда параллельно вычисляются результирующие строки, результирующие столбцы или группы строк и столбцов. Ниже будет показано, как запрограммировать такие параллельные вычисления.
Сначала рассмотрим параллельное вычисление строк матрицы с. Его можно запрограммировать с помощью оператора со (от "concurrent" — "параллельный"):
со [i=0 to n-1] { # параллельное вычисление строк for [j = 0 to n-1] { c[i,j] = 0.0; for [k = 0 to n-1]
c[i,j] = c[i,j] + a[i,k]*b[k,j]; } }
Между этой программой и ее последовательным вариантом есть лишь одно синтаксическое различие — во внешнем цикле оператор for заменен оператором со. Но семантическая разница велика: оператор со определяет, что его тело для каждого значения индекса i будет выполняться параллельно (если не в действительности, то, по крайней мере, теоретически, что зависит от числа процессоров).
Другой способ выполнения параллельного умножения матриц состоит в параллельном вычислении столбцов матрицы с.
Его можно запрограммировать следующим образом.
со [j = 0 to n-1] { #параллельное вычисление столбцов for [i = 0 to n-1] { c[i,j] = 0.0; for [k = 0 to n-1]
c[i,j] = c[i,j] + a[i,k]*b[k,j]; } }
1.4. Итеративный параллелизм: умножение матриц 29
Здесь два внешних цикла (по i и по j) поменялись местами. Если тела двух циклов независимы и приводят к вычислению одинаковых результатов, их можно безопасно менять местами, как это было сделано здесь. (Это и другие аналогичные преобразования программ рассматриваются в главе 12.)
Программу для параллельного вычисления всех промежуточных произведений можно составить несколькими способами. Используем один оператор со для двух индексов.
со [i = 0 to n-1, j = 0 to n-1] { # все строки и с[i,j] = 0.О; #все столбцы
for [k = 0 to n-1]
c[i,j] = c[i,j] + a[i,k]*b[k,j]; }
}
Тело оператора со выполняется параллельно для каждой комбинации значений индексов i и j, поэтому программа задает п2 процессов. (Будут ли они выполняться параллельно, зависит от конкретной реализации.) Другой способ параллельного вычисления промежуточных произведений состоит в использовании вложенных операторов со.
со [i = 0 to n-1] { # строки параллельно, затем со [j = 0 to n-1] { # столбцы параллельно c[i,j] = 0.0; for [k = 0 to n-1]
c[i,j] = c[i,j] + a[i,k]*b[k,j]; } }
Здесь для каждой строки (внешний оператор со) и затем для каждого столбца (внутренний оператор со) задается по одному процессу. Третий способ написать эту программу — поменять местами две строки последней программы. Результат всех трех программ будет одинаковым: выполнение внутреннего цикла для всех п2 комбинаций значений i и j. Разница между ними — в задании процессов, а значит, и во времени их создания.
Заметим, что все параллельные программы, приведенные выше, были получены заменой оператора for на со.
Но мы сделали это только для индексов i и j. А как быть с внутренним циклом по индексу k? Нельзя ли и этот оператор заменить оператором со? Ответ — "нет", поскольку тело внутреннего цикла как читает, так и записывает переменную с [ i, j ]. Промежуточное произведение — цикл for с переменной k — можно вычислить, используя двоичный параллелизм, но для большинства машин это непрактично (см. упражнения в конце главы).
Другой способ определить параллельные вычисления, показанные выше, — использовать декларацию (объявление) process вместо оператора со. В сущности, process — это оператор со, выполняемый как "фоновый". Например, первая параллельная программа из показанных выше — та, что параллельно вычисляет строки результата, — может быть записана следующим образом. process row[i = 0 to n-1] { # строки параллельно for [j = 0 to n-1] { c[i,j] = 0.0; for [k = 0 to n-1]
c[i,j] = c[i,j] + a[i,k]*b[k,j]; } }
Здесь определен массив процессов — row [ 1 ], row [ 2 ] и т.д. — по одному для каждого значения индекса i. Эти п процессов создаются и начинают выполняться, когда встречается данная строка описания. Если за декларацией процесса следуют операторы, то они выполняются параллельно с процессом, тогда как операторы, записанные после оператора со, не выполняются до его завершения. Декларации процесса, в отличие от операторов со, не могут быть вложены в другие декларации или операторы. Декларации процессов и операторы со подробно описаны в разделе 1.9.
30 Глава 1 Обзор области параллельных вычислений
В программах, приведенных выше, для каждого элемента, строки или столбца результирующей матрицы использовано по одному процессу. Предположим, что число процессоров в системе меньше п (так обычно и бывает, особенно когда п велико). Остается еще очевидный способ полного использования всех процессоров: разделить матрицу на полосы (строк или столбцов) и для каждой полосы создать рабочий процесс.
В частности, каждый рабочий процесс вычисляет результаты для элементов своей полосы. Предположим, что есть Р процессоров и п кратно р (т.е. п делится на Р без остатка). Тогда при использовании полос строк оабочие пооиессы можно запоогоаммивовать слелуюшим обоазом.
Отличие этой программы от предыдущей состоит в том, что п строк делятся на Р полос, по п/Р строк каждая. Для этого в программу добавлены операторы, необходимые для определения первой и последней строки каждой полосы. Затем строки полосы указываются в цикле (по индексу i), чтобы вычислить промежуточные произведения для этих строк.
Итак, существенным условием распараллеливания программы является наличие независимых вычислений, т.е. вычислений с непересекающимися множествами записи. Для произведения матриц независимыми вычислениями являются промежуточные произведения, поскольку каждое из них записывает (и читает) свой элемент с [ i, j ] результирующей матрицы. Поэтому можно параллельно вычислять все промежуточные произведения, строки, столбцы или полосы строк. И, наконец, параллельные программы можно записывать, используя операторы со или объявления process.
1.5. Рекурсивный параллелизм: адаптивная квадратура
Программа считается рекурсивной, если она содержит процедуры, которые вызывают сами себя — прямо или косвенно. Рекурсия дуальна итерации в том смысле, что рекурсивные программы можно преобразовать в итеративные и наоборот. Однако каждый стиль программирования имеет свое применение, поскольку одни задачи по своей природе рекурсивны, а другие — итеративны.
В теле многих рекурсивных процедур обращение к самой себе встречается больше одного раза. Например, алгоритм quicksort часто используется для сортировки. Он разбивает массив на две части, а затем дважды вызывает себя: первый раз для сортировки левой части, а второй — для правой. Многие алгоритмы для работы с деревьями и графами имеют подобную структуру.
Рекурсивную программу можно реализовать с помощью параллелизма, если она содержит несколько независимых рекурсивных вызовов.
Два вызова процедуры (или функции) явля ются независимыми, если их множества записи не пересекаются. Это условие выполняется, если: 1) процедура не обращается к глобальным переменным или только читает их; 2) аргументы и результирующие переменные (если есть) являются различными переменными. Например, если процедура не обращается к глобальным переменным и имеет только параметры-значения, то любой ее вызов будет независимым. (Хорошо, если процедура читает и записывает только локальные переменные, тогда каждый экземпляр процедур имеет локальную копию переменных.) В соответствии с этими требованиями можно запрограммировать и алгоритм быстрой сортировки. Рассмотрим еще один интересный пример.
Каждая итерация вычисляет площадь малой фигуры по правилу трапеций и добавляет ее к общему значению площади. Переменная width — ширина каждой трапеции. Отрезки перебираются слева направо, поэтому правое значение каждой итерации становится левым значением следующей итерации.
Второй способ аппроксимации интеграла— использовать парадигму "разделяй и властвуй" и переменное число интервалов. В частности, сначала вычисляют значение m — середину отрезка между а и Ь. Затем аппроксимируют площадь трех областей под кривой, определенной функцией f (): от а до т, от m до b и от а до Ь. Если сумма меньших площадей равна большей площади с некоторой заданной точностью EPSILON, то аппроксимацию можно считать достаточной. Если нет, то большая задача — от а до Ь — делится на две подзадачи — от а до m и от m до Ь, и процесс повторяется. Этот способ называется адаптивной квадратурой, поскольку алгоритм адаптируется к форме кривой. Его можно запрограммировать так.
Интеграл функции f (х) от а до Ь аппроксимируется таким вызовом функции:
area = quad(a, b, f(a), f(b), (f(a)+f(b))*(b-a)/2);
В функции снова используется правило трапеции. Значения функции f () в крайних точках отрезка и приближенная площадь этого интервала передаются в каждый вызов функции quad, чтобы не вычислять их больше одного раза.
Итеративную программу нельзя распараллелить, поскольку тело цикла и считывает, и записывает значение переменной area. Тем не менее в рекурсивной программе вызовы функции quad независимы при условии, что вычисление функции f (х) не дает побочных эффектов. В частности, аргументы функции quad передаются по значению, и в ее теле нет присваи-
32 Глава 1. Обзор области параллельных вычислений
вания глобальным переменным. Таким образом, для задания параллельного выполнения рекурсивных вызовов функции можно использовать оператор со.
со larea = quad(left, mid, fleft, fmid, larea); // rarea = quad(mid, right, fmid, fright, rarea); oc
Это единственное изменение, необходимое для того, чтобы сделать данную программу рекурсивной. Поскольку оператор со не заканчивается до тех пор, пока не будут завершены оба вызова функций, значения переменных larea и rarea вычисляются до того, как функция quad возвратит их сумму.
В операторах со программ умножения матриц содержатся списки инструкций, выполняемых для каждого значения счетчиков (i и j). В операторе со, приведенном выше, содержатся два вызова функций, разделенных знаками //. Первая форма оператора со используется для выражения итеративного параллелизма, вторая — рекурсивного.
Итак, программу с несколькими рекурсивными вызовами функций можно легко преобразовать в параллельную рекурсивную программу, если вызовы независимы. Однако существует чисто практическая проблема: параллельно выполняемых операций может стать слишком много. Каждый оператор со в приведенной выше программе создает два процесса, по одному для каждого вызова функции. Если глубина рекурсии велика, то появится большое число процессов, возможно, слишком большое для параллельного выполнения. Решение этой проблемы состоит в сокращении, или отсечении, дерева рекурсии при достаточном количестве процессов, т.е. переключении'с параллельных рекурсивных вызовов на последовательные.Эта тема рассматривается в упражнениях и далее в этой книге.
