Bc.JaroslavRyba Efektivn´ınásoben´ıˇr´ıdkýchmatic Diplomovápráce

(1)

(2)

(3)

Diplomov´ a pr´ ace

Efektivn´ı n´ asoben´ı ˇ r´ıdk´ ych matic

Bc. Jaroslav Ryba

Katedra softwarov´eho inˇzen´yrstv´ı

Vedouc´ı pr´ace: doc. Ing. Ivan ˇSimeˇcek, Ph.D.

7. kvˇetna 2019

(4)

(5)

Podˇ ekov´ an´ı

Chtˇel bych podˇekovat svému vedouc´ımu diplomové práce doc. Ing. Ivanu Simeˇˇ ckovi, Ph.D., za odborné veden´ı, za pomoc a rady pˇri zpracován´ı této práce. Dˇekuji také Mgr. Drahoslavˇe Rybové za pomoc pˇri gramatické kontrole práce.

(6)

(7)

Prohl´ aˇ sen´ı

Prohlaˇsuji, ˇze jsem pˇredloˇzenou práci vypracoval(a) samostatnˇe a ˇze jsem uvedl(a) veˇskeré pouˇzité informaˇcn´ı zdroje v souladu s Metodickým pokynem o etické pˇr´ıpravˇe vysokoˇskolských závˇereˇcných prac´ı.

Beru na vˇedom´ı, ˇze se na moji práci vztahuj´ı práva a povinnosti vyplývaj´ıc´ı ze zákona ˇc. 121/2000 Sb., autorského zákona, ve znˇen´ı pozdˇejˇs´ıch pˇredpis˚u.

V souladu s ust.§46 odst. 6 tohoto zákona t´ımto udˇeluji nevýhradn´ı oprávnˇen´ı (licenci) k uˇzit´ı této moj´ı práce, a to vˇcetnˇe vˇsech poˇc´ıtaˇcových program˚u, jeˇz jsou jej´ı souˇcást´ı ˇci pˇr´ılohou, a veˇskeré jejich dokumentace (dále souhrnnˇe jen

”D´ılo“), a to vˇsem osobám, které si pˇrej´ı D´ılo uˇz´ıt. Tyto osoby jsou oprávnˇeny D´ılo uˇz´ıt jakýmkoli zp˚usobem, který nesniˇzuje hodnotu D´ıla, a za jakýmkoli

´

uˇcelem (vˇcetnˇe uˇzit´ı k výdˇeleˇcným úˇcel˚um). Toto oprávnˇen´ı je ˇcasovˇe, teri- toriálnˇe i mnoˇzstevnˇe neomezené. Kaˇzdá osoba, která vyuˇzije výˇse uvedenou licenci, se vˇsak zavazuje udˇelit ke kaˇzdému d´ılu, které vznikne (byt’ jen zˇcásti) na základˇe D´ıla, úpravou D´ıla, spojen´ım D´ıla s jiným d´ılem, zaˇrazen´ım D´ıla do d´ıla souborného ˇci zpracován´ım D´ıla (vˇcetnˇe pˇrekladu), licenci alespoˇn ve výˇse uvedeném rozsahu a zároveˇn zpˇr´ıstupnit zdrojový kód takového d´ıla ale- spoˇn srovnatelným zp˚usobem a ve srovnatelném rozsahu, jako je zpˇr´ıstupnˇen zdrojový kód D´ıla.

V Praze dne 7. kvˇetna 2019 . . . .

(8)

Cesk´ˇ e vysok´e uˇcen´ı technick´e v Praze Fakulta informaˇcn´ıch technologi´ı

c 2019 Jaroslav Ryba. Vˇsechna pr´ava vyhrazena.

Tato práce vznikla jako ˇskoln´ı d´ılo na ˇCeském vysokém uˇcen´ı technickém v Praze, Fakultˇe informaˇcn´ıch technologi´ı. Práce je chránˇena právn´ımi pˇredpisy a mezinárodn´ımi úmluvami o právu autorském a právech souvisej´ıc´ıch s právem autorským. K jej´ımu uˇzit´ı, s výjimkou bezúplatných zákonných licenc´ı a nad rámec oprávnˇen´ı uvedených v Prohláˇsen´ı na pˇredchoz´ı stranˇe, je nezbytný sou- hlas autora.

Odkaz na tuto pr´aci

Ryba, Jaroslav. Efektivn´ı násoben´ı ˇr´ıdkých matic. Diplomov´a práce. Praha:

Cesk´ˇ e vysok´e uˇcen´ı technick´e v Praze, Fakulta informaˇcn´ıch technologi´ı, 2019.

(9)

Abstrakt

Tato diplomová práce se zabývá implementac´ı základu knihovny pro práci s ˇr´ıdkými maticemi. Dále je obsaˇzena implementace a optimalizace masivnˇe paralelizovaného násoben´ı matic na GPU za vyuˇzit´ı technologie CUDA.

Práce také slouˇz´ı k poskytnut´ı vhledu do problematiky násoben´ı matic, ˇr´ıdkých matic a efektivn´ı implementace algoritm˚u v technologii CUDA obecnˇe.

Kl´ıˇcová slova Násoben´ı matic, ˇr´ıdké matice, CUDA, optimalizace, knihovna, GPU, C++.

Abstract

This master’s thesis deals with implementation of the basis of sparse matrix library. It also contains implementation and optimalisations of massively parallel matrix multiplication on GPU in the CUDA technology.

This work is also intended to give basic level understanding of the matrix multiplication, sparse matrices and efficient implementation of algorithms under the limitations of the CUDA technology.

Keywords Matrix multiplication, sparse matrices, CUDA, optimization, library, GPU, C++.

(10)

(11)

Obsah

Uvod´ 1

1 Anal´yza 3

1.1 Z´akladn´ı pojmy . . . 3

1.2 Form´aty ˇr´ıdk´ych matic . . . 5

1.3 Algoritmy n´asoben´ı hust´ych matic . . . 11

1.4 Reˇserˇse existuj´ıc´ıch knihoven . . . 13

2 N´avrh 19 2.1 Interface . . . 19

2.2 Form´aty matic . . . 19

2.3 Sekvenˇcn´ı algoritmy . . . 20

2.4 Z´akladn´ı struktura paraleln´ıch algoritm˚u . . . 26

2.5 Pˇrechod mezi high level a low level funkcemi . . . 27

2.6 N´avrh tˇr´ıd . . . 29

3 Implementace 31 3.1 Implementace z´akladn´ı funkcionality . . . 31

3.2 Uˇzivatelsk´y interface . . . 31

3.3 Implementace CUDA . . . 32

3.4 Vyuˇzit´ı prostˇredk˚u CUDA . . . 33

3.5 Optimalizace . . . 33

4 Testov´an´ı 39 4.1 Testovac´ı hardware . . . 39

4.2 Nastaven´ı kompilace . . . 39

4.3 Unit testy . . . 40

4.4 Vˇecn´a spr´avnost implementace . . . 40

4.5 Testovac´ı data . . . 41

4.6 Vstupn´ı form´at matic . . . 41

(12)

4.7 Sekvenˇcn´ı verze . . . 41

4.8 Efektivita optimalizac´ı RC verze . . . 46

4.9 Efektivita optimalizac´ı CR verze . . . 48

4.10 Porovn´an´ı s cuSPARSE . . . 49

Z´avˇer 51

Literatura 53

5 Seznam pouˇzit´ych zkratek 55

6 Obsah pˇriloˇzen´eho CD 57

(13)

Seznam obr´ azk˚ u

1.1 Diagon´aln´ı form´at, pˇrevzato z [1] . . . 8

1.2 Vzorov´a matice z tabulky 1.1 ve form´atu ELL . . . 9

1.3 Vzorov´a matice z tabulky 1.1 v hybridn´ım form´atu . . . 10

2.1 Vysokoúrovˇnová funkce násoben´ı matic . . . 28

2.2 Z´akladn´ı n´avrh tˇr´ıd pro matice . . . 30

4.1 Graf z´avislosti doby bˇehu na velikosti matice s hustotou 10 % . . . 42

4.2 Graf závislosti doby bˇehu na hustotˇe matice o velikosti 1000×1000 42 4.3 Graf závislosti doby bˇehu na hustotˇe matice formátu CRS . . . 44

4.4 Graf z´avislosti doby bˇehu na hustotˇe matice form´atu ELL . . . 44

4.5 Graf z´avislosti doby bˇehu na hustotˇe matice form´atu BCRS . . . . 45

4.6 Graf z´avislosti doby bˇehu na hustotˇe matice form´atu DIA . . . 46

4.7 Graf efektivity optimalizac´ı RC verze . . . 47

4.8 Graf efektivity optimalizac´ı RC verze na BCRS . . . 47

4.9 Graf efektivity optimalizac´ı CR verze . . . 48

4.10 Porovn´an´ı efektivity optimalizovan´ych algoritm˚u s cuSPARSE . . . 49

(14)

(15)

Seznam tabulek

1.1 Vzorov´a hust´a matice . . . 5

1.2 Vzorov´a matice z tabulky 1.1 ve form´atu COO . . . 6

1.3 Vzorov´a matice z tabulky 1.1 ve form´atu CRS . . . 7

1.4 Vzorov´a matice z tabulky 1.1 ve form´atu CCS . . . 7

1.5 Vzorov´a matice z tabulky 1.1 ve form´atu BCRS . . . 7

4.1 Porovn´an´ı doby bˇehu v´ypoˇctu v ms . . . 49

(16)

(17)

Uvod ´

