Ітераційні методи для розріджених систем лінійних рівнянь: Схеми зберігання, Базові операції на розріджених матрицях
Схеми зберігання розріджених матриць
Основною метою є збереження тільки ненульових елементів матрицею і можливість виробляти над нею найпростіші матричні операції.
Нехай Nz - кількість ненульових елементів матриці A. Розглянемо найбільш поширені схеми зберігання, більш детальна інформація наводиться в [4].
Найпростішою схемою зберігання розріджених матриць є так званий «координатний» формат. Структура вихідної матриці відображається на три одновимірних масиву: (1) - масив, що містить значення ненульових елементів матриці А в довільному порядку; (2) цілочисельний масив, що містить їх рядкові індекси; і (3) цілочисельний масив, що містить їх столбцовую індекси. Всі три вектора мають довжину Nz.
Приклад 1. Матриця
Може бути представлена як:
У наведеному прикладі елементи перераховані в довільному порядку, проте зазвичай вони перераховуються по рядках або стовпцях. Якщо елементи перераховуються по рядках, то вектор JC, який містить надлишкову інформацію, може бути замінений масивом покажчиків на початок кожного рядка. Це дасть помітну економію при зберіганні. Нова структура даних матиме три масиву, які будуть виконувати наступні функції:
- масив АА містить ненульові значення aij, які зберігаються через підрядник, від рядка 1 до n. Довжина масиву АА є Nz.
- цілочисельний масив JA містить столбцовую індекси елементів aij в матриці А. Довжина масиву JА є Nz.
Таким чином, представлена вище матриця А може бути збережена як:
Такий формат є найбільш популярним для зберігання розріджених матриць загального вигляду. Він носить назву розрідженого сатиричного формату або 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.
Тому для матриці з прикладу вище ці два масиви матимуть вигляд:
Зірочка вказує на використовуваний елемент. Зауважимо, що 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 =
Більш загальною схемою, популярною для векторних машин є так званий формат. Ця схема передбачає, що кількість діагоналей Nd невелика. Для зберігання потрібні два прямокутних масиву розміром n * Nd кожен.
У першому, COEF, подібному DIAG, зберігаються ненульові елементи матриці А. Ненульові елементи кожного рядка можуть зберігатися в рядку масиву COEF (1: n, 1: Nd), доповнюючи цей рядок нулями, якщо необхідно. Разом з COEF, цілочисельний масив JCOEF (1: n, 1: Nd) повинен зберігати індекс стовпця кожного елемента COEF.
Для матриці А з прикладу 2, схема зберігання Ellpack-Itpack набуде вигляду:
COEF =
Базові операції для розріджених матриць
Множення матриці на вектор - важлива операція, яка потрібна в більшості ітераційних алгоритмів розв'язання розріджених лінійних систем. У цьому розділі показано, як вона може бути реалізована для схем зберігання, наведених раніше.
У наступному прикладі коду на Fortran 90 представлений головний цикл множення матриці на вектор для матриць зберігаються в розрідженому столбцовую форматі:
Зверніть увагу, що кожна ітерація циклу обчислює різні компоненти результуючого вектора. Це вигідно, тому що кожен з цих компонентів може бути обчислений незалежно. Очевидно, якщо матриця зберігається за стовпцями, то повинен бути використаний наступний код:
У кожній ітерації циклу, твір j-го стовпця додається до результату, який, як передбачається, був спочатку дорівнює нулю. Зауважимо також, що зовнішній цикл вже не придатний для паралелізації. В якості альтернативного розпаралелювання можна спробувати розділити векторну операцію у внутрішньому циклі. Внутрішній цикл має всього кілька операцій, так що це навряд чи дасть значне поліпшення продуктивності. Це порівняння показує, що структури даних, можливо, доведеться змінити, щоб поліпшити продуктивність при роботі з потужними комп'ютерами. Розглянемо тепер усноженіе матриці на вектор для матриць, що зберігаються в діагональному форматі.
Тут кожен з діагоналей множиться на вектор х. і результат додається до вектору y. Знову передбачається, що вектор y заповнений нулями на початку циклу. З точки зору розпаралелювання і / або векторизації краще використовувати наведений вище код, ймовірно. З іншого боку, він носить недостатньо загальний характер.
Рішення нижне- або верхньо-трикутної системи є ще одним важливим "ядром" в обчисленнях розріджених матриць. Наступний фрагмент коду показує звичайний блок рішення ніжнетреугольной системи Lx = у для формату зберігання CSR.
На кожному етапі обчислюється результат поточного множення х на i -ю рядок і віднімається з y (i). Це дає значення х (i). Функція Dotproduct обчислює скалярний добуток двох довільних векторів U (k1: k2) і V (k1: k2). Вектор AL (k1: k2) є i-м рядком матриці L в розрідженому формату і X (JAL (k1: k2)) є вектором компонентів X зібраних в короткий вектор, який відповідає столбцовую індексам елементів в рядку AL (k1: k2).