Sparse BLAS

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

Собственно функции по большей части будут те же самые, но к данным, которые им передаются, несколько другая адресация. В частности, поскольку у нас есть не все элементы (каких-то элементов может не быть), то у нас индексы k₁, k₂ ,…, k_z не обязательно лежат последовательно (т.е. 1, 2, 3, 4… ). Могут лежать в любом другом порядке — 1, 10, 15 и т.д. Соответственно, когда мы заступаемся к элементам нашего массива, грубо говоря ужатого, то нам надо хранить не только элементы этого массива, но ещё и индексы, по которым эти элементы лежат.

(Вопрос из аудитории: "Могут ли элементы лежат не в порядке возрастания?" Лектор отвечает, что это ничему не противоречит).

Все функции, которые были в BLAS, остаются теми же самыми, но к некоторым из них можно добавить индекс "i", что означает "indexed", то есть передаётся ещё один параметр, который хранит индексы элементов. Не только сами элементы, но ещё и индексы.

Так, например есть функция ?axpy(), которую вы использовали в прошлый раз, в которую добавляется индекс "i" (?axpyi()), функция ?dot() – скалярное произведение, в которую также добавляется индекс - ?doti() и т.д.

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

Существуют несколько способов именования функции. Есть стандартный – "mkl_", есть упрощённый интерфейс (для Фортрана – все индексы начинаются с 1, в С – все индексы начинаются с 0).

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

Операции те же самые: умножение матрицы на вектор, матрицы на матрицу, решение системы линейных уравнений с одной правой частью, с множеством правых частей.

Теперь о различных способах хранения матриц. Если у нас есть матрица (как на слайде), то один из форматов называется CSR (Compressed Space Row), то есть "сжатый по рядам". Хранится с помощью четырёх массивов. В первом массиве (values) подписываются все элементы, которые для нас имеют значение, во втором массиве (columns) подписываются колонки, в которых лежат эти значения, в массиве PointerB записываются элементы в следующем формате – если мы из какого-то j-ого элемента этого массива вычтем его первый элемент и добавим единицу, то получим индекс элемента из массива values, который указывает на начало строки. Массив PointerB показывает аналогично на последний элемент строки. Для j-ой строки мы берём j-ый элемент массива PointerB, вычитаем от него первый элемет и прибавляем единицу. Получаем индекс, который в массиве values указывает на последний элемент в j-ой строке.

Интересен и формат Представления матриц, когда записывается для определённой матрицы значения элементов в массиве values, в массиве rows указывается, в каком ряду лежат соответствующие элементы и в массиве columns указывается, в каких колонках лежат эти элементы.

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

Диагональный формат. На слайде мы видим матрицу с диагональным преобладанием. С помощью диагонального формата мы можем представить в виде двух массивов. Distance задаёт расстояние между диагональю и главной диагональю. Притом положительное значение означает, что диагональ лежит выше основной диагонали, отрицательное – ниже основной диагонали. В массиве values записываются значения, которые лежат на соответствующих диагоналях. Звёздочками обозначаются те элементы, которые добавлены в диагональ, так как в матрице все диагонали разной длины, кроме основной, поэтому добавляем элементы. То есть, для примера, для диагонали, лежащей на три элемента ниже основной диагонали, мы ставим в distance "-3". Звездочками обозначаются незначимые элементы – можно прописывать нулями. Диагонали записываются в качестве столбцов массива distance.

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

Block Compressed Space Row – вариант блочного распределения матрицы. То есть у нас если матрица состоит из нескольких прямоугольных блоков, то мы можем сжимать отдельно эти блоки и откладывать в массивы.

Также используется Skyline Storage Format – контурное представление матриц.