Jedn´ım ze základn´ıch algebraických problém˚u je násoben´ı matic. Tato operace je základn´ım stavebn´ım kamenem mnoha v praxi významných algoritm˚u, napˇr´ıklad fyzikáln´ıch simulac´ı, hledán´ı nejkratˇs´ı cesty. ˇCasto se také vyuˇz´ıvá v poˇc´ıtaˇcové grafice.

Reˇˇ sen´ı tohoto problému se vˇenovalo jiˇz znaˇcné úsil´ı a za posledn´ı stolet´ı se podaˇrilo sn´ıˇzit teoretickou ˇcasovou sloˇzitost vynásoben´ı dvou matic z na- ivn´ıho kubického algoritmu (O(n³)) aˇz na úroveˇn Coppersmith-Winogradova algoritmu (O(n^2.372)) a dosáhnout i znaˇcných zlepˇsen´ı efektivity implementac´ı tˇechto algoritm˚u.

Mnoho problém˚u vˇsak vede na násoben´ı ˇr´ıdkých matic, tedy matic s vel- kým mnoˇzstv´ım nulových prvk˚u. Klasické algoritmy násoben´ı matic nejsou schopné vyuˇz´ıt této vlastnosti matic a jejich ˇcasová sloˇzitost je nezávislá na mnoˇzstv´ı nenulových prvk˚u. Pˇri vyuˇzit´ı této vlastnosti je vˇsak moˇzné uˇsetˇrit znaˇcnou ˇcást výpoˇcetn´ıch prostˇredk˚u a dosáhnout tak mnohem vyˇsˇs´ı rychlosti výpoˇct˚u.

Vzhledem k neustálému zvyˇsován´ı výkonu grafických karet (GPU) a po- malému r˚ustu výkonu CPU se v posledn´ıch letech stále ˇcastˇeji výpoˇcty pˇresou- vaj´ı na GPU a masivnˇe paralelizuj´ı. Správnˇe implementovaný a optimalizo- vaný algoritmus tak m˚uˇze teoreticky z´ıskat pˇr´ıstup k ˇrádovˇe sto i tis´ıcinásob- ným výpoˇcetn´ım prostˇredk˚um. Tento trend je taktéˇz pos´ılen relativn´ı jednoduchost´ı implementace kód˚u za pomoc specializovaných technologi´ı (CUDA, openCV).

Tato práce se zamˇeˇruje na výˇse zm´ınˇený problém násoben´ı ˇr´ıdkých matic a vyuˇzit´ı výpoˇct˚u na grafické kartˇe pro jeho ˇreˇsen´ı. C´ılem práce je poskytnout ˇ

ctenáˇri vhled do základn´ı problematiky násoben´ı matic, problematiky ˇr´ıdkých matic a operac´ı s nimi a souˇcasného stavu existuj´ıc´ıch knihoven pro práci s ˇr´ıdkými maticemi.

Dále je souˇcást´ı práce navrˇzen´ı a implementace kostry knihovny pro práci s ˇr´ıdkými maticemi, která m˚uˇze v budoucnu slouˇzit jako základ pro testován´ı

(18)

Uvod´

algoritm˚u a pˇr´ıpadnˇe i jako základ vˇetˇs´ıch projekt˚u. V neposledn´ı ˇradˇe je také souˇcást´ı práce implementace a optimalizace nˇekolika implementac´ı násoben´ı ˇr´ıdkých matic v technologii CUDA (a integrace do vytvoˇrené kostry knihovny).

V posledn´ı ˇcásti práce jsou vytvoˇrené kódy testovány, vyhodnoceny a po- rovnány s efektivitou souˇcasné verze

”konkurenˇcn´ı“ knihovny pro n´asoben´ı ˇr´ıdk´ych matic na GPU.

Tato práce m˚uˇze slouˇzit jako studijn´ı materiál pro zájemce o problematiku ˇr´ıdkých matic, praktická ˇcást pak m˚uˇze být vyuˇzita jako základ pro vytvoˇren´ı v praxi pouˇzitelné implementace efektivn´ı knihovny pro práci s ˇr´ıdkými maticemi.

C´ılem prvn´ı kapitoly (Analýza) je seznámit ˇctenáˇre se základn´ı problematikou násoben´ı matic a algoritmy pouˇz´ıvanými pro násoben´ı matic v hustých formátech. Dále se v této kapitole také ˇctenáˇr doˇcte o formátech ukládán´ı ˇr´ıdkých matic, jejich specifikách a silných stránkách. V neposledn´ı ˇradˇe se zde také nacház´ı reˇserˇse souˇcasných knihoven pro práci s ˇr´ıdkými maticemi a informace o tom, jaké algoritmy a jaký interface tyto knihovny pouˇz´ıvaj´ı. Tyto informace budou v dalˇs´ıch ˇcástech vyuˇzity k návrhu nové knihovny.

Ve druhé kapitole (Návrh) je nejdˇr´ıve zd˚uvodnˇen výbˇer formát˚u ˇr´ıdkých matic pro implementaci, dále jsem navrˇzeny a popsány sekvenˇcn´ı algoritmy a datové struktury pro práci s vybranými formáty a navrˇzeny dva zp˚usoby paralelizace algoritm˚u násoben´ı matic. Je zde také popsán základ tˇr´ıdn´ı struktury navrˇzené knihovny a integrace paraleln´ıch algoritm˚u s knihovnou.

Ve tˇret´ı kapitole (Implementace) je popsán pr˚ubˇeh implementace knihovny a paralelizovaného násoben´ı. Je zde popsána ˇskála vyuˇzitých postup˚u pro op- timalizaci CUDA (paralelizovaných) kód˚u a racionalizace tˇechto technik. Jsou zde zm´ınˇena i nˇekterá omezen´ı této technologie a jak se s nimi vypoˇrádat.

V posledn´ı kapitole (Testován´ı) je popsán zp˚usob testován´ı správnosti a efektivity implementac´ı, vˇcetnˇe vyuˇzitého hardwaru a parametr˚u kompilace.

Dále se zde nacház´ı výsledky testován´ı efektivity hlavn´ıch optimalizaˇcn´ıch technik a verz´ı paralelizovaných algoritm˚u. V neposledn´ı ˇradˇe se zde taktéˇz nacház´ı srovnán´ı efektivity s knihovnou cuSPARSE.

(19)

Kapitola 1

Anal´ yza

C´ılem této kapitoly je seznámen´ı s problematikou, kterou se tato práce dále zabývá (násoben´ı matic, ˇr´ıdké matice), a s pˇr´ıstupem

”konkurenˇcn´ıch“ aplikac´ı k ˇreˇsen´ı t´eto problematiky.

1.1 Z´ akladn´ı pojmy

V této sekci jsou definovány nˇekteré významné term´ıny vyuˇz´ıvané dále v práci.

Nen´ı-li uvedeno jinak, je matematický základ ˇcerpán z [2], definice jsou vˇsak parafrázovány, upraveny a zjednoduˇseny pro úˇcely této práce.

Definice 1.1. Necht’ m, n ∈ N. Uspoˇrádaný soubor n·m ˇc´ısel zapsaný do tabulky o n ˇrádc´ıch a m sloupc´ıch nazýváme matice typu n×m. Matici obvykle znaˇc´ıme takto:

A=







a11 a12 . . . a1m

a₂₁ a₂₂ . . . a2m

... ... . .. ... a_n1 a_n2 . . . a_nm





 ,

kde a_ij jsou prvky matice. ˇC´ıslu i ˇr´ıkáme ˇrádkový a ˇc´ıslu j sloupcový index. [2]

Definice 1.2. Necht’m, n, i∈N. MaticiA typun×1 budeme nazývatslup- cový vektoro délcen. Matici typu 1×mbudeme nazývatˇrádkový vektor o délce m, zkr´acenˇe pak vektor o délce m, a jeho prvky a1j;j ≤ m m˚uˇzeme zkrácenˇe psátaj

Definice 1.3. Necht’ n∈N aA, B jsou vektory o d´elce n. Skal´arn´ı souˇcin vektor˚u c=A·B definujeme takto:c=^Pⁿ_k=1a_k·b_k.

(20)

1. Anal´yza

Definice 1.4. Necht’ m, n ∈ N, α ∈ R a A je matice typu n×m. Souˇcin matice A s re´aln´ym ˇc´ıslem α, znaˇc´ıme C =α·A, je matice typu n×m, pro kterou plat´ı

∀i, j∈N, i≤n, j≤m:c_ij =αa_ij.

Definice 1.5. Necht’m, n∈NaA, Bjsou matice typun×m.Souˇcet matic C = A+B je matice typu n×m, pro kterou plat´ı∀i, j ∈ N, i≤ n, j ≤ m : cij =aij +bij.

Definice 1.6. Necht’m, n∈NaAje matice typun×m.Transpozic´ı matice A, znaˇc´ımeA^T, naz´yv´ame matici typum×n, pro niˇz plat´ı

∀i, j ∈N, i≤m, j≤n:a^T_ij =aji.

Definice 1.7. Necht’ m, l, n∈N,A je matice typu n×l a B je matice typu l×m. Maticovým souˇcinem C = A·B nazýváme takovou matici typu n×m, pro kterou plat´ı

∀i, j∈N, i≤n, j ≤m:C_ij =

l

X

k=1

a_ik·b_kj.

Definice 1.8. Necht’ m, n, l ∈ N a A je matice typu n×m. Diagonálou matice Anazveme takový vektorD délkyl, pro který ∃c∈Ztak, ˇze plat´ı

(∀i∈N, i≤l:d_i =a_i,i+c, c≥0∧c+l=m)∨

(∀i∈N, i≤l:d_i=ai−c,i, c≤0∧c+l=n) . Hodnotucpak nazývámeoffsetem diagonály.

