Рассмотрим проблемы, возникающие при разработке высокоэффективных программ для современных вычислительных систем на примере элементарной задачи перемножения двух квадратных матриц. Поскольку эта задача чрезвычайно проста, то она особенно наглядно демонстрирует, что разработка высокоэффективных приложений представляет собой весьма нетривиальную задачу. С другой стороны, для этой задачи достаточно просто отслеживать производительность компьютера на реальном приложении. Для оценки этой величины мы будем измерять время выполнения фиксированного числа операций, которые необходимо выполнить для решения задачи. Для перемножения двух квадратных матриц размерности N нужно выполнить (2*N - 1)*N*N арифметических операций.

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

Конечно, это не совсем корректный тест, поскольку используются только операции умножения и сложения, но его результаты весьма поучительны. Будем выполнять перемножение достаточно больших матриц (1000x1000) и сравнивать полученную производительность с той, которую декларируют производители компьютеров. Тестирование будем выполнять на 3-х различных вычислительных системах: двухпроцессорной системе SUN Ultra60, Linux-кластере с узлами Pentium III 500 Mhz и двухпроцессорной системе Compaq Alpha DS20E.

Для компьютера Alpha производительность одного процессора оценивается в 1000 Mflops, для SUN Ultra60 - 800 Mflops, для Pentium III 500 Mhz - 500 Mflops. Будем называть эти величины пиковой или потенциальной производительностью и посмотрим, в какой мере мы можем достичь этих показателей на реальном приложении.

Для решения нашей задачи любой программист написал бы на языке Фортран77 что-то подобное приведенному ниже:

      program matmult
      integer i, j, k, n, nm
      parameter (n=1000)
      real*8 a(n,n), b(n,n), c(n,n), s
      real*8 time, dsecnd, op, mf
      op = (2.d0*n-1)*n*n
с  Блок начальной инициализации
      do 1 i = 1,n
      do 1 j = 1,n
      a(i,j) = dble(i)
      b(i,j) = 1./dble(j)
   1  continue
с  Вычислительный блок
      time = dsecnd()
      do 2 i = 1,n
      do 2 j = 1,n
      s =0.0D0
      do 3 k = 1,n
   3  s = s + a(i,k)*b(k,j)
      c(i,j) = s
   2  continue
с  Блок печати
      time = dsecnd() - time
      mf = op/(time*1000000.0)
      write(*,10) c(1,1),c(1,n),c(n,1),c(n,n)
  10  format(2x,2f16.6)
      write(*,*) ' time calculation: ', time,
     *'  mflops: ',mf
      end
с  Функция таймер
      double precision function dsecnd()
      real tarray(2)
      dsecnd = etime(tarray)
      return
      end

Откомпилируем программу в стандартном режиме:
f77 -o matmult matmult.f
и посмотрим результат на разных системах.

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

Откомпилируем программу в режиме оптимизации
на Ultra60 и Alpha: fast -O4
на Linux: -O.

Как и следовало ожидать, автоматическая оптимизация в различной степени ускорила решение задачи (в соответствии с обычной практикой в 2-3 раза), но все равно производительность удручающе мала. Причем наибольшую производительность демонстрирует как раз не самая быстрая машина (Pentium III под Linux).

Обратим внимание на то, что в операторе, помеченном меткой 3, выполняется выборка элементов строки из отстоящих далеко друг от друга элементов массива (массивы в Фортране размещаются по столбцам). Это не позволяет буферизовать их в быстрой кэш-памяти. Перепишем вычислительную часть программы следующим образом:

      do 2 i = 1,n
      do m = 1,n
      row(m) = a(i,m)
      end do
      do 2 j = 1,n
      c(i,j) =0.0d0
      do 3 k = 1,n
   3  c(i,j) = c(i,j) + row(k)*b(k,j)
   2  continue

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

Результаты заметно улучшились для компьютера Ultra60 и особенно для Alpha, которые обладают достаточно большой кэш-памятью. Это подтверждает наше предположение о важности эффективного использования кэш-памяти.

На эффективное использование кэш-памяти нацелены поставляемые с высокопроизводительными системами оптимизированные математические библиотеки BLAS. До недавнего времени это обстоятельство было одним из основных аргументов в конкурентной борьбе фирм-производителей высокопроизводительных систем с компьютерами "желтой" сборки. В настоящее время ситуация изменилась. Написана и бесплатно распространяется самонастраивающаяся библиотека BLAS (ATLAS), которая может устанавливаться на любом компьютере, под любой операционной системой.

В частности, в библиотеке BLAS есть готовая процедура перемножения матриц. Заменим в программе весь вычислительный блок на вызов соответствующей процедуры:
CALL DGEMM('N', 'N', N, N, N, ONE, A, N, B, N, ZERO, C, N)

Приведем (табл. 3.5) результаты для программы, в которой для перемножения матриц используется подпрограмма DGEMM из библиотеки ATLAS. Заметим, что эффективность этой подпрограммы не уступает, а в некоторых случаях и превосходит эффективность подпрограммы из стандартных библиотек, поставляемых с программным обеспечением компьютеров (Sun Performance Library для Ultra60 и Common Extended Math Library для Alpha).

Итак, в конце концов, мы получили результаты для производительности, близкие к тем, которые декларируются производителями вычислительных систем. По сравнению с первым вариантом расчета (табл. 3.2) ускорение работы программы для систем Ultra60, Pentium III и Alpha составило соответственно 96, 31 и 48 раз!

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

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

Максимального ускорения работы программы можно добиться при использовании параллельной версии подпрограммы перемножения матриц из пакета ScaLAPACK совместно с библиотекой ATLAS. Приведем в таблице 3.7 данные по производительности при выполнении этой программы на 2-х процессорах.

Итак, в окончательном варианте мы смогли решить задачу на компьютере Ultra60 за 1.62 сек., или в 161 раз быстрее по сравнению с первоначальным расчетом, а на компьютере Alpha за 1.3 сек., или в 83 раза быстрее. Несколько замечаний следует сказать по поводу Linux-кластера. Двухпроцессорный вариант на нем оказался менее эффективным, чем на других системах, и, по-видимому, нет смысла решать эту задачу на большем числе процессоров. Это связано с медленной коммуникационной средой. Доля накладных расходов на пересылки данных слишком велика по сравнению с вычислительными затратами. Однако, как мы отмечали ранее, эта доля уменьшается при увеличении размеров матрицы. В самом деле, при увеличении параметра N до 4000 с приемлемой эффективностью задача решается даже на 4-х процессорах. Суммарная производительность на 4-х процессорах составила 1231 Mflops, или по 307 Mflops на процессор, что, конечно, весьма неплохо.

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