Перемножение плотных неособенных матриц (последовательный вещественный вариант)

Перемножение плотных неособенных матриц (последовательный вещественный вариант)

Перемножение матриц - одна из базовых задач в алгоритмах линейной алгебры, широко применяется в большом количестве разных методов. Здесь мы рассмотрим умножение [math]C = AB[/math]   плотных неособенных матриц (последовательный вещественный вариант), то есть тот вариант, где никак не используются ни специальный вид матрицы, ни ассоциативные свойства операции сложения [1] .

1.2 Математическое описание алгоритма

Исходные данные: плотная матрица [math]A[/math] (элементы [math]a_[/math] ), плотная матрица [math]B[/math] (элементы [math]b_[/math] ).

Вычисляемые данные: плотная матрица [math]C[/math] (элементы [math]c_[/math] ).

[math] \begin c_ = \sum_^ a_ b_, \quad i \in [1, m], \quad j \in [1, l]. \end [/math]

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

1.3 Вычислительное ядро алгоритма

Вычислительное ядро перемножения плотных неособенных матриц можно составить из множественных (всего их [math]l[/math] ) вычислений умножения матрицы [math]A[/math] на столбцы матрицы [math]B[/math] , или (при более детальном рассмотрении), из множественных (всего их [math]ml[/math] ) скалярных произведений строк матрицы [math]A[/math] на столбцы матрицы [math]B[/math] :

в режиме накопления или без него, в зависимости от требований задачи.

1.4 Макроструктура алгоритма

Как уже записано в описании ядра алгоритма, основную часть умножения матриц составляют множественные (всего [math]ml[/math] ) вычисления скалярных произведений строк матрицы [math]A[/math] на столбцы матрицы [math]B[/math]

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

1.5 Схема реализации последовательного алгоритма

Для всех [math]i[/math] от [math]1[/math] до [math]m[/math] и для всех [math]j[/math] от [math]1[/math] до [math]l[/math] выполняются

Особо отметим, что вычисления сумм вида [math]\sum_^ a_ b_[/math] производят в режиме накопления прибавлением к текущему (временному) значению вычисляемого элемента матрицы [math]c_[/math] произведений [math]a_ b_[/math] для [math]k[/math] от [math]1[/math] до [math]n[/math] , c возрастанием [math]k[/math] , вначале все элементы инициализируются нулями. При суммировании "по убыванию" общая схема принципиально не отличается и потому нами не рассматривается. Другие порядки выполнения суммирования приводят к изменению параллельных свойств алгоритма и будут рассматриваться нами в отдельных описаниях.

1.6 Последовательная сложность алгоритма

Для умножения двух квадратных матриц порядка [math]n[/math] (т.е. при [math]m=n=l[/math] ) в последовательном (наиболее быстром) варианте требуется:

  • по [math]n^3[/math] умножений и сложений.

Для умножения матрицы размером [math]m[/math] строк на [math]n[/math] столбцов на матрицу размером [math]m[/math] строк на [math]n[/math] столбцов в последовательном (наиболее быстром) варианте требуется:

  • по [math]mnl[/math] умножений и сложений.

При этом использование режима накопления требует совершения умножений и сложений в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает затраты во времени, требуемом для выполнения умножения матриц.

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

1.7 Информационный граф

Опишем граф алгоритма как аналитически, так и в виде рисунка.

Граф алгоритма умножения плотных матриц состоит из одной группы вершин, расположенной в целочисленных узлах трёхмерной области, соответствующая ей операция [math]a+bc[/math] .

Естественно введённые координаты области таковы:

  • [math]i[/math] — меняется в диапазоне от [math]1[/math] до [math]m[/math] , принимая все целочисленные значения;
  • [math]j[/math] — меняется в диапазоне от [math]1[/math] до [math]l[/math] , принимая все целочисленные значения;
  • [math]k[/math] — меняется в диапазоне от [math]1[/math] до [math]n[/math] , принимая все целочисленные значения.

Аргументы операции следующие:

  • [math]a[/math] :
    • при [math]k = 1[/math] константа [math]0[/math] ;
    • при [math]k \gt 1[/math] — результат срабатывания операции, соответствующей вершине с координатами [math]i, j, k-1[/math] ;

    Результат срабатывания операции является:

    • при [math]k \lt n[/math] - промежуточным данным алгоритма;
    • при [math]k = n[/math] - выходным данным [math]c_[/math] .

    Интерактивное изображение графа алгоритма без входных и выходных данных для случая перемножения двух квадратных матриц порядка 3 и 4

    1.8 Ресурс параллелизма алгоритма

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

    • по [math]n[/math] ярусов умножений и сложений (в каждом из ярусов — [math]n^2[/math] операций).

    Для умножения матрицы размером [math]m[/math] строк на [math]n[/math] столбцов на матрицу размером [math]n[/math] строк на [math]l[/math] столбцов в последовательном (наиболее быстром) варианте требуется:

    • по [math]n[/math] ярусов умножений и сложений (в каждом из ярусов — [math]ml[/math] операций).

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

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

    1.9 Входные и выходные данные алгоритма

    Входные данные: матрица [math]A[/math] (элементы [math]a_[/math] ), матрица [math]B[/math] (элементы [math]b_[/math] )).

    Объём входных данных: [math]mn+nl[/math]

    Выходные данные: матрица [math]C[/math] (элементы [math]c_[/math] ).

    Объём выходных данных: [math]ml[/math]

    1.10 Свойства алгоритма

    Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является квадратичным или билинейным (отношение кубической или трилинейной к линейной).

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

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

    2 Программная реализация алгоритма

    2.1 Особенности реализации последовательного алгоритма

    В простейшем варианте алгоритм умножения матриц на Фортране можно записать так:

    При этом для реализации режима накопления переменная [math]S[/math] должна быть двойной точности.

    2.2 Локальность данных и вычислений 2.2.1 Локальность реализации алгоритма 2.2.1.1 Структура обращений в память и качественная оценка локальности

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

    Исходя из исходного кода, можно увидеть, что всего встречается 6 разных видов фрагментов для отдельных массивов (эти виды выделены на рис.2 зеленым). Стоит сделать два уточнения: 1) профиль для результирующего массива С всегда в 2 раза больше, поскольку к элементам этого массива всегда происходят по два обращения подряд; 2) если во внутреннем цикле используется элементы массива, то обращение к этому элементу выносится за пределы внутреннего цикла (в цикле используется скалярная переменная вместо него).

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

    Рассмотрим подробнее каждый из данных 6 фрагментов на примере массива С.

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

    Если рассмотреть фрагмент более подробно (рис.4), можно увидеть, что каждый перебор является на самом деле последовательным перебором, при этом к каждому элементу происходит обращение дважды.

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