Definice 1.9. Rekneme, ˇˇ ze matice A = (a_i,j) je pásovou matic´ı, pokud existuj´ı nezáporné konstanty p, q(nazývané levý a pravý polo-pás) takové, ˇze ai,j 6= 0 pouze pokud i−p≤j≤i+q.[3]

Definice 1.10. Necht’m, n∈N,Aje matice typun×m.Hustotou matice Anazveme c∈R, pro nˇeˇz plat´ı

c= |{a_ij;i, j∈N∧i≤n∧j ≤m∧a_ij 6= 0}|

m·n .

Neexistuje oficiáln´ı definice ˇr´ıdké matice. Pro naˇse úˇcely vˇsak postaˇc´ı tato pragmatická definice od J. H. Wilkinsona:

Definice 1.11. Necht’m, n∈N, o maticiAtypun×mˇrekneme, ˇze jeˇr´ıdká, pokud m˚uˇzeme vyuˇz´ıt faktu, ˇze ˇcást jej´ıch hodnot je nulová.

Definice 1.12. Necht’ f, g jsou funkce na pˇrirozen´ych ˇc´ıslech. ˇRekneme, ˇzeg jeasymptotickou horn´ı mez´ı f, znaˇc´ımef =O(g) pr´avˇe tehdy, kdyˇz:

∃c, n₀;c∈R⁺, n₀ ∈N:∀n∈N, n≥n₀ :f(n)≤c·g(n).

(21)

1.2. Formáty ˇr´ıdkých matic Definice 1.13. Necht’ f, g jsou funkce na pˇrirozených ˇc´ıslech a x ∈ N. Vý- poˇcetn´ı sloˇzitost´ı implementace algoritmu I nazveme takovou funkcif(x), pro niˇz plat´ıg =O(f) a g(x) je funkc´ı z´avislosti poˇctu provedených operac´ı na velikosti vstupu.

Definice 1.14. Necht’f, gjsou funkce na pˇrirozených ˇc´ıslech ax∈N.Pamˇe- t’ovou nároˇcnost´ı implementace algoritmu I nazveme takovou funkci f(x) pro niˇz plat´ıg=O(f) ag(x) je funkc´ı závislosti mnoˇzstv´ı vyuˇzité pamˇeti na velikosti vstupu.

Definice 1.15. ˇCasovou efektivitouE implementace algoritmu I (téˇz zkrácenˇe efektivitou algoritmu, nebo výkonem algoritmu) na hardwaru H se vstupem V rozum´ıme mnoˇzstv´ı operac´ı typu T, kter´e pˇri spuˇstˇen´ıI na H pr˚umˇernˇe probˇehne za jednotku ˇcasu.

Pro úˇcely této práce (násoben´ı ˇr´ıdkých matic) pak budouT z 1.15 operace s reálnými ˇc´ısly (plovouc´ıˇrádovou ˇcárkou), konkrétnˇe jejich sˇc´ıtán´ı a násoben´ı.

Tuto hodnotu z´ıskáme jako pod´ıl provedených operac´ı a ˇcasu spotˇrebovaného na proveden´ı výpoˇctu. Pro tuto veliˇcinu budeme pouˇz´ıvat jednotku FLOPS (poˇcet operac´ı s ˇrádovou ˇcárkou za sekundu).

1.2 Form´ aty ˇ r´ıdk´ ych matic

V této kapitole pop´ıˇsi nˇekteré základn´ı zp˚usoby ukládán´ı ˇr´ıdkých matic. Pro ilustraci ukládán´ı bude vyuˇzita jednoduchá matice 4x4, kterou m˚uˇzete vidˇet v tabulce 1.1.

Za desetilet´ı vývoje vzniklo nespoˇcet r˚uzných formát˚u ukládán´ı matic.

C´ılem této ˇcásti tedy nen´ı popis vˇsech, ale nˇekolika vybraných. Formáty byly vyb´ırány dle jejich významnosti (rozˇs´ıˇren´ı) a relevantnosti k dalˇs´ımu obsahu této práce.

Základn´ı ukládán´ı matice v hustém formátu je zˇrejmé a nebude v této práci dále rozeb´ıráno.

Tabulka 1.1: Vzorov´a hust´a matice







0 1.1 0 2.0

2.3 0 0 2.4

0 0 1.0 0

0 0 0 0.4







1.2.1 COO

Informace v t´eto kapitole vych´az´ı z [4].

(22)

1. Anal´yza

Nejintuitivnˇejˇs´ım zp˚usobem ukládán´ıˇr´ıdkých matic je pomoc´ı pozic a hodnot nenulových prvk˚u (ˇrádek, sloupec, hodnota). Tento formát se obvykle oznaˇcuje jako COO (coordinate format - souˇradnicový formát). Pro lepˇs´ı práci s pamˇet´ı (cache) je vhodné ukládat tyto hodnoty pomoc´ı tˇr´ı pol´ı. Dle zp˚usobu vyuˇzit´ı se také m˚uˇze hodit hodnoty seˇradit dle pozice.

Dˇr´ıve zm´ınˇená vzorová matice (1.1) se dá t´ımto zp˚usobem zapsat tak, jak je vidˇet v tabulce 1.2.

Tabulka 1.2: Vzorov´a matice z tabulky 1.1 ve form´atu COO

ˇr´adek 0 0 1 1 2 3

sloupec 1 3 0 3 2 3

hodnota 1.1 2.0 2.3 2.4 1.0 0.4

Lze snadno nahl´ednout, ˇze pro maticin×ms nnz nenulov´ymi hodnotami k uspoˇren´ı m´ısta dojde, je-li

nnz·(valueSize+ 2·indexSize)< n·m·valueSize,

kde valueSize je poˇcet bit˚u k uloˇzen´ı hodnoty a indexSize je poˇcet bit˚u na uloˇzen´ı indexu. Pokud dále pˇredpokládámevalueSize=indexSize(napˇr´ıklad obˇe hodnoty 32, coˇz je standardn´ı velikost), dosáhneme úspory, je-li ^nnz_n·m < ¹₃. 1.2.2 CRS

Definice v t´eto kapitole ˇc´asteˇcnˇe pˇrevzata z [5] a rozˇs´ıˇrena o doplˇnuj´ıc´ı informace z [6] a [4].

Compressed row storage (CRS) lze oznaˇcit za inkrementáln´ı vylepˇsen´ı COO. Hodnoty prvk˚u i sloupcové indexy jsou ukládány stejným zp˚usobem (2 vektory), zde je vˇsak jiˇz pˇr´ımo vyˇzadováno jejich seˇrazen´ı (dle ˇrádk˚u a sloupc˚u v tomto poˇrad´ı d˚uleˇzitosti). ˇRádkový vektor pak obsahuje informaci o poˇctu prvk˚u v jednotlivých ˇrádc´ıch. Pro urychlen´ı pˇr´ıstupu k prvk˚um (aby nemusel být tento vektor procházen od zaˇcátku pˇri ˇcten´ı prvku) je pak tato informace uloˇzena ve formˇe prefixového souˇctu poˇct˚u prvk˚u (tedy kolikátý prvek je prvn´ı v tomto ˇrádku, poˇc´ıtáno od nuly), tyto hodnoty budou dále nazývány offset ˇrádku. Je standardem (pro zjednoduˇsen´ı algoritm˚u) rozˇs´ıˇrit vektor offset˚u ˇrádk˚u o jeden prvek, kterým je poˇcet nenulových prvk˚u v matici, tedy offset n+1. (neexistuj´ıc´ıho) ˇrádku.

Pamˇet’ová nároˇcnost tohoto formátu (se stejným znaˇcen´ım a pˇredpoklady jako u COO) je 2nnz +n+ 1, coˇz je u matice bez prázdných ˇrádk˚u hodnota menˇs´ı nebo rovna pamˇet’ové nároˇcnosti COO (pokud bereme v úvahu, ˇze i u COO je nezbytné uloˇzit poˇcet nenulových prvk˚u). ˇC´ım hustˇs´ı je matice, t´ım vˇetˇs´ı výhody CRS pˇrináˇs´ı, k rovnosti dojde u matice s pˇresnˇe jedn´ım prvkem na ˇrádek.

(23)

1.2. Formáty ˇr´ıdkých matic Tabulka 1.3: Vzorová matice z tabulky 1.1 ve formátu CRS

ofset ˇr´adku 0 2 4 5 6

sloupec 1 3 0 3 2 3

hodnota 1.1 2.0 2.3 2.4 1.0 0.4 1.2.3 CCS

Z´akladn´ı popis pˇrebr´an z [6].

Compressed column storage (CCS), také nazývaný Harwell-Boeing˚uv for- mát ˇr´ıdkých matic, je analogický formát k CRS, kde m´ısto po ˇrádc´ıch je vektor procházen po sloupc´ıch a jsou zaznamenávány ˇrádkové indexy hodnot. Tento formát se dá také chápat jako CRS transponované matice.

Jak lze z výˇse uvedeného popisu pochopit, vlastnosti tohoto formátu jsou stejné jako u CRS a tento formát m˚uˇze být vhodný pro násoben´ı s matic´ı v CRS formátu.

Tabulka 1.4: Vzorov´a matice z tabulky 1.1 ve form´atu CCS ofset sloupce 0 1 2 3 6

ˇr´adek 1 0 2 0 1 3

hodnota 2.3 1.1 1.0 2.0 2.4 0.4

1.2.4 BCRS

Definice vyuˇz´ıv´a informace z [5] s rozˇs´ıˇren´ım o [6] a [4].

Block compressed row storage (BCRS, nˇekdy také BSR), je nadstavba nad CRS, která vyuˇz´ıvá vlastnosti matic maj´ıc´ıch nenulové hodnoty shlukované (prostorovˇe) bl´ızko u sebe.

