Итерационные методы для разреженных систем линейных уравнений: Схемы хранения, Базовые операции на разреженных матрицах
Автор:Yousef Saad
Перевод: Дяченко Т.Ф.
Источник:http://www-users.cs.umn.edu/~saad/PS/iter1.pdf
Схемы хранения разреженных матриц
Основной целью является хранение только ненулевых элементов матрицей и возможность производить над ней простейшие матричные операции.
Пусть Nz — количество ненулевых элементов матрицы A. Рассмотрим наиболее распространенные схемы хранения, более подробная информация приводится в [4].
Самой простой схемой хранения разреженных матриц является так называемый «координатный» формат. Структура исходной матрицы отображается на три одномерных массива: (1) — массив, содержащий значения ненулевых элементов матрицы А в произвольном порядке; (2) целочисленный массив, содержащий их строковые индексы; и (3) целочисленный массив, содержащий их столбцовые индексы. Все три вектора имеют длину Nz.
Пример 1. Матрица
Может быть представлена как:
| AA | 12 | 9 | 7 | 5 | 1 | 2 | 11 | 3 | 6 | 4 | 8 | 10 |
| JR | 5 | 3 | 3 | 2 | 1 | 1 | 4 | 2 | 3 | 2 | 3 | 4 |
| JC | 5 | 5 | 3 | 4 | 1 | 4 | 4 | 1 | 1 | 2 | 4 | 3 |
В приведенном примере элементы перечислены в произвольном порядке, однако обычно они перечисляются по строкам или столбцам. Если элементы перечисляются по строкам, то вектор JC, который содержит избыточную информацию, может быть заменен массивом указателей на начало каждой строки. Это даст заметную экономию при хранении. Новая структура данных будет иметь три массива, которые будут выполнять следующие функции:
- массив АА содержит ненулевые значения aij, сохраняемые построчно, от строки 1 до n. Длина массива АА есть Nz.
- целочисленный массив JA содержит столбцовые индексы элементов aij в матрице А. Длина массива JА есть Nz.
- целочисленный массив содержит указатели входа для каждой строки в массивах АА и JА. Таким образом, значение IA(i) показывается позиция, где начинается i-тая строка в массивах АА и JА. Длина IA есть n+1, где IA(n+1) содержит число IA(1)+ Nz, то есть адрес в А и JА начала фиктивной строки с номером n+1.
Таким образом, представленная выше матрица А может быть сохранена как:
| AA | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
| JА | 1 | 4 | 1 | 2 | 4 | 1 | 3 | 4 | 5 | 3 | 4 | 5 |
| IA | 1 | 3 | 6 | 10 | 12 | 13 |
Такой формат является наиболее популярным для хранения разреженных матриц общего вида. Он носит название разреженного строчного формата или Compressed Sparse Row (CSR). Эта схема более предпочтительна, чем координатная, так как часто является более удобной для нескольких важных опе¬раций над разреженными матрицами: сложения, умножения, перестановок строк и столбцов, транспонирования, решения линейных систем с разреженными матрицами коэффициентов как прямыми, так и итерационными методами. С другой стороны, координатная схема более проста, наглядна и гибка, она часто используется в качестве «входного» формата для вычислительных библиотек, предназначенных для работы с разреженными матрицами.
Существует несколько модификаций данной схемы хранения, наиболее очевидной является разреженная столбцовая схема или Compressed Sparse Column (CSC).
Другая распространенная схема учитывает тот факт, что диагональные элементы большинства матриц ненулевые и/или используются чаще других. Модифицированная строчная схема (Modified Sparse Row, MSR) использует только два массива: массив значении АА и целочисленный массив JА. В первых n позициях АА содержит диагональные элементы исходной матрицы по порядку. Элемент АА(n+1) не заполняется (или же несет дополнительную информацию о матрице). Начиная с позиции n+2 записываются ненулевые элементы исходной матрицы по строкам, исключая диагональные. Для каждого элемента АА(k) элемент JA(k) показывает столбцовый индекс в исходной матрице. На n+1 позициях матрицы JA размещаются указатели входа для каждой строки матрицы АА и JA.
Поэтому для матрицы из примера выше эти два массива будут иметь вид:
| AA | 1 | 4 | 7 | 11 | 12 | * | 2 | 3 | 5 | 6 | 8 | 9 | 10 |
| JА | 7 | 8 | 10 | 13 | 14 | 14 | 4 | 1 | 4 | 1 | 4 | 5 | 3 |
Звездочка указывает на неиспользуемый элемент. Заметим, что JA(n)=JA(n+1)=14, показывая, что последняя строка является нулевой, так как диагональный элемент уже описан ранее.
Диагонально структурированные матрицы — это матрицы, чьи ненулевые элементы расположены вдоль небольшого числа диагоналей. Эти диагонали могут храниться в прямоугольном массиве DIAG(1:n, 1:Nd), где Nd — число диагоналей. Смещение каждой диагонали вычисляется относительно главной диагонали, которая должна быть известна. Смещения хранятся в массиве IOFF(1:Nd). Таким образом, элемент ai,i+ioff(j) исходной матрицы будет храниться в позиции (i,j) в массиве DIAG. Порядок, в каком хранятся диагонали, в общем случае, неважен, особенно если большинство операций сосредоточенно около главной диагонали, ее имеет смысл хранить в первом столбце. Также следует отметить, что все диагонали, кроме главной, имеют менее n элементов, так что не все элементы матрицы DIAG будут использоваться.
Пример 2. Трехдиагональная матрица
DIAG =
| * | 1 | 2 |
| 3 | 4 | 5 |
| 6 | 7 | 8 |
| 9 | 10 | * |
| 11 | 12 | * |
| -1 | 0 | 2 |
Более общей схемой, популярной для векторных машин является так называемый формат. Эта схема предполагает, что количество диагоналей Nd невелико. Для хранения требуются два прямоугольных массива размером n*Nd каждый.
В первом, COEF, подобном DIAG, хранятся ненулевые элементы матрицы А. Ненулевые элементы каждой строки могут сохраняться в строке массива COEF(1:n,1:Nd), дополняя эту строку нулями, если необходимо. Вместе с COEF, целочисленный массив JCOEF(1:n,1:Nd) должен хранить индекс столбца каждого элемента COEF.
Для матрицы А из примера 2, схема хранения Ellpack-Itpack примет вид:
COEF =
| 1 | 2 | 0 |
| 3 | 4 | 5 |
| 6 | 7 | 8 |
| 9 | 10 | 0 |
| 11 | 12 | 0 |
| 1 | 3 | 1 |
| 1 | 2 | 4 |
| 2 | 3 | 5 |
| 3 | 4 | 4 |
| 4 | 5 | 5 |
Базовые операции для разреженных матриц
Умножение матрицы на вектор - важная операция, которая требуется в большинстве итерационных алгоритмов решения разреженных линейных систем. В этом разделе показано, как она может быть реализована для схем хранения, приведенных ранее.
В следующем примере кода на Fortran 90 представлен главный цикл умножения матрицы на вектор для матриц хранящихся в разреженном столбцовом формате:
DO 1=1, N Kl = IA(I) K2 = IA(I+1)-1 Y(I) = D0TPR0DUCT(A(K1:K2),X(JA(K1:K2))) ENDDO
Обратите внимание, что каждая итерация цикла вычисляет различные компоненты результирующего вектора. Это выгодно, потому что каждый из этих компонентов может быть вычислен независимо. Очевидно, если матрица хранится по столбцам, то должен быть использован следующий код:
DO J=l, N Kl = IA(J) K2 = IA(J+1)-1 Y(JA(K1:K2)) = Y(JA(K1:K2))+X(J)*A(Kl:K2) ENDDO
В каждой итерации цикла, произведение j-го столбца добавляется к результату, который, как предполагается, был первоначально равен нулю. Заметим также, что внешний цикл уже не пригоден для параллелизации. В качестве альтернативного распараллеливания можно попытаться разделить векторную операцию во внутреннем цикле. Внутренний цикл имеет всего несколько операций, так что это вряд ли даст значительное улучшение производительности. Это сравнение показывает, что структуры данных, возможно, придется изменить, чтобы улучшить производительность при работе с мощными компьютерами. Рассмотрим теперь усножение матриці на вектор для матриц, хранящихся в диагональном формате.
DO J=l, N JOFF = IOFF(J) DO 1=1, N Y(I) = Y(I) +DIAG(I,J)*X(J0FF+I) ENDDO ENDDO
Здесь каждый из диагоналей умножается на вектор х, и результат добавляется к вектору y. Снова предполагается, что вектор y заполнен нулями в начале цикла. С точки зрения распараллеливания и/или векторизации лучше использовать приведенный выше код, вероятно. С другой стороны, он носит недостаточно общий характер.
Решение нижне- или верхне- треугольной системы является еще одним важным "ядром" в вычислениях разреженных матриц . Следующий фрагмент кода показывает обычный блок решения нижнетреугольной системы Lx = у для формата хранения CSR.
X(l) = Y(l) DO I = 2, N Kl = IAL(I) K2 = IAL(I+1)-1 X(I)=Y(I)-D0TPRDDUCT(AL(K1:K2),X(JAL(K1:K2))) ENDDO
На каждом этапе вычисляется результат текущего умножения х на i-ю строку и вычитается из y(i). Это дает значение х(i). Функция Dotproduct вычисляет скалярное произведение двух произвольных векторов U(k1: k2) и V(k1: k2). Вектор AL (k1: k2) является i-й строкой матрицы L в разреженном формата и X (JAL (k1: k2)) является вектором компонентов X собранных в короткий вектор, который соответствует столбцовым индексам элементов в строке AL (k1: k2).