P˚uvodn´ı matice je nejdˇr´ıve rozloˇzena na podmatice (bloky) ˇrádovˇe menˇs´ı velikosti (napˇr´ıklad 2x2) a následnˇe je tato matice blok˚u uloˇzena pomoc´ı CRS, pˇriˇcemˇz blok s pouze nulovými hodnotami je brán jako nulový element. Pˇri ukládán´ı hodnot je po ˇrádc´ıch zapsán celý blok, nikoliv jen jeho nenulové elementy.

Pro velikost bloku 2x2 n´am tedy pro naˇsi vzorovou matici (1.1) vznikne z´apis v tabulce 1.5.

Tabulka 1.5: Vzorová matice z tabulky 1.1 ve formátu BCRS ofset ˇrádku 0 2 3

sloupec 0 1 1

hodnota 0 1.1 2.3 0 0 2.0 0 2.4 1.0 0 0 0.4

(24)

1. Anal´yza

Adresovat je pak nutné dvoufázovˇe (vyhledat adresu bloku a následovnˇe pˇriˇc´ıst pozici v bloku), coˇz ovˇsem na modern´ım hardwaru nepˇridává prakticky ˇzádnou výpoˇcetn´ı sloˇzitost. Pˇr´ınosy ukládán´ı matic v tomto formátu jsou velmi závislé na struktuˇre matice a vhodnˇe zvolené velikosti bloku. ˇC´ım v´ıce je matice tvoˇrená shluky hodnot, t´ım lepˇs´ıch výsledk˚u lze dosáhnout.

Aˇckoliv tato metoda ukládán´ı m˚uˇze (kv˚uli ukládán´ı i nulových element˚u) vést k navýˇsen´ı pamˇet’ových nárok˚u oproti CRS, je v praxi ˇcasto vyuˇz´ıvána pro lepˇs´ı vyuˇzit´ı cache a obecnˇe lepˇs´ı výkonnost nˇekterých algoritm˚u.

1.2.5 DIA

Definice parafr´azov´ana z [1].

Diagonal format (DIA), také známý jako compressed diagonal storage (CDS), je vhodný, pokud jsou nenulové hodnoty soustˇredˇeny do malého mnoˇz- stv´ı diagonál. Matice je uloˇzena za pomoci dvou pol´ı: pole hodnot a fset˚u diagonál od hlavn´ı diagonály. Diagonály nad a pod hlavn´ı diagonálou maj´ı negativn´ı respektivˇe pozitivn´ı offsety. Do pole hodnot jsou hodnoty ukládány postupnˇe od diagonály s nejniˇzˇs´ım offsetem k nejvyˇsˇs´ımu a v rámci diagonál dle sloupcového indexu. Ukládány jsou jsou i nulové hodnoty a diagonály jsou doplnˇeny jako by zaˇc´ınaly na prvn´ım a konˇcily na posledn´ım ˇrádku.

Ukázku tohoto formátu m˚uˇzete vidˇet na obrázku 1.1. Neexistuj´ıc´ı prvky (*) mohou být nahrazeny napˇr´ıklad 0.

Obrázek 1.1: Diagonáln´ı formát, pˇrevzato z [1]

A=







1 7 0 0 0 2 8 0 5 0 3 9 0 6 0 4







DIA

hodnoty * * 5 6 1 2 3 4 7 8 9 *

offsety diagon´al -2 0 1

Výhodou tohoto formátu je vysoká úspornost pro specifické typy matic a efektivn´ı paralelizace nˇekterých maticových algoritm˚u (napˇr. násoben´ı matice vektorem). Nevýhodou je vysoký nár˚ust nároˇcnosti (pamˇet’ové i výpoˇcet- n´ı) pro matice nemaj´ıc´ı vhodnou strukturu a s t´ım spojená omezená kategorie matic, na nˇeˇz lze formát efektivnˇe aplikovat.

1.2.6 ELL

Definice parafr´azov´ana z [1].

(25)

1.2. Formáty ˇr´ıdkých matic ELLPACK (ELL) je dalˇs´ı specializovaný formát. Podm´ınka tohoto formátu je vˇsak volnˇejˇs´ı neˇz u DIA. M´ısto n´ızkého poˇctu diagonál zde poˇzadujeme, aby byly prvky rozloˇzeny mezi ˇrádky co nejrovnomˇernˇeji. Pˇresnˇeji chceme, aby ˇrádek s nejv´ıce nenulovými elementy byl co nejkratˇs´ı.

Matice o n ˇrádc´ıch s nejv´ıce m nenulovými prvky na ˇrádek je uloˇzena pomoc´ı dvou hustých n×m matic. V prvn´ı z tˇechto matic jsou ukládány nenulové hodnoty v poˇrad´ı dle index˚u sloupc˚u, tedy jako kdyby z p˚uvodn´ı matice byly vypuˇstˇeny nuly a ostat´ı hodnoty

”sraˇzeny“ doleva. Do druhé matice jsou ukládány p˚uvodn´ı sloupcové indexy tˇechto prvk˚u. Obˇe tyto matice jsou zprava doplnˇeny nulami (aby byla zachována stejná délka vˇsech ˇrádk˚u).

Obrázek 1.2: Vzorová matice z tabulky 1.1 ve formátu ELL

hodnoty=







1.1 2.0 2.3 2.4 1.0 0 0.4 0







indexy=





 1 3 0 3 2 0 3 0







Dˇr´ıve zm´ınˇená vzorová matice (1.1) se dá v ELL zapsat tak, jak je vidˇet na obrázku 1.2. I v tomto pˇr´ıpadˇe mohou být samozˇrejmˇe pro vyˇsˇs´ı efektivitu matice rozepsány do dvou 1D pol´ı.

ELL se dá povaˇzovat za rozˇs´ıˇren´ı DIA. Zat´ımco u DIA je pozice prvku implicitn´ı (dle pozice diagonály), u ELL je explicitnˇe urˇcená, coˇz umoˇzˇnuje jej vyuˇz´ıt pro mnohem ˇsirˇs´ı ˇskálu matic.

Tento formát je také velmi vhodný pro vektorové architektury a paralelizaci a efektivn´ı práci s cache.

1.2.7 Hybridn´ı formát Tato ˇcást ˇcerpá z [1].

Hybridn´ı formát (HYB) je kombinace ELL a COO. C´ılem vyuˇzit´ı tohoto formátu je zvýˇsit pouˇzitelnost ELL na obecnˇejˇs´ı pˇr´ıpady a zároveˇn zachovat (ˇcásteˇcnˇe) jeho efektivitu.

Prvn´ıchN prvk˚u kaˇzdého ˇrádku se uloˇz´ı pomoc´ı ELL, zbývaj´ıc´ı prvky se uloˇz´ı pomoc´ı COO. C´ılem je vytvoˇrit jakési

”jádro“ matice, které bude uloˇzené v ELL a bude obsahovat pouze malé mnoˇzstv´ı nulových hodnot.

Sloˇzitˇejˇs´ım problémem je urˇcen´ı vhodného N. Tato hodnota závis´ı na struktuˇre matice a výkonnosti algoritm˚u nad COO a ELL. V´ıme-li tento pomˇer výkonnost´ı, m˚uˇzeme tuto hodnotu jiˇz pˇresnˇe urˇcit výpoˇctem z histogramu poˇct˚u nenulových element˚u na ˇrádek.

Pˇr´ıklad (vzorové) matice uloˇzené v hybridn´ım formmátu (s N = 1) si m˚uˇzete prohlédnout na obrázku 1.3.

(26)

1. Anal´yza

Obrázek 1.3: Vzorová matice z tabulky 1.1 v hybridn´ım formátu

hodnotyELL=





 1.1 2.3 1.0 0.4







indexyELL=





 1 0 2 3







ˇr´adek 0 1

sloupec 3 3

hodnota 2.0 2.4 1.2.8 Dynamick´e form´aty

Tato ˇc´ast ˇcerp´a z [7] a [8].

Vˇetˇsina formát˚u ukládán´ı ˇr´ıdkých matic je velmi neefektivn´ı pro zmˇeny matice. Zvláˇstˇe neefektivn´ı jsou operace pˇridáván´ı a odeb´ırán´ı hodnot, které u vˇetˇsiny formát˚u ˇcasto vynut´ı pˇrepis celé matice (coˇz je operace se sloˇzitost´ı O(nnz)). Dynamick´e formáty vycház´ı ze základn´ıch formát˚u pro ukládán´ı ˇr´ıdkých matic (nejˇcastˇeji CRS, BCRS a COO) a snaˇz´ı se odstranit tento problém pˇri zachován´ı n´ızké pamˇet’ové nároˇcnosti.

Dynamické formáty mohou být realizovány pomoc´ı:

• Spojových seznam˚u (jednosmˇerných, obousmˇerných i cyklických) - tento zp˚usob umoˇzˇnuje velmi rychlé pˇridáván´ı a odeb´ırán´ı prvk˚u a je snadno implementovatelný (ˇcasto pouˇz´ıván pro COO).

• Stromových struktur (napˇr´ıklad trie) - tento zp˚usob je vhodný pro ˇcasté ˇcten´ı hodnot matice s náhodným pˇr´ıstupem, je vˇsak významnˇe kom- plikovanˇejˇs´ı na implementaci a zajiˇstˇen´ı efektivity algoritm˚u pracuj´ıc´ıch s celou matic´ı

• Pole pol´ı (pole pointer˚u) - pole hodnot (index˚u) je rozdˇeleno na ˇcásti (obvykle samostatné ˇrádky/sloupce), které jsou ukládány do samostatnˇe alokovaných pol´ı a m´ısto offset˚u ˇrádk˚u je uchováván celý pointer, coˇz umoˇzˇnuje minimalizovat dopad pˇridán´ı prvk˚u pouze na jeden ˇrádek (O(nnz_m)) m´ısto celé matice (ˇcasto pouˇz´ıván pro CRS, BCRS).

• Rozvolnˇen´ı formát˚u - mezi uloˇzené hodnoty jsou pˇridány výplˇnové prvky (nulové), které jsou pˇrepsány v pˇr´ıpadˇe pˇridán´ı nových prvk˚u a v pˇr´ıpadˇe odeb´ırán´ı prvk˚u jsou prvky nahrazeny výplˇnovými hodnotami.

Pˇr´ıpadnˇe mohou být vyuˇzity kombinace nˇekolika z tˇechto technik (napˇr´ı- klad dále zm´ınˇený Eigen vyuˇz´ıvá 3. a 4.).

Výhodami tˇechto formát˚u jsou flexibilita pro vkládán´ı a odeb´ırán´ı prvk˚u, pˇr´ıpadnˇe rychlejˇs´ı pˇr´ıstup k nim. Nevýhodami pak jsou obt´ıˇznost zachován´ı

(27)

1.3. Algoritmy násoben´ı hustých matic pamˇet’ové lokality (vyuˇzit´ı cache), obt´ıˇznost vektorizace/paralelizace, vyˇsˇs´ı pamˇet’ová nároˇcnost a obecnˇe niˇzˇs´ı efektivita pro algoritmy pracuj´ıc´ı systema- ticky po ˇrádc´ıch/sloupc´ıch (mezi které patˇr´ı i násoben´ı matic).

1.3 Algoritmy n´ asoben´ı hust´ ych matic

1.3.1 Klasick´y algoritmus

Základn´ı (a pro mnohé pˇr´ıpady i nejlepˇs´ı) postup vyplývá ze vzorce pro násoben´ı matic a m˚uˇzete jej vidˇet v algoritmu 1.

Algoritmus 1 Klasick´e n´asoben´ı matic

1: procedure NaiveMultiply(A, B) . Vyn´asob A*B

2: C←zeroes(A.height, B.width) .Vytvoˇr 2D pole na ukládán´ı výsledk˚u a nastav jej na nuly

3: fori←0 to A.heightdo

4: for j←0 to B.width do

5: fork←0 to A.width do

6: C[i][j]←C[i][j] +A[i][k]∗B[k][j]

7: return C . Matice s v´ysledky

Pro zvýˇsen´ı efektivity je samozˇrejmˇe moˇzné vhodnˇe zmˇenit poˇrad´ı prochá- zen´ı a pˇr´ıpadnˇe vyuˇzit transpozice.

Výhodami tohoto algoritmu jsou jednoduchost jeho implementace, snad- nost paralelizace, vyuˇzitelnost pro specializované formáty ukládán´ı (ˇr´ıdké matice) a vysoká efektivita pro relativnˇe malé matice.

Jeho nevýhodou je pak jeho asymptotická sloˇzitost, která je Θ(n³), a z toho plynouc´ı pomalost (v˚uˇci dále zm´ınˇeným algoritm˚um) pˇri aplikaci na velké matice (ˇrádovˇe 1000×1000 prvk˚u).

1.3.2 Strassen˚uv algoritmus

Strassen˚uv algoritmus [9] je zaloˇzen na principu rozdˇel a panuj (rekurzivn´ı pˇr´ıstup). Asymptotick´a sloˇzitost tohoto algoritmu je O(n^log²⁷) ≈ O(n^2.808).

D˚ukaz správnosti algoritmu je po dosazen´ı do výsledného vzorce zˇrejmý, a pro- to nebude rozvádˇen.

Chceme-li vypoˇc´ıtat maticov´y souˇcinC =A·B, rozdˇel´ıme nejdˇr´ıve kaˇzdou matici na 4 stejnˇe velk´e podmatice

C=

"

A11 A12

A21 A22

#

·

"

B11 B12

B21 B22

# .

(28)

1. Anal´yza

Následnˇe je nutné vypoˇc´ıtat hodnoty 7 mezivýsledk˚u M₁=(A₁₁+A₂₂)·(B₁₁+B₂₂), M2=(A21+A22)·B11,

M3=A11·(B12−B22), M4=A22·(B21−B11), M₅=(A₁₁+A₁₂)·B₂₂,

M₆=(A₂₁−A₁₁)·(B₁₁+B₁₂), M₇=(A₁₂−A₂₂)·(B₂₁+B₂₂).

Nakonec se za pouˇzit´ı tˇechto pomocných matic vypoˇc´ıtá výsledná matice C=

"

M₁+M₄−M₅+M₇ M₃+M₅ M2+M4 M1−M2+M3+M6

# .

Celkovˇe je tedy nezbytné pro jednu úroveˇn Strassenova algoritmu provést 18 souˇct˚u (Θ(n²)) a 7 násoben´ı (Θ(n³)). Pro dostateˇcnˇe velké matice je jedna

´

uroveˇn Strassenova algoritmu oproti klasické variantˇe rekurzivn´ıho násoben´ı matic, která vyuˇz´ıvá 4 souˇcty a 8 násoben´ı, tedy aˇz o 12.5 % efektivnˇejˇs´ı.

Nevýhodou tohoto algoritmu je n´ızká efektivita na menˇs´ıch matic´ıch (kv˚uli vˇetˇs´ımu poˇctu sˇc´ıtán´ı a celkové vˇetˇs´ı komplexnosti) a zvýˇsená pamˇet’ová nároˇc- nost. I bez pˇrepisován´ı vstupn´ıch matic lze vˇsak algoritmus provést s pamˇe- t’ovou sloˇzitost´ıO(²ⁿ₃ ) [10].

1.3.3 Winogradova varianta Strassenova algoritmu

Winogradova varianta Strassenova algoritmu [11] je algoritmus zaloˇzený na stejném principu jako p˚uvodn´ı Strassen˚uv algoritmus, který vˇsak m´ısto 18 vyuˇz´ıvá pouze 15 souˇct˚u. Aˇckoliv toto zlepˇsen´ı nemá ˇzádný vliv na asymptotickou sloˇzitost, v praktické implementaci se projev´ı drobným zlepˇsen´ı efektivity algortimu.

Cel´e sch´ema algoritmu je (znaˇcen´ıA, B, C pˇrevzato z vysvˇetlen´ı Strasse- nova algoritmu):

S1 =A21+A22, M1 =S2·S6, V1=M1+M2, S₂ =S₁−A₁₁, M₂ =A₁₁·B₁₁, V₂=V₁+M₄, S₃ =A₁₁−A₂₁, M₃ =A₁₂·B₂₁, C₁₁=M₂+M₃, S4 =A12−S2, M4 =S3·S7, C12=V1+M5+M6, S₅ =B₁₂−B₁₁, M₅ =S₁·S₅, C₂₁=V₂−M₇, S₆ =B₂₂−S₅, M₆ =S₄·B₂₂, C₂₂=V₂+M₅, S7 =B22−B12, M7 =A22·S8,

S8 =S6−B21.

Napˇr´ıklad vyuˇzit´ım 4 iterac´ı tohoto algoritmu (pro matice 16384×16384) lze prakticky dosáhnout sn´ıˇzen´ı poˇctu aritmetických operac´ı o 41.3 % oproti klasickému algoritmu [10].

(29)

1.4. Reˇserˇse existuj´ıc´ıch knihoven Pamˇet’ová nároˇcnost z˚ustává stejná jako u p˚uvodn´ıho Strassenova algoritmu[11].

1.3.4 Dalˇs´ı algoritmy

V souˇcasn´e dobˇe existuje jiˇz ˇrada algoritm˚u s niˇzˇs´ı asymptotickou sloˇzitost´ı neˇz Strassen˚uv algoritmus (O(n^w=2.808)). V roce 1978 vznikl Pan˚uv algoritmus (w = 2.796), o rok pozdˇeji algoritmus od Bimi a spol. (w = 2.78).

Následovala ˇrada dalˇs´ıch algoritm˚u, které pˇrekonávaly do té doby nejlepˇs´ı dosaˇzené výsledky. Z dalˇs´ıch známých algoritm˚u stoj´ı za zm´ınku napˇr´ıklad Coppersmith-Winograd˚uv algoritmus (w = 2.376) z roku 1989 (pozdˇeji sn´ı- ˇ

zeno naw= 2.372), který je dodnes nejrychlejˇs´ım známým algoritmem.[12]

Aˇckoliv tyto algoritmy jsou asymptoticky rychlejˇs´ı, matice, pro nˇeˇz jsou tyto algoritmy oproti klasickému (a také Strassenovu) algoritmu výhodné, jsou natolik velké, ˇze nejsou zpracovatelné na souˇcasném hardwaru. Zvláˇstˇe pak mnohé z novˇejˇs´ıch algoritm˚u jsou pouze teoretického rázu a v praxi (napˇr´ıklad kv˚uli nepˇresnostem pˇri výpoˇctech u ˇc´ısel s plovouc´ı desetinnou ˇcárkou) ne- pouˇzitelné.

1.3.5 Vyuˇzit´ı subkubick´ych algoritm˚u pro ˇr´ıdk´e matice

Pro implementaci Strassenova algoritmu (a jeho variant), coˇz je jedin´y v sou- ˇ

casnosti pouˇz´ıvaný subkubický algoritmus, je nutné dokázat matici horizon- tálnˇe a vertikálnˇe rozdˇelit. Tato operace je u husté matice jednoduchá (nemus´ı probˇehnout ˇzádné ˇcten´ı, staˇc´ı upravit zp˚usob indexován´ı do matice a meze), na druhou stranu u ˇr´ıdkých matic se jedná o operaci velmi sloˇzitou, s ˇcasovou nároˇcnost´ıO(nnz), tedy pˇreˇcten´ı celé matice a jej´ı znovuzapsán´ı do nˇekolika nových, coˇz znaˇcnˇe zmenˇsuje pˇr´ınos Strassenova algoritmu i u sekvenˇcn´ıch implementac´ı.

Efektivita násoben´ı matic na grafické kartˇe (masivnˇe paralelizovaná) je obvykle mnohem v´ıce omezena vyuˇzit´ım pamˇeti, neˇz výkonem výpoˇcetn´ıch jed- notek. Pro masivn´ı paralelizaci je tedy tento algoritmus, který znaˇcnˇe zvyˇsuje poˇcet pamˇet’ových operac´ı a mnoˇzstv´ı vyuˇzité pamˇeti, naneˇstˇest´ı jeˇstˇe ménˇe vhodný, neˇz pro sekvenˇcn´ı verzi a proto nebude jeho praktická implementace uvaˇzována.

1.4 Reˇ serˇ se existuj´ıc´ıch knihoven

1.4.1 NIST Sparse BLAS

Informace v této sekci vycházej´ı z uˇzivatelského manuálu [13] a zkoumán´ı kód˚u dostupných z [14].

(30)

1. Anal´yza

Touto knihovnou se budu nejv´ıce zabývat nejen kv˚uli jej´ı rozˇs´ıˇrenosti, ale zároveˇn protoˇze jej´ı kód je open-source, a lze ji proto snadno a d˚ukladnˇe analyzovat.

1.4.1.1 Interface

Celou funkcionalitu t´eto knihovny lze rozdˇelit do 4 kategori´ı:

• Level 1 - operace typuvektor×vektor

• Level 2 - operace typumatice×vektor

• Level 3 - operace typumatice×matice

• Ostatn´ı - inicializace, destrukce, getry, setry

Z výˇse uvedených kategori´ı je pro tuto práci nejvýznamnˇejˇs´ı Level 3, protoˇze právˇe ten obsahuje násoben´ı matic. Dále stoj´ı za zm´ınku také ˇzivotn´ı cyklus matic.

Práce s ˇr´ıdkými maticemi sestává z tˇechto tˇr´ı krok˚u:

• Vytvoˇren´ı matice a z´ısk´an´ı handleru k n´ı (int)

• Pouˇzit´ı handleru jako argumentu pˇri v´ypoˇctech s matic´ı

• Kdyˇz matice jiˇz nen´ı potˇrebná, je explicitnˇe (uˇzivatelem) zavolána ˇcistic´ı funkce, která uvoln´ı vyuˇz´ıvané prostˇredky

Pro nás nejvýznamnˇejˇs´ı funkce násoben´ı matic má pak základn´ı tvary:

C ←α·A·B+C, C ←α·A^T ·B+C, a pro matice komplexn´ıch ˇc´ısel

C ←α·A^H ·B+C,

kde A je ˇr´ıdká matice, B, C jsou husté matice a αje reálné ˇc´ıslo (float). Starˇs´ı verze podporovaly nav´ıc koeficientβ pro sˇc´ıtán´ı sC. Tato funkce pak existuje ve 4 variantách pro r˚uzné typy prvk˚u matic (jednoduchá, dvojitá pˇresnost atd.).

Násoben´ı dvou ˇr´ıdkých matic knihovna neumoˇzˇnuje. Pˇri volán´ı funkc´ı jsou ˇr´ıdké matice reprezentovány handlerem, husté matice 1D polem (pointer), zp˚usobem, jak jsou uloˇzeny (po ˇrádc´ıch, sloupc´ıch - enum), a velikost´ı jejich hlavn´ı dimenze (int). Výbˇer varianty s transpozic´ı, ˇci konjungovanou transpozic´ı je reprezentován enumem.

(31)

1.4. Reˇserˇse existuj´ıc´ıch knihoven Tato knihovna vnitˇrnˇe (nedeklarováno v hlaviˇckových souborech) repre- zentuje matice pomoc´ı formátu BCRS, konkrétn´ı parametry se nedaj´ı z vnˇej- ˇsku ovlivnit. Dále také umoˇzˇnuje specializované algoritmy a ukládán´ı pro matice s význaˇcnými vlastnostmi (symetrické, triangulárn´ı).

Knihovna NIST Sparse BLAS vznikala p˚uvodnˇe v jazyce Fortran, coˇz je poznat i ze zp˚usobu práce s maticemi a z interface funkc´ı. Jak bylo vˇsak jiˇz zm´ınˇeno v [15], tento inferface má pro vyuˇzit´ı v C++ mnoho problém˚u a pro

´

uˇcely této práce by jej bylo vhodné vylepˇsit. Zejména se jedná o ˇzivotn´ı cyklus matic. Zde by bylo vhodné vyuˇz´ıt moˇznost´ı C++ práce s objekty a drˇzet se obvyklejˇs´ıch standard˚u pro jejich vytváˇren´ı a niˇcen´ı (konstruktory, destruktory). Dále pak veˇskeré parametry týkaj´ıc´ı se jedné matice (poˇrad´ı prvk˚u, délka dimenze, data) by mˇely být sjednoceny do jednoho objektu. Bude-li pod- porováno v´ıce typ˚u prvk˚u, mˇely by nav´ıc být funkce pˇret´ıˇzeny (nebo vyuˇzity ˇsablony) m´ısto nˇekolika r˚uzných jmen funkc´ı.

1.4.1.2 Algoritmy

Algoritmy této knihovny vycház´ı z pˇr´ımé implementace klasického algoritmu násoben´ı matic, samozˇrejmˇe s odliˇsným zp˚usobem adresován´ı (ˇr´ıdké matice obvykle nemaj´ı pˇr´ımý pˇr´ıstup k prvk˚um). Ve zdrojových kódech jsem nenalezl ˇ

zádné specifické triky a techniky slouˇz´ıc´ı k zlepˇsen´ı efektivity kromˇe rozdˇelen´ı funkc´ı dle parametr˚u na mnoho jednoduˇzˇs´ıch funkc´ı, m´ısto vyuˇzit´ı vˇetven´ı uvnitˇr jádra funkc´ı.

1.4.2 Eigen

Informace o této knihovnˇe pocház´ı z dokumentace dostupné na [16] a studia kódu dostupného z [17].

Eigen je open-source knihovna napsaná v C++, zabývaj´ıc´ı se lineárn´ı al- gebrou. Pro nás významnˇe se vˇenuje mimo jiné i ˇr´ıdkým matic´ım a operac´ım s nimi.

1.4.2.1 Interface

Interface této knihovny je mnohem bliˇzˇs´ı objektovému programován´ı (a modern´ım standard˚um C++). S maticemi se tedy pracuje pro objekty obvyklým zp˚usobem (konstruktory, destruktory, metody).

Z operac´ı s ˇr´ıdkými maticemi umoˇzˇnuje knihovna (mimo jiné) sˇc´ıtán´ı, násoben´ı skalárem, transpozici (a z´ıskán´ı hermitiánu matice pro komplexn´ı ˇ

c´ısla, vice o hermitianu napˇr´ıklad [18]), násoben´ı dvou matic a operace apli- kované po prvku (tedy výsledek z´ıskán jako C_i,j =op(A_i,j, B_i,j)).

Významnou zmˇenou oproti NIST Sparse BLAS je, ˇze výsledkem operace se dvˇema ˇr´ıdkými maticemi je znovu ˇr´ıdká matice, nikoliv matice hustá.

Pro ukládán´ıˇr´ıdkých matic je vyuˇz´ıván upravený formát CRS/CCS. Úpra- va spoˇc´ıvá v ponecháván´ı volných m´ıst pro pˇridáván´ı prvk˚u. Nav´ıc v nˇekterých

(32)

1. Anal´yza

m´ıstech vyuˇz´ıvá nekomprimovaného formátu, ve kterém jsou jednotlivé ˇrádky (pˇr´ıpadnˇe sloupce) ukládány do samostatných pol´ı (ve formátu odpov´ıdaj´ıc´ım jednomu ˇrádku/sloupci CRS/CCS) a na nˇe jsou uloˇzené odkazy (coˇz umoˇzˇnuje jednoduˇzˇs´ı a rychlejˇs´ı pˇridáván´ı prvk˚u do matice).

1.4.2.2 Algoritmy

I zde je vyuˇzit základn´ı algoritmus násoben´ı matic. Výsledná matice je vytvá- ˇrena sekvenˇcnˇe po sloupc´ıch. Celý algoritmus je významnˇe komplikován t´ım, ˇze výsledná matice je také v ˇr´ıkém formátu. Napˇr´ıklad je nutné odhadovat mnoˇzstv´ı nenulových prvk˚u ve výsledném sloupci (provádˇeno velmi jedno- duchým odhadem na základˇe pˇredpokladu, ˇze max 1 index nenulového prvku se liˇs´ı). Algoritmus vyˇzaduje mnoho dynamických alokac´ı (a realokac´ı) a velmi sloˇzitˇe paralelizovatelný.

1.4.3 cuSPARSE

Informace v této sekci vycház´ı z dokumentace dostupné na [19].

Knihovna cuSPARSE obsahuje sadu low-level funkc´ı pro práci s ˇr´ıdkými maticemi na GPU a je urˇcená pro jazyky C a C++.

Konkrétn´ı algoritmy pouˇz´ıvané v cuSPARSE nebudou rozeb´ırány, protoˇze kódy knihovny nejsou veˇrejnˇe k dispozici. Tato knihovna se vˇsak (kv˚uli po- dobné problematice, kterou se zabývá) dá vyuˇz´ıt pro porovnán´ı efektivity s kódy implementovanými bˇehem této diplomové práce.

1.4.3.1 Interface

Interface této knihovny je znaˇcnˇe podobný s NIST Sparse BLAS. V cuSPARSE se nacház´ı m´ırnˇe obecnˇejˇs´ı funkce typu

C ←αop(A)·op(B) +β·C,

(cusparse<t>csrmm2) kde op() znaˇc´ı identitu, transpozici, nebo hermiti´an.

Kromˇe matice A (a op(A)), kter´a je v ˇr´ıdkém formátu, jsou vˇsechny ostatn´ı matice husté aα, β jsou reálná ˇc´ısla.

Tato knihovna ovˇsem nav´ıc obsahuje i funkci n´asoben´ı dvou ˇr´ıdk´ych matic.

Konkr´etnˇe se jedn´a o funkce typu

C ← ·op(A)·op(B) (cusparse<t>csrgemm) a

C←α·op(A)·op(B) +β·D

(cusparse<t>csrgemm2). V tˇechto dvou pˇr´ıpadech jsou vˇsechny matice v CRS formátu. Tyto funkce vˇsak nejsou zcela implementovány, takˇzeop() m˚uˇze re- prezentovat pouze identitu (tedy bez transpozic) a nejsou podporována ˇzádná

(33)

1.4. Reˇserˇse existuj´ıc´ıch knihoven pokroˇcilejˇs´ı ukládán´ı matic, které knihovna v jiných ˇcástech umoˇzˇnuje (sy- metrické, trojúheln´ıkové aj.). Pˇred volán´ım tˇechto funkc´ı je nav´ıc nutné volat dalˇs´ı funkce pro výpoˇcet správné alokace maticeC a v druhém pˇr´ıpadˇe i alo- kaci bufferu (jehoˇz úˇcel nen´ı pˇresnˇe specifikován). Tyto dvˇe funkce budou dále pouˇzity k porovnán´ı efektivity s kódy vzniklými jako souˇcást této práce.

Matice jsou pˇredávány dle zp˚usobu obvyklého v jazyce C, tedy po jednot- livých pol´ıch a promˇenných. Toto vˇsak vede k znaˇcnˇe velkému mnoˇzstv´ı argu- ment˚u (napˇr´ıklad 27 pro csrgemm2), coˇz zvláˇstˇe v kombinaci s dˇr´ıve zm´ınˇeným v´ıcefázovým volán´ım funkc´ı nen´ı pˇr´ıliˇs uˇzivatelsky pˇr´ıvˇetivé.

(34)

(35)

Kapitola 2

N´ avrh

V prvn´ı fázi návrhu je nezbytné urˇcit, s jakými formáty matic bude implementace pracovat a jaký interface budou m´ıt metody pro práci s maticemi.

Dále je nutné navrhnout základn´ı sekvenˇcn´ı algoritmy pro práci s jednotlivými formáty a urˇcit, jak tyto algoritmy paralelizovat. Nakonec je také nezbytné vy- tvoˇrit základn´ı návrh tˇr´ıd.

2.1 Interface

Základ interface funkce násoben´ı (a sˇc´ıtán´ı) jsem se rozhodl ˇcerpat z knihovny NIST Sparse BLASS, ovˇsem za vyuˇzit´ı nˇekterých prvk˚u pouˇz´ıvaných knihovnou Eigen (a obecnˇe objektovým programován´ım). Nav´ıc jsem se vzhledem k charakteru prototypu a malým pˇr´ınos˚um opaku rozhodl podporovat pouze ukládán´ı hodnot ve formˇe ˇc´ısel s plovouc´ı ˇrádovou ˇcárkou s jednodu- chou pˇresnost´ı (tedy v C++ float).

Bude tedy prov´adˇena operace:

C←α·A·B+C, C←α·A^T ·B+C.

Argumenty operace budou konstantn´ı pointery na maticeA,B, hodnotaα (float) a pointer na výstupn´ı hustou maticiC, jehoˇz nulovost znamená, ˇze má být alokována nová matice. Dále bude volán´ı funkce obsahovat dvˇe výbˇerové hodnoty (enum), které urˇcuj´ı proveden´ı transpozice maticeA pˇred výpoˇctem a zda má být výpoˇcet proveden paralelnˇe.

2.2 Form´ aty matic

Z interface operace násoben´ı je zˇrejmé, ˇze mus´ı být pˇr´ıtomný formát pro husté matice.

Z formát˚u pro ukládán´ı ˇr´ıdkých matic byly dále vybrány:

(36)

2. N´avrh

• CRS - nejˇcastˇeji pouˇz´ıvaný, základn´ı formát ukládán´ı ˇr´ıdkýc matic

• BCRS - zaloˇzený na CRS, vhodný pro blokové matice

• DIA - prakticky nejefektivnˇejˇs´ı formát, má vˇsak velmi omezené pouˇzit´ı (pro matice s malým mnoˇzstv´ım diagonnál)

• ELL - nadstavba nad DIA (a CRS), obecnˇe povaˇzován za efektivn´ı a znaˇcnˇe obecný (pro matice s rovnomˇernˇe rozloˇzenými prvky na ˇrádky) Formát CCS nebyl specificky pˇridán, protoˇze je prakticky shodný s CRS transponované matice. Formát COO je pro samotné násoben´ı nevhodný a pa- mˇet’ovˇe nároˇcný, proto byl také z vnitˇrn´ı reprezentace vyˇrazen, aˇckoliv pˇrevody mezi n´ım a ostatn´ımi formáty se ukázaly být nezbytné pro z´ıskán´ı testovac´ıch matic (viz sekci testován´ı). Dynamické formáty byly vynechány z d˚uvodu n´ızké efektivity algoritm˚u násoben´ı (tedy opaku c´ıle této práce) a toho, ˇze tato práce nec´ıl´ı na snadné dynamické mˇenˇen´ı matic. Hybridn´ı formát nebyl pˇridán kv˚uli jeho sloˇzitosti a jeho pˇridán´ım by byl významnˇe pˇrekroˇcen zamýˇslený rozsah práce.

2.3 Sekvenˇ cn´ı algoritmy

C´ılem je implementovat (a paralelizovat) funkci α·A^(T⁾·B +C, kde A, B jsou dvˇe ˇr´ıdké matice uloˇzené ve stejném formátu aCje hustá matice, slouˇz´ıc´ı zároveˇn k uloˇzen´ı výsledku. K implementaci této funkce (a ovˇeˇren´ı funkˇcnosti) je zapotˇreb´ı 4 d´ılˇc´ıch funkc´ı. Tyto funkce jsou:

• Hust´a matice→ˇr´ıdk´a matice (vytvoˇren´ı)

• Transpozice

• Násoben´ı (s akumulac´ı do výsledkové matice)

• R´ıdk´ˇ a matice→ hustá matice (pro kontroln´ı výsledky a testován´ı) Z nich posledn´ı je obvykle velmi jednoduchá (a zˇrejmá), proto zde nebude specificky popisována. Jediné násoben´ı pak má (asymptoticky) vyˇsˇs´ı výpoˇcetn´ı sloˇzitost neˇz sloˇzitost pamˇet’ovou, a proto je jej vhodné paralelizovat. Vˇsechny ostatn´ı funkce maj´ı stejnou pamˇet’ovou a výpoˇcetn´ı sloˇzitost (O(nnz)), a proto jejich paralelizace pomoc´ı GPU (kv˚uli pˇrenos˚um dat) nepˇrinese ˇzádné zrych- len´ı.

Abych pˇredeˇsel exponenciáln´ımu nár˚ustu poˇctu nutných impelementac´ı násoben´ı, rozhodl jsem se implementovat pouze násoben´ı matic ve stejném formátu (u BCRS nav´ıc se stejnou velikost´ı bloku).

(37)

2.3. Sekvenˇcn´ı algoritmy 2.3.1 CRS

2.3.1.1 Datov´a struktura

Z definice formátu vyplývá, ˇze pro uloˇzen´ı matice typun×msnnznenulovými elementy budou potˇreba nejménˇe 3 pole. Pole offset˚u (int) o délce n+ 1 (na pozici nje uloˇzenonnz), pole index˚u sloupc˚u (int) o délcennz a pole hodnot (float) také o délcennz. Aby bylo moˇzné reprodukovat p˚uvodn´ı matici, je dále nezbytné m´ıt uloˇzenou nejen hodnotu n, ale im.

2.3.1.2 Pˇrevod z hust´e matice

Jak m˚uˇzete vidˇet na algoritmu 2, pˇrevod je pomˇernˇe intuitivn´ı (algoritmus pˇresto pˇrepsán, protoˇze se bude hodit pˇri popisu dalˇs´ıch formát˚u). Za zm´ınku stoj´ı pˇredevˇs´ım ˇrádek 4. Tuto operaci lze odebrat a vyuˇz´ıt hodnoty count z hlavn´ıho for cyklu. Problém tohoto pˇr´ıstupu vˇsak je, ˇze bychom museli vyuˇz´ıvat dynamické alokace pol´ı (taktéˇz O(n· m), ovˇsem pravdˇepodobnˇe s vyˇsˇs´ımi konstantami), takˇze výsledný kód by byl pomalejˇs´ı a zároveˇn sloˇzi- tˇejˇs´ı.

Algoritmus 2 Hust´a matice →CRS

1: procedure CRS(M :DenseM atrix)

2: C:CRS .Tvoˇren´a matice

3: C.height←M.height;C.width←M.width

4: nnz←count nnz(M) . O(n·m)

5: C.offsets←zeroes(C.height)

6: C.indexes, C.values←zeroes(nnz)

7: count←0

8: fori←0 to M.heightdo

9: C.offsets[i]←count

10: for j←0 to M.width do

11: if notEmpty(M[i][j])then

12: C.values[count]←M[i][j]

13: C.indexes[count+ +]←j

14: C.offsets[M.height]←nnz

15: return C

2.3.1.3 Transpozice

Je zˇrejmˇe neˇz´adouc´ı prov´adˇet tuto transpozici pˇres hustou matici (O(n·m)).

Jako druhá moˇznost se jev´ı pˇreˇrazen´ı prvk˚u. Obvyklé ˇrad´ıc´ı algoritmy maj´ı sloˇzitotO(nnz∗lognnz), coˇz je významnˇe lepˇs´ı, avˇsak ne ideáln´ı. M˚uˇzeme se vˇsak inspirovat z counting sortu. Je podstatné, ˇze hodnoty jsou jiˇz seˇrazené dle ˇrádk˚u. Staˇc´ı tedy, abychom je counting sortem seˇradili podle sloupc˚u, a d´ıky

(38)

2. N´avrh

stabilnosti tohoto algoritmu z´ıskáme poˇzadované seˇrazen´ı, a tedy i transpono- vanou matici. Pole offset˚u z´ıskáme pˇr´ımo pˇri provádˇen´ı counting sortu (pre- fixový souˇcet poˇct˚u jednotlivých hodnot), jen jej oproti standardn´ı variantˇe mus´ıme posunout o jedno pole doprava, tedy tak, aby

∀x∈N∪ {0}, x≤m:offset_x =

x−1

X

k=0

nnz_k,

kdennz_k je poˇcet nenulov´ych prvk˚u ve sloupci k.

Sloˇzitost celé transpozice je (ze sloˇzitosti counting sortu)O(nnz+m) a za pˇredpokladu, ˇze nnz ≥ m (tedy, ˇze v kaˇzdém sloupci je pr˚umˇernˇe alespoˇn jeden prvek) je tato sloˇzitostO(nnz), coˇz je optimáln´ı sloˇzitost, protoˇze celou matici mus´ıme pˇrepsat.

2.3.1.4 N´asoben´ı

Algoritmus násoben´ı vycház´ı z klasického postupu násoben´ı (viz algoritmus 1, dále indexy znaˇceny dle tohoto algoritmu). D˚uleˇzité vˇsak je, v jakém poˇrad´ı budeme matice procházet. Formát CRS je vhodný pouze k procházen´ı po ˇrádc´ıch. M˚uˇzeme vˇsak vyuˇz´ıt vlastnosti, ˇze transpozic´ı uloˇzené matice z´ıskáme stejnou matici uloˇzenou ve formátu CCS (zˇrejmé z definice CRS/CCS), který naopak dobˇre pracuje po sloupc´ıch. Dále budeme zp˚usob procházen´ı násoben´ı maticA·B znaˇcit dle uˇzitého formátu, tedy napˇr´ıklad CRS·CCS.

Z definice maticov´eho souˇctu jsou zˇrejm´e tyto dva fakty:

• Mezi sebou n´asob´ıme vˇsechny prvkya_i,k s prvkyb_k,j, tedy k-t´y sloupec z maticeA s k-t´ym ˇr´adkem z maticeB

• Hodnotu prvku cij matice C z´ıskáme jako skalárn´ı souˇcin i-t´eho ˇrádku maticeA s j-t´ym sloupcem maticeB

Tyto dva pohledy nás vedou na dva logické zp˚usoby násoben´ı, coˇz jsou CCS·CRS (dále také CR) a CRS·CCS (dále také RC). Výhodou prvn´ıho je, ˇze nemus´ıme porovnávat hodnoty index˚u a jediná provádˇená operace je tedy samotné násoben´ı. Tento zp˚usob se tud´ıˇz zˇrejmˇe hod´ı pro sériovou verzi a vyˇzaduje nejménˇe provedených operac´ı. Druhý zp˚usob vyˇzaduje porovnáván´ı index˚u a je významnˇe nároˇcnˇejˇs´ı na výpoˇcetn´ı prostˇredky, na druhou stranu vˇsak umoˇzˇnuje samostatný výpoˇcet hodnot jednotlivých prvk˚u matice. Tato vlastnost umoˇzˇnuje CRS ·CCS masivnˇe paralelizovat, aniˇz bychom museli vyuˇz´ıt nákladnou synchronizaci zápis˚u do výsledné matice (která by pohltila témˇeˇr vˇsechny výhody paralelizace), nebo zápis do v´ıce pomocných matic a následnou redukci.

(39)

2.3. Sekvenˇcn´ı algoritmy 2.3.2 BCRS

Datová struktura BCRS je totoˇzná s CRS. Rozd´ılem je, ˇze jsou zde indexy blo- kových sloupc˚u (nikoliv samotných prvk˚u) a bude m´ıt tedy pouze délkunnzb, coˇz je poˇcet nenulových blok˚u. Nav´ıc je tˇreba uloˇzit velikost bloku. Pro jednoduchost (a malé výhody opaku) jsem se rozhodl podporovat pouze ˇctvercové bloky, takˇze tato velikost m˚uˇze být uloˇzena za pomoc´ı jediné hodnoty bs.

Algoritmus 3 Hust´a matice →BCRS

1: procedure CRS(M :DenseM atrix, bs:int)

2: C:BCRS .Tvoˇren´a matice

3: C.height←M.height;C.width←M.width;bs2←bs·bs

4: max len←C.width∗C.height/bs2

5: C.offsets←zeroes(C.height)

6: C.indexes←zeroes(max len);C.values←zeroes(max_len∗bs2)

7: count←0

8: fori←0 to M.height/blokdo

9: C.offsets[i]←count

10: for j←0 to M.width do

11: if notEmpty(M[i∗bs: (i+ 1)∗bs][j∗bs: (j+ 1)∗bs])then

12: C.values[count∗bs2 : (count+ 1)∗bs2]←

13: M[i∗bs: (i+ 1)∗bs][j∗bs: (j+ 1)∗bs]

14: C.indexes[count+ +]←j

15: C.offsets[M.height]←nnz

16: return C

Jak si m˚uˇzeme vˇsimnout na algoritmu 3, lze tento pˇrevod provést prakticky stejnˇe jako u CRS (viz algoritmus 2). Rozd´ıly jsou ve velikosti alokovaných pol´ı (po dokonˇcen´ı pˇrevodu mohou být pole zmenˇsena), a dále pak na ˇrádc´ıch 11, kde m´ısto testován´ı jednoho prvku testujeme celý blok a na ˇrádku 12, kde kop´ırujeme celý blok hodnot.

2.3.2.3 Transpozice

I zde m˚uˇzeme vyuˇz´ıt jiˇz hotovou transpozici z CRS. Celý postup m˚uˇze pro- bˇehnout stejnˇe, m´ısto prvk˚u vˇsak pracujeme s celými bloky prvk˚u. Samotné bloky mus´ı být také transponovány zp˚usobem obvyklým pro husté matice, coˇz lze provést napˇr´ıklad pˇri jejich pˇresunu v counting sortu.

(40)

2. N´avrh

2.3.2.4 N´asoben´ı

Vynásoben´ım dvou blok˚u vˇzdy vznikne nový blok hodnot (který m˚uˇze být pˇr´ımo ukládán do husté matice, pˇr´ıpadnˇe ukládán bˇehem výpoˇctu po hod- notách). I zde tedy lze násobit dvˇema zp˚usoby, tedy BCCS·BCRS (CR) aBCRS·BCCS(RC). Prvky uvnitˇr blok˚u vˇsak nemá velký smysl násobit jinak, neˇz pomoc´ı nezávislých skalárn´ıch souˇct˚u (ˇrádek×sloupec) a postupného ukládán´ı, protoˇze bloky samy o sobˇe jsou husté, a tedy nevzniká dodateˇcná nároˇcnost s porovnáván´ım index˚u.

2.3.3 ELL

Pro uloˇzen´ı matice (n×m) v ELL form´atu je nezbytné uchovat dvˇe matice velikosti n×ell width, tedy matici hodnot a matici sloupcov´ych index˚u, kdeell widthje nejvˇetˇs´ı poˇcet nenulových prvk˚u na ˇrádek. Pˇrebyteˇcné m´ısto v ˇrádc´ıch s ménˇe nenulovými prvky m˚uˇze být pro jednoduchost práce s matic´ı vyplnˇeno nulami. Dále je samozˇrejmˇe nezbytné také uloˇzit ell width a ˇs´ıˇrku a výˇsku p˚uvodn´ı matice.

Algoritmus 4 Hust´a matice→ ELL

1: procedureCRS(M :DenseM atrix)

2: C :ELL . Tvoˇren´a matice

3: C.height←M.height;C.width←M.width

4: C.ell width←max(count nnz(M.rows)) . O(n·m)

5: C.indexes, C.values←zeroes(C.ell width∗C.height)

6: fori←0 to M.heightdo

7: row index←0

8: forj←0 to M.widthdo

9: if notEmpty(M[i][j]) then

10: C.values[i∗C.ell width+row index]←M[i][j]

11: C.indexes[i∗C.ell width+row index]←j;

12: row index+ +;

13: return C

Aˇckoliv je struktura tohoto formátu znaˇcnˇe rozd´ılná od CRS, i zde lze pˇrevod provést podobným zp˚usobem. Pseudokód tohoto algoritmu si m˚uˇzete prohlédnout na algoritmu 4.