Objektově orientovaná implementace skyline systému

ˇ – Technicka´ univerzita Ostrava VSB Fakulta elektrotechniky a informatiky Katedra informatiky

Objektoveˇ orientovana´ ´ implementace skyline systemu ´ an´ ´ ı rˇ´ıdkych uklad ´ matic ´ rska´ prace ´ Bakalaˇ

2008

Martin Stachonˇ

´ ı Zadan´ Téma: Diplomant: Vedouc´ı: Akad. rok: Obor:

Objektovˇe orientovaná implementace skyline systému ukládán´ı rˇ´ıdkych ´ matic Martin Stachonˇ Doc. Mgr. V´ıt Vondrák, Ph.D. 2007 / 2008 1103R021 Poˇc´ıtaˇcová matematika

´ ´ ı Zasady pro vypracovan´ ´ V nároˇcnych inˇzenyrsk uloh´ ach se velmi cˇ asto objevuj´ı rozsáhlé systémy soustav ´ ´ ych ´ ˚ Tyto matice lineárn´ıch rovnic jejichˇz matice obsahuj´ı velké mnoˇzstv´ı nulovych prvku. ´ ´ se nazyvaj´ bude ´ ı rˇ´ıdké a je pochopitelné, zˇ e kl´ıcˇ ovou roli v efektivnosti rˇ eˇsen´ı tˇechto uloh hrát rychlost pˇri práci s tˇemito maticemi. Dalˇs´ı neménˇe vyznamnou roli pak hraje i ve´ likost pamˇeti, kterou tato struktura pˇri bˇezˇ nych ´ maticovych ´ operac´ıch zab´ırá. Jednou z velmi pouˇz´ıvanych ´ struktur je tzv. skyline systém ukládán´ı matic, ktery´ spoleˇcnˇe se Sloanovym ´ algoritmem pro pˇreˇc´ıslován´ı matic je velmi hojnˇe pouˇz´ıván pˇredevˇs´ım tam, kde je zapotˇreb´ı rychlá LU cˇ i Choleského faktorizace. Práce se dá rozloˇzit do následuj´ıc´ıch ˚ bodu: ˚ algoritmus pro sn´ızˇ en´ı profilu matic • Sloanuv • Objektovˇe orientovany´ návrh implementace skyline systému ukládán´ı matic • Implementace vyˇ ´ se uvedeného v jazyce C++

Prohlaˇsuji, zˇ e jsem tuto bakalárˇ skou práci vypracoval samostatnˇe. Uvedl jsem vˇsechny literárn´ı prameny a publikace, ze kterych ´ jsem cˇ erpal.

V Ostravˇe 2. kvˇetna 2008

.............................

Zde bych rád podˇekoval Ing. Pavle Kabel´ıkové a Doc. Mgr. V´ıtu Vondrákovi, Ph.D. za poskytnut´ı literatury, konzultac´ı a pˇr´ıkladu˚ pro testován´ı. Pˇri tvorbˇe této práce jsem mnohokrát pouˇzil GNU Octave, dˇekuji vˇsem, kteˇr´ı se pod´ılej´ı na vyvoji tohoto volnˇe ´ sˇ iˇritelného softwaru.

Abstrakt Tato bakalárˇ ská práce se zabyv´ ´ a Skyline systémem ukládán´ı rˇ´ıdkych ´ matic a jeho objek˚ tovˇe orientovanou implementac´ı v programovac´ım jazyce C++. Dále popisuje Sloanuv ´ celem sn´ızˇ en´ı profilu a jeho implementaci v proalgoritmus pro pˇreˇc´ıslován´ı matice za uˇ gramovac´ım jazyce C++. ˇ a´ slova: Sloanuv ˚ algoritmus, rˇ´ıdké, matice, skyline, C++ Kl´ıcov

Abstract This thesis describes the Skyline system of sparse matrix storage and its object-oriented implementation in C++ programming language. Also, Sloan’s algorithm for matrix reordering for lower matrix profile is described and implemented in C++. Keywords: Sloan’s algorithm, sparse, matrix, skyline, C++

Seznam pouˇzitych ´ zkratek a symbolu˚ OOSol

–

RCM

–

Object Oriented SOLvers, objektovˇe orientovaná knihovna v programovac´ım jazyce C++ pro numerické vypoˇ ´ cty Reverse Cuthill–McKee, algoritmus pro pˇreˇc´ıslován´ı vrcholu˚

1

Obsah 1

´ Uvod

2

Skyline matice 2.1 Základn´ı pojmy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ˇ ıdké matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 R´ 2.3 Skyline matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7 7 9 10

3

Systém ukládán´ı Skyline matic 3.1 Datové struktury . . . . . . . . . . . . . . . . . . . . . 3.2 Pamˇet’ová nároˇcnost . . . . . . . . . . . . . . . . . . . ˚ 3.3 Pˇr´ıstup k prvkum . . . . . . . . . . . . . . . . . . . . . 3.4 Sestaven´ı skyline systému . . . . . . . . . . . . . . . . 3.5 Implementace skyline systému pro knihovnu OOSol .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

12 12 13 13 13 14

Grafy 4.1 Základn´ı pojmy z teorie grafu˚ . . . . . . . . 4.2 Souvislost mezi grafem a strukturou matice 4.3 Ukládán´ı grafu pomoc´ı seznamu sousedu˚ . 4.4 Sestaven´ı seznamu sousedu˚ . . . . . . . . . ´ ˇ a struktura . . . . . . . . 4.5 Koˇrenová urov nov´

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

18 18 19 20 20 20

Sloanuv ˚ algoritmus ˚ eru . . . . . . . . . . . . . . . . . . . 5.1 Urˇcen´ı pseudoprumˇ 5.2 Pˇreˇc´ıslován´ı vrcholu˚ . . . . . . . . . . . . . . . . . . . . 5.3 Hodnoty parametru˚ W1 a W2 . . . . . . . . . . . . . . . 5.4 Implementace sloanova algoritmu pro knihovnu OOSol

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

22 22 22 24 25

4

5

6

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

6

Testován´ı a vysledky ´

28

7

Závˇer

33

8

Reference

34

Pˇr´ılohy

34

A Zdrojové kody ´

35

2

Seznam tabulek 1 2

Testovac´ı matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Vysledn´ e profily . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ´

31 31

3

´ u˚ Seznam obrazk 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

ˇ e kˇr´ızˇ ky reprezentuj´ı nenulovy´ prvek. Vystup pˇr´ıkazu spy. Sed´ ´ Ukázka matice se skyline profilem. . . . . . . . . . . . . . . . . . Ukázka matice se symetrickym ´ skyline profilem. . . . . . . . . . Hierarchie tˇr´ıd reprezentuj´ıc´ıch matice v knihovnˇe OOSol . . . Pˇr´ıklad grafu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Stavy vrcholu˚ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S´ıt’ matice ctverec1 . . . . . . . . . . . . . . . . . . . . . . . . . S´ıt’ matice ctverec3 . . . . . . . . . . . . . . . . . . . . . . . . . S´ıt’ matice krychle1 . . . . . . . . . . . . . . . . . . . . . . . . . S´ıt’ matice krychle2 . . . . . . . . . . . . . . . . . . . . . . . . . S´ıt’ matice krychle3 . . . . . . . . . . . . . . . . . . . . . . . . . S´ıt’ matice krychle4 . . . . . . . . . . . . . . . . . . . . . . . . . S´ıt’ matice UF871 - helikoptéra Commanche . . . . . . . . . . . . ˚ Puvodn´ ı profil matice krychle4 . . . . . . . . . . . . . . . . . . Profil matice krychle4 po pˇreˇc´ıslován´ı sloanovym ´ algoritmem Profil matice krychle4 po pˇreˇc´ıslován´ı algoritmem RCM . . . ˚ Puvodn´ ı profil matice ctverec3 . . . . . . . . . . . . . . . . . . Profil matice ctverec3 po pˇreˇc´ıslován´ı sloanovym ´ algoritmem

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

10 10 11 17 19 24 29 29 29 29 29 30 30 30 32 32 32 32

4

Seznam algoritmu˚ 1 2 3 4

˚ matice ve skyline systému Pˇr´ıstup k prvkum ´ ˇ e struktury . . . Sestaven´ı koˇrenové urov nov´ ˚ eru grafu . . . . . . . . Nalezen´ı pseudoprumˇ ˚ algoritmus pro pˇreˇc´ıslován´ı vrcholu˚ Sloanuv

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

14 21 23 27

5

´ ´ Seznam vypis ´ u˚ zdrojoveho kodu 1 2 3 4

Soubor OOSmatricesSkylineStatic.h . Soubor OOSmatricesSkylineStatic.cpp Soubor OOSgraph.cpp . . . . . . . . . Soubor OOSgraph.h . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

36 39 58 65

6

1

´ Uvod

Ve vˇedecké a technické praxi se, napˇr´ıklad pˇri rˇ eˇsen´ı soustav parciálnˇe diferenciáln´ıch ˚ cˇ asto objevuj´ı rˇ´ıdké matice, tedy matice s mnoha rovnic metodou koneˇcnych prvku, ´ nulovymi prvky. Pro efektivn´ı uloˇzen´ı tˇechto rozsáhlych matic a manipulaci s nimi je ´ ´ vyhodn´ e a u rozsáhlejˇs´ıch problému˚ nutné, pouˇz´ıvat efektivn´ı datové struktury. V této ´ bakalárˇ ské práci pop´ısˇ i jednu z tˇechto struktur, takzvany´ Skyline storage, ktery´ je vyhodn y´ ´ pro rˇ´ıdké matice s nenulovymi prvky soustˇredˇenymi okolo diagonály. ´ ´ Jelikoˇz ne kaˇzdá rˇ´ıdká matice má uspoˇra´ dán´ı vhodné pro Skyline storage, pouˇz´ıvaj´ı se pˇreˇc´ıslovac´ı algoritmy, které pˇreuspoˇra´ daj´ı rˇ a´ dky a sloupce matice tak, aby se sn´ızˇ il jejich profil a t´ım pádem i pamˇet’ potˇrebná pro jejich uloˇzen´ı. Jedn´ım z nejefektivnˇejˇs´ıch je algoritmus, jehoˇz autorem je S. W. Sloan a tento algoritmus zde pop´ısˇ i. Souˇca´ st´ı této práce je i implementace vyˇ ´ se uvedeného v jazyce C++ pro knihovnu OOSol.

7

2 2.1

Skyline matice ´ Zakladn´ ı pojmy

Definujme nˇekolik základn´ıch pojmu˚ tykaj´ ´ ıc´ıch se matic, které pak budou v této práci dále pouˇzity. Definice 2.1 Matic´ı rozum´ıme tabulku matematických výrazu˚ o m rˇa´ dc´ıch a n sloupc´ıch.   a11 a12 . . . a1n  a21 a22 . . . a2n    A= . ..  ..  .. . .  am1 am2 . . . amn aij pak oznaˇcujeme jako prvek matice A na i-tém rˇa´ dku a j-tém sloupci. V této prác´ı se budeme zabyvat pouze reálnými maticemi, to jest takovymi, jejichˇz ´ ´ prvky jsou pouze reálná cˇ ´ısla, aij ∈ R Poznámka 2.1 V matematice byv´ ´ a zvykem pouˇz´ıvat indexy zaˇc´ınaj´ıc´ı od jedniˇcky. V po˚ zpusobu ˚ pisu C++ programu se kvuli indexace pol´ı v tomto jazyce pouˇz´ıvaj´ı indexy zaˇc´ınaj´ıc´ı od nuly. ˇ Definice 2.2 Ctvercov a´ matice je taková, která má stejný poˇcet rˇa´ dku˚ jako sloupcu, ˚ m=n Definice 2.3 Symetricka´ matice je matice, která je cˇ tvercová a nav´ıc plat´ı aij = aji pro vˇsechna 1 ≤ i ≤ n, 1 ≤ j ≤ n ´ Definice 2.4 i-ty´ rˇadek rA (i) matice A je vektor rA (i) =

ai1 , ai2 , . . . , ain

Definice 2.5 j-ty´ sloupec sA (i) matice A je vektor  a1j  a2j  sA (i) =  .  ..

    

amj ´ ı prvek aij matice A je takový prvek, pro který plat´ı i = j. Definice 2.6 Diagonaln´ ´ matice A je vektor vˇsech diagonáln´ıch prvku˚ Definice 2.7 Diagonala diag(A) = a11 , a22 , . . . , akk kde k = min{m, n}

8

´ ´ jsou prvky nad diagonálou, tedy takové aij , pro které Definice 2.8 Horn´ı trojuheln´ ıkova´ cˇ ast plat´ı i < j. ´ ´ jsou prvky pod diagonálou, tedy takové aij , pro které Definice 2.9 Doln´ı trojuheln´ ıkova´ cˇ ast plat´ı i > j. Pˇr´ıklad 2.1 Ukaˇzme si tyto pojmy na konkrétn´ı matici. 

 1 2 3 A= 4 5 6  7 8 9

Tato matice má 3 sloupce a 3 rˇ a´ dky, je cˇ tvercová. Prvek v druhém rˇ a´ dku, tˇret´ım sloupci: a2,3 = 6 Tˇret´ı rˇ a´ dek matice : rA (3) = 7 8 9   1 A  4  Prvn´ı sloupec matice : s (1) = 7 Diagonála matice : diag(A) = 1 5 9   2 3 ´ 6  Horn´ı trojuheln´ ıková cˇ a´ st U =   ´ Doln´ı trojuheln´ ıková cˇ a´ st L =  4 7 8

 

Definice 2.10 Symetrickou permutac´ı cˇ tvercové matice o rozmˇeru m permutaˇcn´ım vektorem p délky m, nazveme zobrazen´ı A˜ = symperm(A, p), takové zˇ e, a ˜i,j = api ,pj ˚ Napˇr´ıklad pro Poznámka 2.2 Permutaˇcn´ı vektor p udává nové poˇrad´ı rˇ a´ dku˚ a sloupcu. matici A z pˇredchoz´ıch pˇr´ıkladu˚ bude symetrická permutace vektorem p = (3, 1, 2) :   9 7 8 A˜ =  3 1 2  6 4 5 Nˇekdy se m´ısto permutaˇcn´ıho vektoru pouˇz´ıvá takzvany´ inverzn´ı permutaˇcn´ı vektor, ktery´ m´ısto poˇrad´ı sloupcu˚ a rˇ a´ dku˚ udává na pozici i novy´ index pro i-ty´ sloupec nebo rˇ a´ dek. Inverzn´ı permutaˇcn´ı vektor pro vektor p je p¯ = (2, 3, 1). Permutaci lze také zavést pomoc´ı permutaˇcn´ı matice P , a permutaci pak lze zapisovat jako násoben´ı matic. Tento pˇr´ıstup se ale v praktické implementaci nepouˇz´ıvá.

9

2.2

ˇ ıdke´ matice R´

Definice 2.11 Jako nenulovy´ prvek oznaˇcujeme takový prvek matice aij , pro který plat´ı aij 6= 0, Naproti tomu Definice 2.12 Jako nulovy´ prvek oznaˇcujeme takový prvek matice aij , pro který plat´ı aij = 0, Definice 2.13 Necht’ je dána matice A o m rˇa´ dc´ıch a n sloupc´ıch. Pak definujeme poˇcet nenulových prvku˚ nnz(A) = ||{aij : 1 ≤ i ≤ m, 1 ≤ j ≤ n, aij 6= 0}|| Poˇcet nulovych ´ prvku je potom roven m · n − nnz(A) ˚ zeme definovat rˇ´ıdkou matici jako matici, která má mnohem ménˇe nenuNyn´ı jiˇz muˇ lovych ´ prvku˚ neˇz je celkovy´ poˇcet prvku˚ v matici. Definice 2.14 Matice A o m rˇa´ dc´ıch a n sloupc´ıch se nazývá rˇı´dka´ matice, jestliˇze nnz(A) m·n Pˇr´ıklad 2.2 Pˇr´ıklad rˇ´ıdké matice s tˇremi nenulovymi prvky ´    S=  

0 0 2 0 0

0 0 0 0 0

0 0 0 0 0

0 1 0 0 3

0 0 0 0 0

    , nnz(S) = 3  

Pˇri zápisu v tabulce se pro pˇrehlednost nuly na m´ıstˇe nulovych ´ prvku nevypisuj´ı. Tak tomu bude nadále i v této práci. ˇ ıdké matice muˇ ˚ zeme vizuálnˇe interpretovat tak, zˇ e matici zobraz´ıme jako obrázek R´ o m · n bodech, pˇriˇcemˇz kaˇzdy´ bod reprezentuje jeden prvek matice a barevnˇe odliˇs´ıme nulové a nenulové prvky. V prostˇred´ı Matlab nebo Octave k tomu lze pouˇz´ıt pˇr´ıkaz spy(A). Pˇr´ıklad takového zobrazen´ı je na obrázku 1. Poznámka 2.3 Pˇri zápisu rˇ´ıdkych ´ matic se, zejména v matematickém softwaru pouˇz´ıvá zápis pomoc´ı tˇr´ı vektoru˚ i, j, c. Vˇsechny tˇri vektory maj´ı stejny´ poˇcet sloˇzek jako poˇcet nenulovych ´ prvku˚ matice. k-ty´ nenulovy´ prvek v matici je pak na rˇ a´ dku ik , sloupci jk a má hodnotu vk . Ostatn´ı prvky matice jsou implicitnˇe nulové. Napˇr´ıklad pro vyˇ ´ se uvedenou matici S : i = (2, 3, 5), j = (4, 1, 4), v = (1, 2, 3) ˚ matice s2,4 = 1, s3,1 = 2, s5,4 = 3 Coˇz odpov´ıdá nenulovym ´ prvkum

10

ˇ e kˇr´ızˇ ky reprezentuj´ı nenulovy´ prvek. Obrázek 1: Vystup pˇr´ıkazu spy. Sed´ ´       S=     

2 1 2 1 1

1 1 2 1 1 1 0

 1 2 0 0 1

 1   1   1 1  2 0 1   0 2 0 1   0 0 2 0  2

Obrázek 2: Ukázka matice se skyline profilem.

2.3

Skyline matice

V této práci se budu dále zabyvat speciáln´ım druhem rˇ´ıdkych matic, takzvané Sky´ ´ line matice. Tyto se vyznaˇcuj´ı t´ım, zˇ e nenulové prvky jsou soustˇredˇeny okolo diagonály. Název je odvozen z toho, zˇ e prvky nad diagonálou pˇripom´ınaj´ı siluetu mˇesta s mrakodrapy - anglicky skyline. Ukázka takové matice je na obrázku 2 . Pro pˇresnˇejˇs´ı popis skyline matice si definujme pojmy : profil rˇ a´ dku, profil sloupce, ˚ vektor sloupcovych ˚ rˇ a´ dkovy´ profil, sloupcovy´ profil, vektor rˇ a´ dkovych ´ profilu, ´ profilu, ˚ zité pro urˇcen´ı profil matice a matice se symetrickym ´ profilem. Tyto pojmy budou duleˇ poˇzadované velikosti pamˇeti, která je potˇrebná pro uloˇzen´ı matice ve skyline systému. Definice 2.15 Profil j-tého sloupce sp (j) matice A rozmˇerech (m, n) definujeme jako sp (j) = j − min{i : aij 6= 0, 1 ≤ i ≤ j} ˚ zeme vyjádˇrit jako poˇcet prvku˚ v j-tém sloupci nad Slovnˇe profil j-tého sloupce muˇ ´ diagonáln´ım prvkem (v horn´ı trojuheln´ ıkové cˇ a´ sti) aˇz do prvn´ıho nenulového (s nejniˇzsˇ´ım rˇ a´ dkovym ´ indexem), vˇcetnˇe. Obdobná je definice rˇ a´ dkového profilu, která odpov´ıdá

11



2 1  3 2 1   1 2 T =  1 1   1

 7 1 2 0 1

   1  1 3   2 3  0 2

Obrázek 3: Ukázka matice se symetrickym ´ skyline profilem. ´ poˇctu prvku˚ v i-tém rˇ a´ dku nalevo od diagonáln´ıho prvku v doln´ı trojuheln´ ıkové cˇ a´ sti aˇz po nenulovy´ prvek s nejniˇzsˇ´ım sloupcovym ´ indexem. Definice 2.16 Profil i-tého rˇa´ dku rp (i) matice A rozmˇerech (m, n) definujeme jako rp (i) = i − min{j : aij 6= 0, 1 ≤ j ≤ i} ˚ Vektor rˇ a´ dkovych ´ profilu˚ je pak vektor, jehoˇz sloˇzky jsou profily jednotlivych ´ rˇ a´ dku, ˚ vektor sloupcovych ´ profilu˚ pak obsahuje profily sloupcu. Definice 2.17 Vektor rˇa´ dkových profilu˚ rh matice o rozmˇerech (m,n) je vektor délky m, rhi = rp (i), 1 ≤ i ≤ m. Vektor sloupcových profilu˚ sh matice je vektor délky n, shi = sp (i), 1 ≤ i ≤ n ˇ adkový profil matice o rozmˇerech (m,n) je roven souˇctu profilu˚ vˇsech rˇa´ dku˚ maDefinice 2.18 R´ m n P P tice. rs = rp (i). Sloupcový profil je roven souˇctu profilu˚ vˇsech sloupcu. ˚ ss = sp (j) i=1

j=1

Definice 2.19 Matice má symetrický profil, pokud jsou si rovny vektory rˇa´ dkových a sloupcových profilu, ˚ rh = sh . Definice 2.20 Profil p matice je roven souˇctu rˇa´ dkového a sloupcového profilu p = rs + ss . Definice 2.21 Prvek aij leˇz´ı uvnitˇr profilu matice, pokud : • i=j • i < j (horn´ı trojuheln´ ´ ıková cˇ a´ st) : i ≥ j − sp (j) • i > j (doln´ı trojuheln´ ´ ıková cˇ a´ st) : j ≥ i − rp (i) Pˇr´ıklad 2.3 Vektory profilu˚ pro matici S na obrázku 2 jsou, rˇ a´ dkovy´ : rh = (0, 0, 1, 2, 2, 4, 3, 0) a sloupˇ adkovy´ profil : rs = 12, sloupcovy´ profil : ss = 13. Celcovy´ sh = (0, 1, 2, 1, 3, 0, 3, 3). R´ kovy´ profil matice p = 12 + 13 = 25. Matice nemá symetricky´ profil. Pˇr´ıklad matice se symetrickym ´ profilem je na obrázku 3.

12

´ uklad ´ an´ ´ ı Skyline matic System

3

Skyline systém vyuˇz´ıvá profilu matice, a ukládá pouze prvky na diagonále, prvky v pro´ filu kaˇzdého sloupce v horn´ı trojuheln´ ıkové cˇ a´ sti a prvky v profilu kaˇzdého rˇ a´ dku doln´ı ´ trojuheln´ ıkové cˇ a´ sti. Nulové prvky matice mimo profil se neukládáj´ı, ale uvnitˇr profilu ano. Pro matici s rozmˇery (m,n) jsou pouˇzity tyto datové struktury :

3.1

Datove´ struktury

• Pole celych ´ cˇ ´ısel diag délky min{m, n} slouˇz´ı k ukládán´ı prvku˚ na diagonále. Diagonáln´ı prvek aii je na i-té pozici v poli. • Celoˇc´ıselná promˇenná luc je rovna sloupcovému profilu ss matice. • Pole reálnych ´ cˇ ´ısel uc délky luc obsahuje postupnˇe (prvky jsou rˇ azeny s rostouc´ım rˇ a´ dkovym ´ indexem) po sloupc´ıch vˇsechny prvky uvnitˇr sloupcového profilu. • Pole celych ´ cˇ ´ısel iuc délky n + 1 obsahuje indexy zaˇca´ tku˚ sloupcu˚ v poli uc. Prvky ´ v profilu j-tého sloupce horn´ı trojuheln´ ıkové cˇ a´ sti matice tedy nalezneme v poli uc na indexech iuc[j] aˇz iuc[j+1]-1. Profil j-tého sloupce sp (j) = iuc[j + 1] − iuc[j] • Celoˇc´ıselná promˇenná llr je rovna rˇ a´ dkovému profilu rs matice. • Pole reálnych ´ cˇ ´ısel lr délky llr obsahuje postupnˇe (prvky jsou rˇ azeny s rostouc´ım sloupcovym ´ indexem) po rˇ a´ dc´ıch vˇsechny prvky uvnitˇr rˇ a´ dkového profilu. • Pole celych ´ cˇ ´ısel ilr délky m + 1 obsahuje indexy zaˇca´ tku˚ rˇ a´ dku˚ v poli lr. Prvky v ´ profilu i-tého rˇ a´ dku doln´ı trojuheln´ ıkové cˇ a´ sti matice tedy nalezneme v poli lr na indexech ilr[i] aˇz ilr[i+1]-1. Profil i-tého rˇ a´ dku rp (i) = ilr[i + 1] − ilr[i] Pˇr´ıklad 3.1 ˇ Ukaˇzme si obsah tˇechto datovych zˇ e indexy ´ struktur na ukázkové matici. Upozornuji, zaˇc´ınaj´ı od nuly. Mˇejme matici   1 6 9  2 7 0    13 3 8 0  A=    14 4 12  15 16 5 Potom obsah datovych ´ struktur v skyline systému bude : • diag = {1, 2, 3, 4, 5} • luc = 7 • uc = {6, 3, 7, 8, 9, 0, 0, 12} • iuc = {0, 0, 1, 2, 4, 8}

13

• llr = 4 • lr = {13, 14, 15, 16} • ilr = {0, 0, 0, 1, 2, 4}

Poznámka 3.1 Tento systém je statický. Jakmile je obsah datovych struktur sestaven v ´ pamˇeti, nˇekteré operace, jako je napˇr´ıklad zmˇena profilu matice pˇridán´ım nového prvku mimo profil, nebo rˇ a´ dková cˇ i sloupcová permutace je nemoˇzné provést bez realokace a znovusestaven´ı datovych ´ struktur. Alternativou by byl dynamický systém, ve kterém by prvky v profilu matice nebyly uloˇzeny v jednom poli, ale kaˇzdy´ rˇ a´ dek a sloupec by byl dynamicky alokované pole a m´ısto pole indexu˚ by se pouˇz´ıvalo pole ukazatelu˚ na tyto alokované pole. Délka jednotlivych ´ pol´ı by se musel uchovávat v samostatném poli, jelikoˇz tuto informaci by uˇz nebylo moˇzné z´ıskat z pole indexu˚ jako u statického systému. V tomto dynamickém systému by byly nˇekteré operace spojené se zmˇenou profilu mnohem ménˇe nákladné, protoˇze pro zmˇenu profilu jednotlivého rˇ a´ dku cˇ i sloupce staˇc´ı realokovat pouze pˇr´ısluˇsné pole a aktualizovat ukazatel, m´ısto realokace celé matice. Tento systém má i své nevyhody : vˇetˇs´ı pamˇet’ová nároˇcnost vlivem fragmentace ´ ˚ pamˇeti a niˇzsˇ´ı efektivita skalárn´ıch operac´ı s matic´ı (m´ısto jednoho pˇr´ımého pruchodu pole je tˇreba procházet vˇsechny pole). Vzhledem k tomu, zˇ e Skyline systém se pouˇz´ıvá vˇetˇsinou v aplikac´ıch, které zachovávaj´ı profil matice (Choleského faktorizace), implementoval jsem v této práci staticky´ systém.

3.2

ˇ ’ova´ naro ´ cnost ˇ Pamet

Pro uloˇzen´ı matice o rozmˇerech (m, n) ve skyline struktuˇre je potˇreba p + min{m, n} slov pamˇeti pro desetinná cˇ ´ısla (uloˇzen´ı prvku˚ matice v profilu a na diagonále) a m + n + 2 ˚ slov pamˇeti pro celá cˇ ´ısla (uloˇzen´ı indexu).

3.3

Pˇr´ıstup k prvkum ˚

Pˇr´ıstup k prvku na i-tém rˇ a´ dku a j-tém sloupci je pˇr´ımoˇcary´ - staˇc´ı zjistit, zda prvek nen´ı na diagonále, pˇr´ıpadnˇe mimo profil, a spoˇc´ıtat index, na kterém je uloˇzen v konkrétn´ım poli. Pˇr´ıstup je popsán v algoritmu 1.

3.4

´ Sestaven´ı skyline systemu

˚ Implementoval jsem dva zpusoby sestaven´ı skyline matice - z plné matice, a z vektoru˚ i, j, v (viz poznámka 2.3). Princip je u obou stejny´ : v prvn´ım kroku jsou urˇceny vektory sloupcovych ´ a rˇ a´ dkovych ´ profilu˚ (urˇcen´ı nenulového prvku s minimáln´ım sloupcovym, ´ pˇr´ıpadnˇe rˇ a´ dkovym ´ indexem). Poté je jiˇz moˇzné alokovat datové struktury a prvky matice um´ıstit na pˇr´ısluˇsná m´ısta v pamˇeti.

14

˚ matice ve skyline systému Algoritmus 1 Pˇr´ıstup k prvkum Require: 0 ≤ i < m and 0 ≤ j < n 1: if i = j then {diagon´ aln´ı prvek} 2: return diag[i] ´ 3: else if i < j then {horn´ı trojuheln´ ıková cˇ a´ st} 4: if j − i > cp (j) then {mimo profil} 5: return 0 6: else 7: return uc[iuc[j + 1] − j + i] 8: end if ´ 9: else {doln´ı trojuheln´ ıková cˇ a´ st} 10: if i − j > rp (i) then {mimo profil} 11: return 0 12: else 13: return lr[ilr[i + 1] − i + j] 14: end if 15: end if

3.5

´ Implementace skyline systemu pro knihovnu OOSol

Implementaci skyline systému jsem provádˇel do jiˇz existuj´ıc´ı knihovny OOSol. Skyline storage a základn´ı operace s matic´ı jsou implementovány v tˇr´ıdˇe OOSspSkylineStatic. Systém ukládán´ı je vyˇ ´ se popsany´ staticky´ skyline systém. Tˇr´ıda OOSspSkylineStatic je zaˇrazena do tˇr´ıdn´ı hierarchie knihovny OOSol, jako potomek abstraktn´ı tˇr´ıdy reprezentuj´ıc´ı rˇ´ıdké matice OOSspMatrix. OOSspMatrix je potomek tˇr´ıdy OOSspMatrix, abstraktn´ı obecné matice. Knihovna jiˇz obsahovala implementaci plné matice v tˇr´ıdˇe OOSfullMatrix a systému ukládán´ı rˇ´ıdké matice Compressed Row Storage v tˇr´ıdˇe OOSspMatrixStatic. Hierarchie tˇr´ıd je na obrázku 4. Poznámka 3.2 TIND je v knihovnˇe OOSol datovy´ typ pro celoˇc´ıselné indexy. Implementovány jsou následuj´ıc´ı metody : • OOSspSkylineStatic(TIND ndim) konstruktor, alokuje datové struktury pro nulovou cˇ tvercovou matici. • OOSspSkylineStatic(TIND nrow, TIND ncol) konstruktor, alokuje datové struktury pro nulovou obdéln´ıkovou matici. • OOSspSkylineStatic(OOSfullMatrix &mat) konstruktor, sestav´ı skyline matici z plné matice mat. • OOSspSkylineStatic(TIND nrow,TIND ncol, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V) konstruktor, sestav´ı skyline matici, která má nenulové prvky dané indexy ve vektorech I a J a jejich hodnoty ve V .

15

• OOSspSkylineStatic(TIND nrow, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V, OOSindexArray& p) konstruktor, sestav´ı skyline matici, která má nenulové prvky dané indexy ve vektorech I a J a jejich hodnoty ve V , permutovanou podle permutaˇcn´ıho vektoru p. • void clean() uvoln´ı vˇsechny asociované datové struktury a zavolá init(). • void init() nastav´ı matici na nulovou o rozmˇerech (0, 0). Datové struktory jsou nastaveny na nulovy´ ukazatel. • void resize(TIND ndim) uvoln´ı datové struktury a realokuje je pro cˇ tvercovou matici o rozmˇeru ndim. Prvky matice jsou nastaveny na nulové hodnoty. • void resize(TIND nrow, TIND ncol) uvoln´ı datové struktury a realokuje je pro obdélnikovou matici o rozmˇerech (nrow, ncol). Prvky matice jsou nastaveny na nulové hodnoty. • double get(TIND i,TIND j) vrac´ı hodnotu prvku matice na rˇ a´ dku i a sloupci j. • void set(TIND i,TIND j,const double& nval) nastav´ı hodnotu prvku matice na rˇ a´ dku i a sloupci j na nval. Pokud je prvek uvnitˇr profilu, je operace triviáln´ı. Pokud je ovˇsem prvek mimo profil, je realokována celá struktura pˇr´ısluˇsného horn´ıho nebo doln´ıho profilu, a jsou pˇrepoˇc´ıtány indexy. • void copy(OOSspSkylineStatic& mat, const double& alpha) nastav´ı matici jako kopii matice mat, skalárnˇe vynásobenou hodnotou alpha (A = αB). • void copyt(OOSspSkylineStatic& mat, const double& alpha) nastav´ı matici jako transponovanou kopii matice mat, skalárnˇe vynásobenou hodnotou alpha (A = αB T ). • void trans() transponuje matici. Tato operace je triviáln´ı, protoˇze staˇc´ı pouze ´ prohodit ukazatele na pole pro horn´ı a doln´ı trojuheln´ ıkovou cˇ a´ st. • void swap(OOSspSkylineStatic& mat) vymˇen´ı ukazatel na matici s ukazatelem na matici mat. • void operator=(OOSspSkylineStatic& mat) nastav´ı tuto matici jako mˇelkou kopii matice mat (tzn. jsou zkop´ırovány pouze ukazatele) a zavolá clean() na matici mat. • void operator*=(const double& alpha) provede na matici operaci skalárn´ıho násoben´ı reálnym ´ cˇ ´ıslem (A = αA). • void operator/=(const double& alpha) provede na matici operaci skalárn´ıho dˇelen´ı reálnym ´ cˇ ´ıslem (A = α1 A).

16

• void operator-=(const double& alpha) provede na matici operaci skalárn´ıho násoben´ı cˇ ´ıslem -1 (A = −A). • double norm 1() vrac´ı sloupcovou normu matice. Sloupcová norma je definována m P jako kAk1 = max |aij |). 1≤j≤n i=1

ˇ adková norma je definována • double norm inf() vrac´ı rˇ a´ dkovou normu matice. R´ n P jako kAk∞ = max |aij |. 1≤i≤m j=1

• double norm s F() vrac´ı Frobeniovu normu matice. Frobeniova norma je definována m P n P jako kAkF = |aij |2 . i=1 j=1

• double norm(const char i) vrac´ı jednu z vyˇ ´ se uvedenych ´ norem, podle hodnoty agrumentu. Pro i = 0 Frobeniovu, pro i = 1 sloupcovou, pro i = 3 rˇ a´ dkovou. ˚ typ • double diagNorm(const type) vrac´ı normu vektoru diagonáln´ıch prvk P u, normy je dán argumentem type. Pro hodnotu 1 se vrac´ı norma ||v||1 = |vi |, 1≤i≤dim v P pro hodnotu 2 norma ||v||2 = vi2 , pro hodnotu 3 norma ||v||∞ = max{|vi | : 1≤i≤dim v

1 ≤ i ≤ dim v} • ostream& operator<<(ostream& os, OOSspSkylineStatic& mat) vyp´ısˇ e matici mat do streamu os. Slouˇz´ı pˇredevˇs´ım k informaˇcn´ımu vypisu napˇr´ıklad do ´ cout • istream& operator>>(istream& is,OOSspSkylineStatic& mat) naˇcte matici ze souborového streamu is. Formát je textovy´ nebo binárn´ı (lze nastavit zavolán´ım metody setText() nebo setBinary()), stejny´ jako pouˇz´ıvá tˇr´ıda OOSspMatrixStatic. • ofstream& operator<<(ofstream& os,OOSspSkylineStatic& mat) naproti tomu zap´ısˇ e matici do souborového streamu os tak, aby ji bylo moˇzné naˇc´ıst ˚ ze byt operátorem >>. Formát muˇ ´ rovnˇezˇ textovy´ nebo binárn´ı. Tˇr´ıda obsahuje nav´ıc metody : • symmetricProfile vrac´ı true, pokud má matice symetricky´ profil. • TIND columnHeight(TIND i) vrac´ı profil i-tého sloupce. • TIND rowHeight(TIND i) vrac´ı profil i-tého rˇ a´ dku. • static int profile(TIND nrow, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V, OOSindexArray& p) statická metoda, slouˇz´ıc´ı k urˇcen´ı

17

Obrázek 4: Hierarchie tˇr´ıd reprezentuj´ıc´ıch matice v knihovnˇe OOSol profilu skyline matice dané indexy prvku˚ v I a J, jejich hodnotami ve V a permutaˇcn´ım vektorem p. Data pro matici nejsou alokována, pouze se poˇc´ıtaj´ı rˇ a´ dkové ˚ ze slouˇzit napˇr´ıklad k vhodnému urˇcen´ı parametru˚ W1 a W2 a sloupcové profily. Muˇ sloanova algoritmu.

18

4

Grafy

˚ algoritmus, stejnˇe jako napˇr´ıklad algoritmus Cuthill–McKee, je zaloˇzen na repreSloanuv zentaci matici neorientovanym ´ grafem. Nové indexy vrcholu˚ jsou urˇceny pˇri postupném ˚ nezbytné pro poprocházen´ı grafu. V následuj´ıc´ıch sekc´ıch uvedu pojmy z teorie grafu, ˚ pis tohoto algoritmu, poté pop´ısˇ i zpusob ukládán´ı grafu v poˇc´ıtaˇcové pamˇeti pomoc´ı ˚ Následuje popis samotného algoritmu a komentárˇ e k implementaci. seznamu sousedu.

4.1

´ Zakladn´ ı pojmy z teorie grafu˚

Definice 4.1 Jednoduchy´ neorientovany´ graf je uspoˇra´ daná dvojice G = (V, E) kde V je mnoˇzina vrcholu˚ a E je mnoˇzina hran. Hrany jsou neuspoˇra´ dané dvojice ruzn´ ˚ ych vrcholu. Poˇcet vrcholu˚ ||V || oznaˇcujeme jako stupenˇ grafu G. Hranu mezi vrcholy u a v ∈ V znaˇc´ıme (u, v) a rˇ´ıkáme, zˇ e vrcholy u a v jsou sousedn´ı. Nejsou povoleny hrany spojuj´ıc´ı vrchol sám se sebou nebo v´ıce neˇz jedna hrana mezi dvˇema danými vrcholy. ˚ zeme graficky znázornit tak, zˇ e vrcholy zobraz´ıme jako body, a vrcholy mezi Graf muˇ ´ ckou. kterymi existuje hrana, spoj´ıme useˇ ´ Definice 4.2 Stupenˇ vrcholu u oznaˇcujeme jako poˇcet hran obsahuj´ıc´ıch tento vrchol. Definice 4.3 Sled v grafu G je posloupnost vrcholu˚ a hran P = (v0 , e1 , v1 , . . . , en , vn ), pro kterou plat´ı ei = (vi−1 , vi ). (Existuje tedy hrana mezi vrcholem a následuj´ıc´ım vrcholem v posloupnosti) Definice 4.4 Cesta v grafu G je sled, pro který plat´ı vi 6= vj pro i 6= j (ˇza´ dný vrchol se neopakuje). Lemma 1 Mˇejme relaci ∼ na mnoˇzinˇe vrcholu˚ libovolného grafu G. Pro dva vrcholy u ∼ v, právˇe kdyˇz existuje v G cesta zaˇc´ınaj´ıc´ı v u a konˇc´ıc´ı ve v. Pak ∼ je relac´ı ekvivalence. ˚ Dukaz je triviáln´ı. Definice 4.5 Tˇr´ıdy ekvivalence relace ∼ popsané v lemma 1 se nazývaj´ı komponenty souvislosti grafu. ´ Definice 4.6 Vzdalenost mezi vrcholy u a v oznaˇcujeme jako poˇcet hran nejkratˇs´ı moˇzné cesty spojuj´ıc´ı tyto vrcholy. Definice 4.7 Pr˚ umˇer grafu je taková dvojice vrcholu˚ (u, v), která ma ze vˇsech moˇzných dvojic vrcholu˚ v grafu nejvˇetˇs´ı vzdálenost

19

Obrázek 5: Pˇr´ıklad grafu

4.2

Souvislost mezi grafem a strukturou matice

˚ zeme reprezentovat pomoc´ı symetrické matice A Jednoduchy´ neorientovany´ graf G muˇ ˚ následuj´ıc´ım zpusobem : • Matice A je cˇ tvercová matice o rozmˇerech (n, n), kde n je stupenˇ grafu G. • Pokud v grafu G existuje hrana mezi vrcholy u a v, pak auv = avu = 1, jinak auv = avu = 0 • Vˇsechny diagonáln´ı prvky jsou rovny jedné aii = 1, 1 ≤ i ≤ n Tuto matici nazyv´ ´ ame matic´ı sousednosti. Pˇr´ıklad 4.1 Na obrázku 5 je pˇr´ıklad jednoduchého neorientovaného grafu. Odpov´ıdaj´ıc´ı matice sousednosti k tomuto grafu je pak : 

1  1 A=  1 1

1 1 1 0

1 1 1 0

 1 0   0  1

Vˇsimnˇeme si, zˇ e stupenˇ vrcholu i je roven nnz(riA ) − 1 = nnz(sA i )−1 ˚ zeme pˇriˇradit graf G Obrácenˇe pak cˇ tvercové matici A se symetrickym ´ profilem muˇ ˚ t´ımto zpusobem : ˚ kde n je poˇcet rˇ a´ dku˚ matice A • Graf G má n vrcholu, • V grafu G je hrana mezi vrcholy i a j, pokud je prvek matice aij = aji , j 6= i nenulovy. ´

20

4.3

´ an´ ´ ı grafu pomoc´ı seznamu sousedu˚ Uklad

Pouˇz´ıvat pro ukládán´ı grafu rˇ´ıdké matice plnou matici sousednosti by bylo neefektivn´ı, ˚ proto se v implementaci sloanova algoritmu pouˇz´ıvá alternativn´ı zpusob reprezentace grafu : seznam sousedu, ˚ pro jehoˇz uloˇzen´ı je potˇreba pouze v+1+2e celoˇc´ıselnych ´ jednotek pamˇeti, kde v je poˇcet vrcholu˚ a e je poˇcet hran grafu (naproti tomu uloˇzen´ı plné matice sousednosti vyˇzaduje v 2 jednotek). Seznam sousedu˚ je do jisté m´ıry podobny´ systému ˚ pouze jejich indexy. ukládán´ı rˇ´ıdkych ´ matic, nepotˇrebujeme ale hodnoty prvku, Seznam sousedu˚ je implementován pomoc´ı dvou pol´ı : • Pole adj délky 2e slouˇz´ı k ukládán´ı indexu˚ sousedn´ıch vrcholu˚ pro vˇsechny vrcholy • Pole iadj je pole indexu˚ pro pole adj. Sousedn´ı vrcholy k vrcholu i najdeme v poli adj na indexech iadj[i] aˇz iadj[i+1]-1. Stupenˇ vrcholu i je tedy roven iadj[i+1]-iadj[i] Kaˇzdá hrana je tedy uchovávána dvakrát pro jednoduchost programu a rychlost pˇr´ıstupu.

4.4

Sestaven´ı seznamu sousedu˚

˚ Stejnˇe jako u sestaven´ı skyline systému jsem implementoval dva zpusoby sestaven´ı seznamu sousedu˚ - z plné matice a z vektoru˚ i, j, v. Sestaven´ı prob´ıhá ve dvou kroc´ıch - nej˚ ktery´ je dán poˇctem nenulovych prve je urˇcen stupenˇ kaˇzdého vrcholu, ´ prvku pˇr´ısluˇsného rˇ a´ dku nebo sloupce mimo diagonálu. V pˇr´ıpadˇe plné matice se procház´ı celá matice a testuje se, zda je prvek nulovy´ nebo nenulovy. ´ Hranu v grafu ukládáme, pokud je alesponˇ jeden z prvku˚ (aij , aji ) nenulovy´ - to znamená, zˇ e sestaven´ı grafu funguje i pro matice s nesymetrickym ´ profilem. Pouˇz´ıvá se ovˇsem profil, ktery´ je maximem profilu pˇr´ısluˇsného rˇ a´ dku a sloupce. Pˇri sestaven´ı grafu z vektoru˚ i, j, v procház´ıme pouze vektory i a j, tˇemito jsou jiˇz dány hrany v grafu (ik , jk ). ˚ je aloPoté, co je urˇcen celkovy´ poˇcet hran v grafu a stupnˇe jednotlivych vrcholu, ´ ˚ ˚ kováno pole adj a v dalˇs´ım kroku je dalˇs´ım pruchodem naplnˇeno indexy vrcholu.

4.5

ˇ a´ struktura Koˇrenova´ urov ´ nov

˚ algoritmus pouˇz´ıvá v prvn´ım kroku takzvanou koˇrenovou urov Sloanuv ´ novou ˇ strukturu. ´ ˇ a´ struktura pˇr´ısluˇs´ıc´ı vrcholu r je posloupnost V neprázdných Definice 4.8 Koˇrenova´ urov nov mnoˇzin vrcholu. ˚ V1 = {r}, Vi obsahuje vˇsechny vrcholy, které soused´ı s vrcholy ve Vi−1 a pˇritom nejsou v zˇ a´ dné z mnoˇzin V1 , . . . , Vi−1 . Vi se oznaˇcuje jako i-tá urove ´ n. ˇ ´ ˇ e struktury pro vrPoˇcet prvku˚ posloupnosti V oznaˇcujeme jako hloubku koˇrenové urov nov´ chol r. Poˇcet prvku˚ mnoˇziny ||Vi || oznaˇcujeme jako sˇ´ırˇku i-té urovnˇ ´ e. max ||Vi ||, maximum sˇ´ırˇek jednot´ ˇ e struktury pro vrchol r. livých urovn´ ´ ı oznaˇcujeme jako sˇ´ırˇ ku koˇrenové urov nov´

21

´ ˇ a struktura tedy na své i-té urovni ´ Koˇrenová urov nov´ obsahuje vrcholy, jejichˇz vzdálenost ´ ˇ e struktury odpov´ıdá od vrcholu r je i − 1. Algoritmus 2 pro sestaven´ı koˇrenové urov nov´ procházen´ı grafu do sˇ´ırˇ ky. Jiˇz navˇst´ıvené vrcholy jsou oznaˇceny v maskovac´ım vektoru, aby se zabránilo zacyklen´ı. ´ ˇ a struktura implementována jako dvˇe pole. Pole rls V C++ je pak koˇrenová urov nov´ obsahuje vˇsechny prvky struktury, pole irls jsou indexy v poli rls, kde zaˇc´ınaj´ı jednot´ livé urovnˇ e. Vytvoˇren´ı struktury zajiˇst’uje metoda OOSgraph::rootLevelStructure ´ ˇ e struktury Algoritmus 2 Sestaven´ı koˇrenové urov nov´ Require: r {poˇca´ teˇcn´ı vrchol} Require: n {stupenˇ grafu} ´ 1: w ← 1 {ˇs´ırˇ ka souˇcasné urovnˇ e} 2: wmax ← w {maxim´ aln´ı sˇ´ırˇ ka} 3: V1 ← {r} 4: d ← 1 {hloubka} 5: m = (0, . . . , 0), dim m = n {maskovac´ı vektor} 6: while w > 0 do 7: w←0 8: for all p in Vd do 9: s ← vrcholy soused´ıc´ı s p 10: for all si in s do 11: if msi = 0 then {dosud nenavˇst´ıveny´ vrchol} 12: w ←w+1 13: msi ← 1 ´ 14: insert si into Vd+1 {pˇr´ıdán´ı vrcholu do souˇcasné urovnˇ e} 15: end if 16: end for 17: end for 18: if w > 0 then 19: d←d+1 20: end if 21: wmax ← max{w, wmax } 22: end while 23: return V, d, wmax

22

5

Sloanuv ˚ algoritmus

˚ algoritmus, popsany´ v cˇ lánku [2], s menˇs´ımi upravami ´ Sloanuv v [3] slouˇz´ı ke sn´ızˇ en´ı ´ ˚ algoritmus je vyprofilu matice. Minimalizace profilu je NP-upln y´ problém [4]. Sloanuv ˚ lepˇsen´ım oproti pˇredchoz´ım algoritmum, jako je RCM (Reverse Cuthill–McKee). Vstupem sloanova algoritmu je graf G, vystupem je permutaˇcn´ı vektor p takovy, ´ ´ zˇ e symetrická permutace odpov´ıdaj´ıc´ı matice vektorem p má pokud moˇzno niˇzsˇ´ı profil. Al˚ eru grafu a pˇreˇc´ıslován´ı vrgoritmus se dá rozdˇelit do dvou cˇ a´ st´ı - urˇcen´ı pseudoprumˇ ˚ cholu.

5.1

ˇ ı pseudoprum ˇ Urcen´ ˚ eru

˚ algoritmus vyˇzaduje dva vrcholy, které urˇcuj´ı poˇca´ tek a konec cˇ ´ıslován´ı. Ideálnˇe Sloanuv ˚ er grafu, urˇcen´ı prumˇ ˚ eru grafu je vˇsak pˇr´ıliˇs vypoˇ jsou tyto vrcholy prumˇ ´ cetnˇe nároˇcné. ˚ algoritmus proto pouˇz´ıvá takzvany´ pseudoprumˇ ˚ jejichˇz Sloanuv ˚ er, coˇz je dvojice vrcholu, ˚ eru grafu. vzdálenost je velmi bl´ızká, nebo rovna vzdálenost´ı vrcholu˚ prumˇ ˚ eru slouˇz´ı heuristicky´ algoritmus 3, ktery´ se snaˇz´ı naj´ıt dvojici K nalezen´ı pseudoprumˇ vrcholu˚ s maximáln´ı vzdálenost´ı a minimáln´ı sˇ´ırˇ kou grafu. Algoritmus vˇetˇsinou najde ˚ er v nˇekolika málo iterac´ıch. V C++ je tento algoritmus implementován v pseudoprumˇ ´ metodˇe OOSgraph::pseudoDiameter. K setˇr´ıdˇen´ı posledn´ı urovnˇ e je pouˇzit insertion ˚ sort, jelikoˇz tato má obvykle nˇekolik málo vrcholu.

5.2 5.2.1

ˇ ıslovan´ ´ ı vrcholu˚ Pˇrec´ Stav vrcholu

˚ algoritmus jeden z tˇechto stavu˚ : postactive, Bˇehem pˇreˇc´ıslován´ı vrcholu˚ pˇriˇrazuje sloanuv active, preactive, inactive. Jejich vyznam jsou: ´ • Vrchol má stav postactive, pokud mu jiˇz bylo pˇriˇrazeno nové cˇ ´ıslo. • Vrcholy, které soused´ı s postactive vrcholy, ale jeˇstˇe jim nebylo pˇriˇrazeno cˇ ´ıslo, jsou active. • Vrcholy, které soused´ı s vrcholy stavu preactive, a pˇritom nejsou postactive ani active, maj´ı stav preactive • Ostatn´ı vrcholy maj´ı stav inactive. Stavy vrcholu˚ tvoˇr´ı vektor σ. Ukázka stavu vrcholu˚ je na obrázku 6. 5.2.2

´ ı stupenˇ vrcholu Momentaln´

Kaˇzdému vrcholu i je také pˇriˇrazen momentáln´ı stupen, ˇ mi ktery´ je definován takto : • Vrchol stavu postactive má mi = 0

23

˚ eru grafu Algoritmus 3 Nalezen´ı pseudoprumˇ 1: s ← vrchol s nejniˇ zsˇ´ım stupnˇem {odhad poˇca´ teˇcn´ıho vrcholu} 2: snew ← 1 3: while snew = 1 do {dokud je volen novy ´ poˇca´ teˇcn´ı bod} 4: snew ← 0 ´ ˇ a struktura z vrcholu s o hloubce d a sˇ´ırˇ ce w 5: V ← koˇrenová urov nov´ 6: setˇrid’ Vd podle stupnˇe vrcholu ´ ´ 7: Q ← Vd tak, zˇ e r neobsahuje dva vrcholy stejné urovnˇ e {zkrácen´ı posledn´ı urovnˇ e} 8: wmin = ∞ 9: for all qi in Q do ´ ˇ a struktura z vrcholu qi o hloubce dq a sˇ´ırˇ ce wq 10: U ← koˇrenová urov nov´ 11: if wq < wmin then {nalezena struktura s menˇs´ı sˇ´ırˇ kou } 12: if dq > d then 13: s ← qi {novy´ poˇca´ teˇcn´ı vrchol} 14: snew ← 1 15: d ← dq 16: break 17: else 18: e ← qi {novy´ koncovy´ vrchol} 19: wmin ← wq 20: end if 21: end if 22: end for 23: end while 24: return s, e ˚ • Momentáln´ı stupenˇ vrcholu se stavem active je roven poˇctu jeho preactive sousedu. • Stupenˇ preactive vrcholu je roven souˇctu inactive a preactive sousedu˚ plus jedna. • Pro vrcholy stavu inactive, je momentáln´ı stupenˇ roven stupni vrcholu plus jedna. Na zaˇca´ tku algoritmu maj´ı vˇsechny vrcholy momentáln´ı stupenˇ roven stupni vrcholu plus jedna, a s postupem algoritmu, kdyˇz jsou oˇc´ıslovávány sousedn´ı vrcholy, momentáln´ı stupenˇ klesá. Na konci algoritmu maj´ı vˇsechny vrcholy nulovy´ momentáln´ı ˇ stupen. Napˇr´ıklad na obrázku 6 jsou momentáln´ı stupnˇe vrcholu˚ m1 = 0, m2 = 0, m3 = 0, m4 = 2, m5 = 2, m6 = 3, m7 = 4, m8 = 3, m9 = 5, m10 = 6, m11 = 4, m12 = 4. 5.2.3

Priorita vrcholu

Pˇri volbˇe dalˇs´ıho vrcholu urˇceného k pˇreˇc´ıslován´ı je rozhoduj´ıc´ı priorita vrcholu. Kaˇzdému vrcholu je pˇriˇrazena priorita, která je dána dvˇema sloˇzkami - momentáln´ım stupnˇem a

24

Obrázek 6: Stavy vrcholu˚ vzdálenost´ı od koncového vrcholu. Priorita vrcholu i je definována jako pi = W1 · δi − W2 · mi ´ ˇ e kde δi je vzdálenost vrcholu od koncového vrcholu (tu zjist´ıme z koˇrenové urov nov´ struktury sestavené z koncového vrcholu), mi je momentáln´ı stupenˇ vrcholu, a W1 , W2 jsou váhové parametry. 5.2.4

Algoritmus

˚ algoritmus pro pˇreˇc´ıslován´ı vrcholu. ˚ Vstupem Algoritmus 4 popisuje samotny´ sloanuv ˚ eru grafu. Vystupem algoritmu je graf G a poˇca´ teˇcn´ı a koncové vrcholy psedoprumˇ algo´ ˚ které by mˇelo sn´ızˇ it profil odpov´ıdaj´ıc´ı matice. Bˇehem ritmu je nové oˇc´ıslován´ı vrcholu, bˇehu algoritmu se udrˇzuje fronta vrcholu˚ q se stavem preactive nebo active urˇcenych ´ k pˇreˇc´ıslován´ı, z nich se vˇzdy vybere ten s maximáln´ı prioritou a ten je oˇc´ıslován. Poté se aktualizuje stav a priorita sousedn´ıch vrcholu˚ a vrcholy, jejichˇz stav se zmˇenil z inactive na preactive jsou pˇridány do fronty. Algoritmus konˇc´ı, kdyˇz jsou pˇreˇc´ıslovány vˇsechny ˚ vrcholy, vystupem je inverzn´ı permutaˇcn´ı vektor η urˇcuj´ıc´ı nové oˇc´ıslován´ı vrcholu. ´

5.3

Hodnoty parametru˚ W1 a W2

Volba hodnot parametru˚ W1 a W2 je zásadn´ı pro kvalitu vysledn´ eho pˇreˇc´ıslovan´ı. Pa´ rametr W1 v prioritn´ı funkci urˇcuje globáln´ı prioritu, jenˇz je nemˇenná bˇehem celého alˇ goritmu a upˇrednostnuje vrcholy vzdálenˇejˇs´ı od koncového vrcholu. Zat´ımco parametr ˇ W2 reprezentuje lokáln´ı prioritu, upˇrednostnuje vrcholy s n´ızkym ´ stupnˇem soused´ıc´ı s jiˇz oˇc´ıslovanymi vrcholy. Pokud zvol´ ı me W = 6 0, W = 0 algoritmus funguje velmi podobnˇe ´ 1 2 jako RCM. Sloan [3] doporuˇcuje na základˇe empirickych ´ testu˚ hodnoty W1 = 1, W2 = 2, která jsou vhodné pro mnoho matic.

25

˚ e typy matic. V kapitole 6 je Optimáln´ı hodnoty ovˇsem mohou byt ´ jiné pro ruzn´ napˇr´ıklad matice, která má s parametry W1 = 1, W2 = 1 niˇzsˇ´ı vysledn y´ profil neˇz pˇri ´ pouˇzit´ı parametru˚ W1 = 1, W2 = 2. ˚ e Pˇri praktickém vyuˇzit´ı Sloanova algoritmu je tedy vhodné nejprve vyzkouˇset ruzn´ W2 pomˇery hodnot W1 a vybrat hodnoty parametru˚ optimáln´ı pro dany´ typ problému, nebo automaticky vybrat z nˇekolika pˇredem urˇcenych ´ hodnot ty, pro které je profil nejmenˇs´ı.

5.4

Implementace sloanova algoritmu pro knihovnu OOSol

˚ algoritmus pro pˇreˇc´ıslován´ı vrcholu˚ jsem implementoval v Sestaven´ı grafu a sloanuv tˇr´ıdˇe OOSgraph. Tato tˇr´ıda reprezentuje jednoduchy´ neorientovany´ graf pomoc´ı seznamu ˚ Tˇr´ıda dále obsahuje metody pro sestaven´ı koˇrenové urov ´ ˇ e struktury, nalesousedu. nov´ ˚ eru a koneˇcnˇe nalezen´ı permutaˇcn´ıho vektoru sloanovym zen´ı pseudoprumˇ ´ algoritmem. Tato implementace dokázˇ e zpracovat grafy s v´ıce neˇz jednou komponentou souvislosti, vˇcetnˇe osamocenych ´ vrcholu˚ (vrcholy stupnˇe nula). Osamocené vrcholy jsou vˇzdy oˇc´ıslovány aˇz jako posledn´ı. Oznaˇcován´ı jednotlivych ´ komponent je rˇ eˇseno pomoc´ı maskovac´ıho vektoru, kterym ´ jsou oznaˇceny jiˇz oˇc´ıslované vrcholy grafu. Fronta vrcholu˚ s preactive a active stavem q je implementována jako pole. Pˇridán´ı a odebrán´ı vrcholu je operace se sloˇzitost´ı O(1), vyhledán´ı prvku s nejvyˇssˇ´ı prioritou má sloˇzitost O(n). Tato implementace je dostaˇcuj´ıc´ı pro grafy s nˇekolika tis´ıci vrcholy, pro rozsáhlejˇs´ı problémy je vhodnˇejˇs´ı pouˇz´ıt sofistikovanˇejˇs´ı datovou strukturu, napˇr´ıklad binárn´ı strom. Seznam privátn´ıch metod tˇr´ıdy OOSgraph : • TIND nodeDegree(TIND i) vrac´ı stupenˇ vrcholu i. • void sortNodeByDegree(TIND* array, TIND start, TIND end) slouˇz´ı ˚ eru. Je k setˇr´ıdˇen´ı vrcholu˚ podle jejich stupnˇe pˇri algoritmu hledán´ı pseudoprumˇ ´ pouˇzit insertion sort, jelikoˇz ve zkrácené posledn´ı urovn´ ı byv´ ´ a obvykle jen nˇekolik ˚ i. málo vrcholu. • void rootLevelStructure(TIND root, TIND* mask, TIND* rls, TIND* irls, TIND &rlslen, TIND &depth, TIND &width ) slouˇz´ı k se´ ˇ e struktury z vrcholu root. staven´ı koˇrenové urov nov´ • number(int W1, int W2, TIND s, TIND* rls, TIND rlslen, int *status, TIND *q, int *priority, int &labeled) implementuje ˚ algoritmus pro jednu komponentu souvislosti grafu. Pole jsou samotny´ sloanuv pˇredávány jako ukazatele, aby se zamezilo jejich nˇekolikanásobné alokaci bˇehem bˇehu sloanova algoritmu pro graf s nˇekolika komponentami souvislosti. Veˇrejné metody jsou : • OOSgraph(OOSmatrix &mat) konstruktor - sestav´ı graf z plné cˇ tvercové matice mat

26

• OOSgraph(TIND nrow, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V) konstruktor - sestav´ı graf z matice dané indexovymi vektory I, J a vektorem ´ hodnot V ˚ algoritmus pro tento graf s hodno• OOSindexArray* sloan() provede sloanuv tami parametru˚ W1 = 1, W2 = 2. Vrac´ı ukazatel na permutaˇcn´ı vektor. ˚ algoritmus se spe• OOSindexArray* sloan(int W1, int W2) provede sloanuv ˚ cifikovanymi hodnotami parametru. ´

27

˚ algoritmus pro pˇreˇc´ıslován´ı vrcholu˚ Algoritmus 4 Sloanuv ˚ er grafu} Require: (s, e) {pseudoprumˇ ´ ˇ a struktura z vrcholu e 1: V ← koˇrenov´ a urov nov´ 2: for all Vi in V do 3: for all sj in Vi do ´ 4: δsi ← i {urˇcen´ı vzdálenosti vrcholu od koncového vrcholu e podle jeho urovnˇ e ´ ˇ e struktuˇre} v koˇrenové urov nov´ 5: end for 6: end for ´ teˇcn´ı priority a stavu} 7: for i = 1 to n do {pˇriˇrazen´ı poˇca 8: pi ← W1 · δi − W2 · mi 9: σi ← inactive 10: end for ˚ eru} 11: insert s into q {pˇreˇc´ıslov´ an´ı zaˇc´ıná od poˇca´ teˇcn´ıho vrcholu psedoprumˇ 12: σs ← preactive 13: l ← 1 {nové cˇ ´ıslo vrcholu} 14: while q not empty do {dokud jsou ve frontˇe vrcholy k pˇreˇc´ıslov´ an´ı} 15: m, takovy´ vrchol v q, pro ktery´ je pqi maximáln´ı 16: delete m from q 17: if σm = preactive then 18: for all j soused´ıc´ı s m do 19: pj ← pj + W2 {coˇz odpov´ıdá sn´ızˇ en´ı momentáln´ıho stupnˇe mj o 1} 20: if σj = inactive then 21: insert j into q 22: σj ← preactive 23: end if 24: end for 25: end if 26: ηi ← l 27: l ←l+1 28: for all i soused´ıc´ı s m do 29: if σj = preactive then 30: pj ← pj + W2 31: σj ← active 32: for all k soused´ıc´ı s j do 33: p j ← p j + W2 34: if σk = inactive then 35: σj ← preactive 36: insert k into q 37: end if 38: end for 39: end if 40: end for 41: end while 42: return η

28

6

´ ı a vysledky Testovan´ ´

Implementaci algoritmu jsem otestoval na nˇekolika rˇ´ıdkych uvád´ım ´ matic´ıch, vysledky ´ ˚ v tabulce 2. Tabulka obsahuje puvodn´ ı profil matice, profil po sloanovˇe pˇreuspoˇra´ dán´ı s parametry W1 , W2 a pro srovnán´ı profil po symetrickém RCM pˇreuspoˇra´ dán´ı, tak jak je implementováno v Octave. Nˇekteré matice jsou vygenerované s´ıtˇe, jiné jsou reálné matice pˇrevzaté z kolekce rˇ´ıdkych ´ matic University of Florida Sparse Matrix Collection dostupné na adrese http://www.cise.ufl.edu/research/sparse/matrices/. U dvou matic je uvedeno grafické znázornˇen´ı pomoc´ı pˇr´ıkazu spy. Vysledky ukazuj´ı, zˇ e pro matice s velmi jednoduchou strukturou, jako je obdelnik ´ ˚ algoritmus a RCM stejny´ vysledn nebo krychle1 dává sloanuv y´ profil. U matic se ´ sloˇzitˇejˇs´ım profilem se jiˇz projevuje vyhoda pˇriˇrazen´ı priority vrcholu˚ (algoritmus RCM ´ ˚ algoritmus znatelnˇe niˇzsˇ´ı propouze procház´ı graf do hloubky), a pro tyto dává sloanuv fil. Za povˇsimnut´ı stoj´ı matice UF52, která pro hodnoty parametru˚ W1 = 1, W2 = 2 má ˚ vyˇssˇ´ı profil neˇz puvodn´ ı matice, zato s parametry W1 = 1, W2 = 1 niˇzsˇ´ı profil neˇz RCM. To ukazuje, zˇ e hodnoty parametru˚ W1 = 1, W2 = 2 nejsou univerzáln´ı pro vˇsechny typy ˚ e hodnoty parametru. ˚ matic a je vhodné pro konkrétn´ı typy matic otestovat ruzn´ Seznam pouˇzitych ´ matic je uveden v tabulce 1.

29

Obrázek 7: S´ıt’ matice ctverec1

Obrázek 8: S´ıt’ matice ctverec3

Obrázek 9: S´ıt’ matice krychle1



30


Obrázek 13: S´ıt’ matice UF871 - helikoptéra Commanche

˚ Obrázek 14: Puvodn´ ı profil matice krychle4

31

sˇ sˇ´ı matice obdelnik

rozmˇer 4000

nnz 27122

ctverec1 ctverec3 krychle1 krychle2 krychle3 krychle4 UF52

60 12054 21 195 1581 12294 28924

370 78518 39 1471 13263 142442 2043492

UF111

1007

8575

UF267

479

1888

UF871

7920

31680

popis obdéln´ık o 200 × 20 vrcholech spo˚ jenych ıkovou s´ıt´ı. ´ trojuheln´ diskretizace 2D cˇ tverce. Viz obr. 7 diskretizace 2D cˇ tverce. Viz obr. 8 diskretizace 3D krychle. Viz obr. 9 diskretizace 3D krychle. Viz obr. 10 diskretizace 3D krychle. Viz obr. 11 diskretizace 3D krychle. Viz obr. 11 STIFFNESS MATRIX FOR OFFSHORE GENERATOR PLATFORM (MSC NASTRAN) - PATTER SYMMETRIC CONNECTION TABLE FROM DTNSRDC, WASHINGTON U 8 STAGE COLUMN SECTION, ALL SECTIONS RIGOROUS ( CHEM. ENG. ), nesymetrická matice, velmi obt´ızˇ ná pro iteraˇcn´ı numerické metody S´ıt’ helikoptéry Commanche, viz obr. 13

Tabulka 1: Testovac´ı matice matice obdelnik ctverec1 ctverec1 ctverec3 ctverec3 krychle1 krychle2 krychle3 krychle4 UF52 UF52 UF111 UF267 UF871

profil 1527960 2150 2150 114487702 114487702 18 24786 1685816 110541680 31373936 31373936 51572 77693 33086588

W1 , W2 1,2 1,2 1,1 1,2 1,1 1,2 1,2 1,2 1,2 1,2 1,1 1,2 1,2 1,2

profil sloan 155020 806 822 2164978 2250146 18 6146 169992 5455392 34510616 28365376 43480 57268 860070

profil RCM 155020 914 914 2469134 2469134 18 6878 245034 8529714 44006432 44006432 47522 81110 1059214

Tabulka 2: Vysledn´ e profily ´

poznámka

obr. 17, 18

obr. 14, 15, 16

32

Obrázek 15: Profil matice krychle4 po pˇreˇc´ıslován´ı sloanovym ´ algoritmem

Obrázek 16: Profil matice krychle4 po pˇreˇc´ıslován´ı algoritmem RCM

˚ Obrázek 17: Puvodn´ ı profil matice ctverec3

Obrázek 18: Profil matice ctverec3 po pˇreˇc´ıslován´ı sloanovym ´ algoritmem

33

7

´ er ˇ Zav

˚ algoritmus se ukázal jako velmi vhodny´ pro pˇreuspoˇra´ dán´ı matic pˇred jejich Sloanuv ˚ zeme takto efektivnˇe uloˇzit i velmi rozsáhlé rˇ´ıdké mauloˇzen´ım ve skyline systému. Muˇ tice, které by bylo nemoˇzné uloˇzit v poˇc´ıtaˇcové pamˇeti jako plné matice. Velikost profilu je závislá na volbˇe hodnot parametru˚ W1 , W2 , jenˇz nelze urˇcit uni˚ e matice. verzálnˇe a jejichˇz optimáln´ı hodnoty se liˇs´ı pro ruzn´ Knihovnu OOSol by bylo moˇzné dále rozˇs´ırˇ it o efektivn´ı implementace nˇekterych ´ ˚ které zachovávaj´ı skyline profil, napˇr´ıklad choleského faktorizace nebo LU algoritmu, rozklad a umoˇznit tak praktické pouˇzit´ı skyline systému pro rˇ eˇsen´ı lineárn´ıch soustav. ˚ algoritmus by bylo moˇzno dále zefektivnit vylepˇsen´ım datové struktury Sloanuv fronty, nab´ızej´ı se i experimenty s prioritn´ı funkc´ı. Algoritmus by bylo moˇzné rozˇs´ırˇ it z neorientovaného grafu na orientovany´ graf reprezentuj´ıc´ı matice s nesymetrickym ´ profilem. V cˇ lánku [4] je popsáno rozˇs´ırˇ en´ı algoritmu na ohodnocené grafy. Martin Stachonˇ

34

8

Reference

ˇ [1] Zdenˇek Dostál, Lineárn´ı algebra, VSB-TUO 2000 [2] S. W. Sloan, An algorithm for profile and wavefront reduction of sparse matrices, International Journal for Numerical Methods in Engineering Vol. 23, John Wiley & Sons, Ltd., 1986 [3] S. W. Sloan, A Fortran program for profile and wavefront reduction, International Journal for Numerical Methods in Engineering Vol. 28, John Wiley & Sons, Ltd., 1989 [4] Gary Kumfert, Alex Pothen, Two improved algorithms for envelope and wavefront reduction, NASA, 1997

35

A

´ Zdrojove´ kody

´ ˚ algoritmus v V pˇr´ıloze uvád´ım zdrojové kody implementuj´ıc´ı skyline storage a sloanuv rámci knihovny OOSol.

36

5

10

15

/********************************************************************/ /* */ /* OOSol2 2005 */ /* Dept. of Applied Mathematics, */ /* Technical University in Ostrava, Czech Rep. */ /* */ /* Filename: OOSmatricesSkylineStatic.h */ /* Created by Martin Stachon, 2008 */ /* */ /* List of updates: */ /* Date Author Description */ /* */ /********************************************************************/ #ifndef OOSMATRICESSKYLINESTATIC_H #define OOSMATRICESSKYLINESTATIC_H #include "OOSmatrices.h" namespace oosol {

20 class OOSspSkylineStatic: public OOSspMatrix { friend class OOSgraph; 25

/** diagonal elements */ double *diag; /** upper triangle columns */ double *uc;

30 /** index */ TIND *iuc;

35

/** number of upper triangle items */ TIND luc; /** lower triangle rows */ double *lr;

40

/** index */ TIND *ilr; /** number of lower triangle items */ TIND llr;

45 public: /** Default constructor - initiates empty matrix */ OOSspSkylineStatic(): OOSspMatrix() {}; 50 /** Constructor - initiates square matrix */ OOSspSkylineStatic(TIND ndim);

37

55

/** Constructor - initates rectangular matrix */ OOSspSkylineStatic(TIND nrow,TIND ncol); /** Constructor - creates sparse matrix from full matrix mat */ OOSspSkylineStatic(OOSfullMatrix &mat);

60

65

70

75

80

85

90

95

100

/** Copy constructor - copies all member data from the matrix mat */ OOSspSkylineStatic(OOSspSkylineStatic& mat): OOSspMatrix(mat) {init(); copy(mat);} /** Constructor - assembles matrix with nonzeros given by vectors of indices I, J and values V */ OOSspSkylineStatic(TIND nrow,TIND ncol, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V); /** Constructor - assembles rectangular matrix with nonzeros given by vectors of indices I, J and values V, rows and columns are permuted with permutation verctor p */ OOSspSkylineStatic(TIND nrow, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V, OOSindexArray& p); /** Destructor - deallocates memory reserved for member data */ virtual ÕOSspSkylineStatic(); /** CLEAN deallocates memory reserved for member data and initiates them by 0 or NULL */ virtual void clean(); /** INIT initiates all member data by 0 or NULL */ virtual void init(); /** RESIZE(NDIM) removes allocated matrix structures and resizes matrix to the square one with ndim rows and columns */ virtual void resize(TIND ndim); /** RESIZE(NROW,NCOL) removes allocated matrix structures and resizes matrix to the one with nrow rows and ncol columns */ virtual void resize(TIND nrow,TIND ncol); /** SPCONVERT converts the nonzero entries described by vectors I,J,V to sparse matrix */ void spconvert(OOSindexArray& I, OOSindexArray& J, OOSvectorD& V, TIND nrow, TIND ncol) {OOSspSkylineStatic mat(nrow,ncol, I,J,V); *this=mat;} /** SLOAN converts the nonzero entries described by vectors I,J,V to sparse matrix, using Sloan’s algorithm to reoder matrix */ OOSindexArray* sloan(OOSindexArray& I, OOSindexArray& J, OOSvectorD& V, TIND nrow); /** GET virtual /** SET virtual

returns value from the i-th row and the j-th column */ double get(TIND i, TIND j); sets value in the i-th row and the j-th column to value nval */ void set(TIND i, TIND j, const double& nval);

/** COPY matrix copy and scale MatS <- alpha*mat */ void copy(OOSspSkylineStatic& mat, const double& alpha=1.0);

38

105

110

115

120

125

130

135

140

145

150

/** COPY transposed matrix copy and scale MatS <- alpha*matˆT */ void copyt(OOSspSkylineStatic& mat, const double& alpha=1.0); /** SWAP swaps matrices MatS <-> mat */ void swap(OOSspSkylineStatic& mat); /* MATRIX OPERATIONS */ /* ================= */ /** OPERATOR= (MatS=mat) assigns member data from the matrix mat to MatS and initiates mat to empty vector */ void operator=(OOSspSkylineStatic& mat); /** OPERATOR*=(mat) cummulative matrix-matrix multiplication MatS<-MatS *mat */ void operator*=(OOSspSkylineStatic& mat) {cmm(mat);} /** OPERATOR*=(alpha) matrix-scalar multiplication MatS<-alpha*MatS */ virtual void operator*=(const double& alpha); /** OPERATOR/=(alpha) matrix-scalar dividing MatS<-MatS/alpha */ virtual void operator/=(const double& alpha); /** OPERATOR˜ matrix transposition MatS <- MatSˆT */ virtual void operator˜() {copyt(*this);}; /** TRANS matrix transposition MatS <- MatSˆT */ virtual void trans(); /** OPERATOR- (unary) changes the signs of the vector entries MatS <(-MatS) */ virtual void operator-();

/* MATRIX NORMS */ /* ============ */ /** NORM_1 returns the row norm */ virtual double norm_1(); /** NORM_2 returns the spectral norm */ /** NORM_INF returns the column norm */ virtual double norm_inf(); /** NORM_F returns the Frobenius norm */ virtual double norm_F(); /** NORM matrix norm ||A||_i i==0 ... ||A||_F (Frobenius norm) i==1 ... ||A||_1 (row norm) i==2 ... ||A||_2 (spectral norm) i==3 ... ||A||_inf (column norm) */ virtual double norm(const char i=0); /** NORM returns a norm of the diagonal, type defines type of the norm type==1 -> l_1 type==2 -> l_2 type==3 -> l_inf <default> */ virtual double diagNorm(const char type=3); /** returns true if the matrix is square and has a symmetric profile */ virtual bool symmetricProfile(); /** skylineProfile returns the size of skyline profile */ TIND skylineProfile() { return luc+llr;

39

155

}

160

/** columnHeight return the column profile for column i */ TIND columnHeight(TIND i) { ASSERT(i>=0 && i
165

/** columnHeight return the column row for row i */ TIND rowHeight(TIND i) { ASSERT(i>=0 && i
170

/** return the profile of Skyline matrix given by vectors I,J,V and permutation vector p*/ static int profile(TIND nrow, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V, OOSindexArray& p); 175

};

/* STREAM OPERATORS OVERLOADING */ /* ============================ */ 180

185

190

/** Overloaded input stream operator */ istream& operator>>(istream& is, OOSspSkylineStatic& mat); /** Overloaded output stream operator */ ostream& operator<<(ostream& os, OOSspSkylineStatic& mat); /** Overloaded input file stream operator */ ifstream& operator>>(ifstream& is, OOSspSkylineStatic& mat); /** Overloaded output file stream operator */ ofstream& operator<<(ofstream& os, OOSspSkylineStatic& mat); } #endif /*OOSMATRICESSKYLINESTATIC_H*/

Vypis 1: Soubor OOSmatricesSkylineStatic.h ´

3

8

/********************************************************************/ /* */ /* OOSol2 2005 */ /* Dept. of Applied Mathematics, */ /* Technical University in Ostrava, Czech Rep. */ /* */ /* Filename: OOSmatricesSkylineStatic.cpp */ /* Created by Martin Stachon, 2008 */ /* */ /* List of updates: */ /* Date Author Description */

40

13

/* */ /********************************************************************/ #include "OOSmatricesSkylineStatic.h" namespace oosol {

18

23

/*******************************************************************/ /* */ /* Method definitions of class OOSmatricesSkylineStatic */ /* */ /*******************************************************************/ /* Constructor - initiates square matrix */ OOSspSkylineStatic::OOSspSkylineStatic(TIND ndim): OOSspMatrix(ndim) {

28 if (row==0 && col==0) init(); else { diag = new double[row]; ASSERT(diag!=NULL); iuc = new TIND[row+1]; ASSERT(iuc!=NULL); ilr = new TIND[row+1]; ASSERT(ilr!=NULL); llr = 0; luc = 0;

33

for (TIND i=0; i
38

43 } }

48

53

/** Constructor - initates rectangular matrix */ OOSspSkylineStatic::OOSspSkylineStatic(TIND nrow, TIND ncol): OOSspMatrix( nrow, ncol) { if (row==0 && col==0) init(); else { TIND i; row = nrow; col = ncol; diag = new double[MIN(row,col)]; ASSERT(diag!=NULL); iuc = new TIND[col+1]; ASSERT(iuc!=NULL); ilr = new TIND[row+1]; ASSERT(ilr!=NULL);

58

63

for (i=0; i<MIN(row,col); i++) { diag[i] = 0; } for (i=0; i<=col; i++) iuc[i] = 0;

41

for (i=0; i<=row; i++) ilr[row] = 0; } 68

73

78

83

} /* Constructor - creates sparse matrix from full matrix mat */ OOSspSkylineStatic::OOSspSkylineStatic(OOSfullMatrix &mat) { TIND i,j; TIND ucpos=0, lrpos=0; bool inSkyline; double v; row = mat.getRows(); col = mat.getCols(); diag = new double[MIN(row,col)]; ASSERT(diag!=NULL); iuc = new TIND[col+1]; ASSERT(iuc!=NULL); ilr = new TIND[row+1]; ASSERT(ilr!=NULL); // copy diagonal for (i=0; i<MIN(row,col); i++) { diag[i] = mat(i,i); }

88

93

98

103

108

113

118

luc = llr = 0; // find number of elements under skyline profile in upper triangle for (j=1; j
42

inSkyline = true; if (inSkyline) uc[ucpos++] = v; } } iuc[col] = ucpos; // store non-zero elements in lower triangle for (i=1; i
123

128

133

138

143

148

153

158

} /* Constructor - assembles matrix with nonzeros given by vectors of indices I, J and values V */ OOSspSkylineStatic::OOSspSkylineStatic(TIND nrow,TIND ncol, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V): OOSspMatrix(nrow,ncol) { TIND i, N=I.getDim(); if (row==0 && col==0) init(); else { row = nrow; col = ncol; diag = new double[MIN(row,col)]; ASSERT(diag!=NULL); iuc = new TIND[col+1]; ASSERT(iuc!=NULL); ilr = new TIND[row+1]; ASSERT(ilr!=NULL); luc = 0; llr = 0; // initialize diagonal for (i=0; i<MIN(row,col); i++) { diag[i] = 0.0; }

163 TIND* ch = new TIND[col]; // column heights TIND* rh = new TIND[row]; // row heights TIND h; 168

for (i=0; i
43

173

178

183

188

193

198

203

208

213

218

223

for (i=0; i=0 && Ii= 0 && Ji < ncol); // in first pass, we place items on diagonal, and determine // the heights of skyline profile if (Ii == Ji) { // diagonal diag[Ii] = V(i); } else if (V(i) != 0.0) { if (Ii<Ji) { // upper triangle h = Ji-Ii; if (h>ch[Ji]) ch[Ji] = h; } else { // lower triangle h = Ii-Ji; if (h>rh[Ii]) rh[Ii] = h; } } } // now we can construct skyline storage iuc[0] = 0; ilr[0] = 0; for (i=1; i<=col; i++) { iuc[i] = iuc[i-1]+ch[i-1]; } for (i=1; i<=row; i++) { ilr[i] = ilr[i-1]+rh[i-1]; } luc = iuc[col]; llr = ilr[row]; uc = new double[luc]; lr = new double[llr]; for (i=0; i
44

delete[] ch; delete[] rh; } 228

}

233

/** Constructor - assembles rectangular matrix with nonzeros given by vectors of indices I, J and values V, rows and columns are permuted with permutation verctor p */ OOSspSkylineStatic::OOSspSkylineStatic(TIND nrow, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V, OOSindexArray& p): OOSspMatrix(nrow) { TIND i, N=I.getDim(); TIND ii, ji;

238 ASSERT(p.getDim() == nrow); ASSERT(I.getDim() == J.getDim()); // inverse permutation vector OOSindexArray ip(nrow); for (i=0; i
243

248

OOSindexArray Ip(N), Jp(N);

253

// renumber indexes for (i=0; i
258

OOSspSkylineStatic(nrow, nrow, Ip, Jp, V); }

263

268

273

/* Destructor - deallocates memory reserved for member data */ OOSspSkylineStatic::ÕOSspSkylineStatic() { delete[] diag; delete[] uc; delete[] lr; delete[] iuc; delete[] ilr; } /* CLEAN deallocates memory reserved for member data and initiates them by 0 or NULL */ void OOSspSkylineStatic::clean() { delete[] diag; delete[] uc; delete[] lr;

45

278

delete[] iuc; delete[] ilr; init(); }

283

288

293

298

/* INIT initiates all member data by 0 or NULL */ void OOSspSkylineStatic::init() { row=0; col=0; diag = NULL; iuc = NULL; ilr = NULL; uc = NULL; lr = NULL; llr = 0; luc = 0; }

/* RESIZE(NDIM) removes allocated matrix structures and resizes matrix to the square one with ndim rows and columns */ void OOSspSkylineStatic::resize(TIND ndim) { ASSERT(ndim>=0);

303

clean(); if (ndim>0) { row = col = ndim; diag = new double[ndim]; ASSERT(diag!=NULL); iuc = new TIND[ndim+1]; ASSERT(iuc!=NULL); ilr = new TIND[ndim+1]; ASSERT(ilr!=NULL);

308

for (TIND i=0; i
313

318

} }

323

328

/** RESIZE(NROW,NCOL) removes allocated matrix structures and resizes matrix to the one with nrow rows and ncol columns */ void OOSspSkylineStatic::resize(TIND nrow,TIND ncol) { ASSERT(nrow>=0 && ncol>=0); clean(); if (nrow>0 || ncol>0) { TIND i; row = nrow;

46

col = ncol; diag = new double[MIN(row,col)]; ASSERT(diag!=NULL); iuc = new TIND[col+1]; ASSERT(iuc!=NULL); ilr = new TIND[row+1]; ASSERT(ilr!=NULL);

333

for (i=0; i<MIN(row,col); i++) { diag[i] = 0; } for (i=0; i<=col; i++) iuc[i] = 0; for (i=0; i<=row; i++) ilr[i] = 0;

338

343 } } 348

/* GET returns value from the i-th row and the j-th column */ double OOSspSkylineStatic::get(TIND i,TIND j) { ASSERT(i>=0 && i=0 && j
353

if (i == j ) { // on diagonal return diag[i]; } else if (i < j) { // upper triangle if ((j-i) > columnHeight(j)) // above skyline profile return 0; else return uc[iuc[j+1]-j+i]; } else { // lower triangle if ((i-j) > rowHeight(i)) // above skyline profile return 0; else return lr[ilr[i+1]-i+j]; }

358

363

368

} 373

/* SET sets value in the i-th row and the j-th column to value nval */ void OOSspSkylineStatic::set(TIND i,TIND j,const double& nval) { TIND diff, ci, ri, ni; ASSERT(i>=0 && i=0 && j
378

383

if (i == j ) { // on diagonal diag[i] = nval; } else if (i < j) { // upper triangle diff = (j-i) - columnHeight(j); if (diff > 0 && nval != 0.0) {

47

388

393

398

403

408

413

418

423

428

433

// above skyline profile, need to reallocate columns double *newuc = new double[luc+diff]; // allocate new array for column items for (ci=0; ci 0 && nval != 0.0) { // above skyline profile, need to reallocate rows double *newlr = new double[llr+diff]; // allocate new array for row items for (ri=0; ri
48

438

llr += diff; } else { lr[ilr[i+1]-i+j] = nval; } }

443

448

453

458

463

} /* COPY matrix copy and scale MatS <- alpha*mat */ void OOSspSkylineStatic::copy(OOSspSkylineStatic& mat, const double& alpha ) { if (this==&mat) *this*=alpha; else if (mat.isempty()) clean(); else if (ABS(alpha)<=OOSZERO) {resize(mat.row,mat.col);} else { clean(); TIND i; row = mat.getRows(); col = mat.getCols(); diag = new double[MAX(row,col)]; ASSERT(diag!=NULL); iuc = new TIND[col+1]; ASSERT(iuc!=NULL); ilr = new TIND[row+1]; ASSERT(ilr!=NULL); llr = mat.llr; luc = mat.luc; uc = new double[luc]; lr = new double[llr]; for (i=0; i<MAX(row,col); i++) { diag[i] = alpha*(mat.diag[i]); }

468

for (i=0; i<=col; i++) iuc[i] = mat.iuc[i]; for (i=0; i<=row; i++) ilr[i] = mat.ilr[i]; for (i=0; i
473

478

} }

483

488

/* COPY transposed matrix copy and scale MatS <- alpha*matˆT */ void OOSspSkylineStatic::copyt(OOSspSkylineStatic& mat, const double& alpha) { if (this==&mat) *this*=alpha; else if (mat.isempty()) clean(); else if (ABS(alpha)<=OOSZERO) {resize(mat.row,mat.col);} else { clean();

49

TIND i; row = mat.col; col = mat.row; diag = new double[MAX(row,col)]; ASSERT(diag!=NULL); iuc = new TIND[col+1]; ASSERT(iuc!=NULL); ilr = new TIND[row+1]; ASSERT(ilr!=NULL); llr = mat.luc; luc = mat.llr; uc = new double[luc]; lr = new double[llr];

493

498

for (i=0; i<MAX(row,col); i++) { diag[i] = alpha*(mat.diag[i]); }

503

for (i=0; i<=col; i++) iuc[i] = mat.ilr[i]; for (i=0; i<=row; i++) ilr[i] = mat.iuc[i]; for (i=0; i
508

513 } }

518

523

/* matrix transposition */ void OOSspSkylineStatic::trans() { if (!isempty()) { TIND tmp; TIND* tmpi; double* tmpp; tmp = row; row = col; col = tmp;

528

tmpi = iuc; iuc = ilr; ilr = tmpi; 533 tmp = llr; llr = luc; luc = tmp; 538

tmpp = uc; uc = lr; lr = tmpp; } }

543

50

548

/* SWAP swaps matrices MatS <-> mat */ void OOSspSkylineStatic::swap(OOSspSkylineStatic& mat) { if (this!=&mat) { OOSspSkylineStatic temp; temp=*this; *this=mat; mat=temp; } }

553

558

563

568

/** returns true if the matrix is square and has a symmetric profile */ bool OOSspSkylineStatic::symmetricProfile() { TIND i; if (row != col) return false; if (llr != luc) return false; for (i=0; i<=row; i++) if (iuc[i] != ilr[i]) return false; return true ; } /** return the profile of Skyline matrix given by vectors I,J,V and permutation vector p*/ static int profile(TIND nrow, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V, OOSindexArray& p) { ASSERT(I.getDim() == J.getDim()); TIND N = I.getDim(); TIND i;

573 int* ch = new int[nrow]; int* rh = new int[nrow];

578

583

588

593

for (int i=0; i
51

// diagonal } else if (V(i) != 0.0) { if (Ii<Ji) { // upper triangle h = Ji-Ii; if (h>ch[Ji]) ch[Ji] = h; } else { // lower triangle h = Ii-Ji; if (h>rh[Ii]) rh[Ii] = h; } }

598

603

608 }

int sum = 0; for (i=0; i
613

618 }

623

628

633

638

643

648

/* MATRIX OPERATIONS */ /* ================= */ /* OPERATOR= (MatS=mat) assigns member data from the matrix mat to MatS and initiates mat to empty vector */ void OOSspSkylineStatic::operator=(OOSspSkylineStatic& mat) { clean(); // Remove allocated arrays of *this // Shallow copy: *this <- mat row=mat.row; col=mat.col; binary=mat.binary; diag = mat.diag; uc = mat.uc; lr = mat.lr; iuc = mat.iuc; ilr = mat.ilr; luc = mat.luc; llr = mat.llr; // Reset mat to empty matrix mat.clean(); }; /* OPERATOR*=(alpha) matrix-scalar multiplication MatS<-alpha*MatS */ void OOSspSkylineStatic::operator*=(const double& alpha) { TIND i;

52

if (alpha!=1.0) { if (alpha!=1.0) { for (i=0; i
653

658

663

} 668

673

/* OPERATOR/=(alpha) matrix-scalar dividing MatS<-MatS/alpha */ void OOSspSkylineStatic::operator/=(const double& alpha) { TIND i; ASSERT(ABS(alpha)>OOSZERO); if (alpha!=1.0) { for (i=0; i
678

683

688

693

698

}

/* OPERATOR- (unary) changes the signs of the vector entries MatS <- (MatS) */ void OOSspSkylineStatic::operator-() { TIND i; for (i=0; i
53

703

lr[i] = -lr[i]; } }

708

713

718

/* MATRIX NORMS */ /* ============ */ /* NORM_1 returns the column norm */ double OOSspSkylineStatic::norm_1() { TIND i,j,ind; if (isempty()) return 0.0; else { OOSvectorD y(col,0.); for (j=0; j
723

728

733 return y(y.max()); } } 738

743

748

753

/* NORM_INF returns the row norm */ double OOSspSkylineStatic::norm_inf() { TIND i,j,ind; if (isempty()) return 0.0; else { OOSvectorD y(col,0.); for (i=0; i
54

y(i) += ABS(get(i,j)); 758

} } return y(y.max()); }

763

768

} /* NORM_F returns the Frobenius norm */ double OOSspSkylineStatic::norm_F() { if (isempty()) return 0.0; else { TIND i; double s;

773 s=0.; // sum upper triangle for (i=0; i
778

783

788

return sqrt(s); } }

793

798

/* NORM matrix norm ||A||_i i==0 ... ||A||_F (Frobenius norm) i==1 ... ||A||_1 (row norm) i==2 ... ||A||_2 (spectral norm) i==3 ... ||A||_inf (column norm) */ double OOSspSkylineStatic::norm(const char i) { ASSERT(i>=0 && i<=3); switch (i) { case 0: return norm_F(); case 1: return norm_1(); case 2: return norm_2(); case 3: return norm_inf(); default: return OOSINF; }

803

808 }

55

813

818

/** NORM returns a norm of the diagonal, type defines type of the norm type==1 -> l_1 type==2 -> l_2 type==3 -> l_inf <default> */ double OOSspSkylineStatic::diagNorm(const char type) { ASSERT(type>=1 && type<=3); TIND i, m=getRows(), n=getCols(), d; double res=0, w, temp;

823

d=MIN(m,n); switch (type) { case 1: for(i=d-1;i>=0;i--) {w=diag[i]; if (w<0.0) res-=w; else res+=w;}; break; case 2: for(i=d-1;i>=0;i--) {w=diag[i]; res+=w*w;}; break; case 3: for(i=d-1;i>=0;i--) {w=diag[i]; if (res<(temp=ABS(w))) res=temp;}; break; default: res=OOSINF; } return res;

828

833 }

838

843

848

/* STREAM OPERATORS OVERLOADING */ /* ============================ */ /* Overloaded output stream operator */ ostream& operator<<(ostream& os, OOSspSkylineStatic& mat) { TIND i, j, lrow=mat.getRows(), lcol=mat.getCols(); double v; if (mat.isempty()) { os << endl << "Empty matrix!" << endl; ASSERT(os.good()); return os; } os << endl << "Matrix type (" << lrow << "," << lcol << ")" << endl; ASSERT(os.good()); for (i=0; i
853

858

}

56

863

/* Overloaded input file stream operator */ istream& operator>>(istream& is,OOSspSkylineStatic& mat) { TIND i, j, k, lrow, lcol, lnnz; double v;

868 if (mat.isBinary()) { is.read((char*)&lrow,sizeof(TIND)); is.read((char*)&lcol,sizeof(TIND)); is.read((char*)&lnnz,sizeof(TIND)); } else { is >> lrow >> lcol >> lnnz; } ASSERT(is.good()&&lrow>=0&&lcol>=0&&lnnz>=0);

873

878 OOSindexArray I(lnnz), J(lnnz); OOSvectorD V(lnnz); if (lrow>0 && lcol>0 && lnnz>0) { if (mat.isBinary()) { ASSERT(is.good()); is.read((char*)I.getInd(),lnnz*sizeof(TIND)); is.read((char*)J.getInd(),lnnz*sizeof(TIND)); is.read((char*)V.getVal(),lnnz*sizeof(double)); } else for (k=0; k> i >> j >> v; I(k)=i; J(k)=j; V(k)=v; }

883

888

893

898

mat.spconvert(I,J,V,lrow,lcol); } return is; }

903

908

913

/* Overloaded output file stream operator */ ofstream& operator<<(ofstream& os,OOSspSkylineStatic& mat) { TIND i, j, k, c, n, lrow=mat.getRows(), lcol=mat.getCols(), lnnz=mat.skylineProfile()+MIN(lrow,lcol); double v; if (mat.isBinary()) { os.write((const char*)&lrow,sizeof(TIND)); os.write((const char*)&lcol,sizeof(TIND)); os.write((const char*)&lnnz,sizeof(TIND));

57

} else os << lrow << ’\t’ << lcol << ’\t’ << lnnz << endl; ASSERT(os.good());

918

OOSindexArray I(lnnz), J(lnnz); OOSvectorD V(lnnz); 923

c=0; TIND rh, ch; // diagonal elements for (i=0; i<MIN(lrow, lcol); i++) { I(c) = i; J(c) = i; V(c) = mat.get(i,i); c++; } // upper triangle elements for (j=0; j=(j-ch);i++) { I(c) = i; J(c) = j; V(c) = mat.get(i,j); c++; } } // lower triangle elements for (i=0; i=(i-rh);j++) { I(c) = i; J(c) = j; V(c) = mat.get(i, j); c++; } } if (mat.isBinary()) { os.write((const char*)I.getInd(),lnnz*sizeof(TIND)); os.write((const char*)J.getInd(),lnnz*sizeof(TIND)); os.write((const char*)V.getVal(),lnnz*sizeof(double)); ASSERT(os.good()); } else { for (k=0; k
928

933

938

943

948

953

958

963

} }

Vypis 2: Soubor OOSmatricesSkylineStatic.cpp ´

58

#include "OOSgraph.h" namespace oosol { 5

10

OOSgraph::OOSgraph(OOSmatrix &mat) { ASSERT((n = mat.getRows()) == mat.getCols()); TIND i,j; TIND adjpos=0; iadj = new TIND[n+1]; ASSERT(iadj!=NULL); for (i=0; i<=n; i++) iadj[i] = 0;

15

ladj = 0; // find number of non-zero elements for (j=0; j
20

25

30

35

40

} 45

OOSgraph::OOSgraph(TIND nrow, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V) { TIND i, N=I.getDim(); ASSERT((N == J.getDim()) && (J.getDim() == V.getDim()));

50 n = nrow; iadj = new TIND[n+1]; TIND* adjlen = new TIND[n];

59

ladj = 0; 55 for (i=0; i
TIND ii, ji; double vi; for (i=0; i= 0) && (ii <= nrow) && (ji >= 0) && (ji <= nrow));

65

if (vi != 0.0) { if (ii < ji) { // count only nodes in upper triangle ladj+=2; // every edge will be stored twice adjlen[ii]++; adjlen[ji]++; } }

70

75

} iadj[0] = 0; for (i=1; i<=n; i++) iadj[i] = iadj[i-1]+adjlen[i-1]; 80 adj = new TIND[ladj]; for (i=0; i
95

100

} delete[] adjlen; }

105

void OOSgraph::sortNodeByDegree(TIND* array, TIND start, TIND end) { int i, j, key, v; for (j=start+1; j<end; j++) { v = array[j];

60

key = nodeDegree(v); for(i=j-1; (i>=start) && (nodeDegree(array[i])>key); i--) { array[i+1] = array[i]; } array[i+1] = v;

110

} } 115

120

125

130

135

140

145

150

155

void OOSgraph::rootLevelStructure( TIND root, // root node TIND* mask, // masking vector TIND* rls, // root level structure TIND* irls, // index to root level structure TIND &rlslen, // length of root level structure TIND &depth, // depth of root level structure TIND &maxwidth) // maximal width of root level structure) { TIND i, j, parent, nb; TIND curwidth = 1; // width of current level TIND lvls, lvle; // start and end of level depth = 1; rls[0] = root; irls[0] = 0; irls[1] = 1; lvle = 0; rlslen = 1; mask[root] = 1; maxwidth = 1; while (1) { lvls = lvle; lvle = lvls+curwidth; curwidth = 0; // for every node in previous depth level for (i=lvls; i
160

depth++; irls[depth] = rlslen;

61

maxwidth = MAX(curwidth, maxwidth); } 165

// reset mask for nodes in root level structure for (i=0; i
170

175

180

185

190

195

200

205

void OOSgraph::pseudoDiameter(TIND* mask, TIND &s, TIND* rls, TIND* irls, TIND &rlslen) { TIND e; // ending node TIND sd; TIND szd; TIND d; TIND i,j; TIND* q; // shrunken last level TIND qdepth; TIND qwidth; TIND minwidth; TIND node; TIND depth; sd = n; // first guess for starting node // find node with smallest degree for (i=0; i
210

while (news) { // sort last level by node degree sortNodeByDegree(rls, irls[depth-1], irls[depth]); // shrink last level - create list with one node of each degree q = new TIND[irls[depth]-irls[depth-1]];

62

215

TIND ql = 1; TIND degree; q[0] = rls[irls[depth-1]]; degree = nodeDegree(q[0]); for (i=irls[depth-1]+1; i degree) { degree = d; q[ql++] = node; } }

220

225

minwidth = n+1; 230

news = false; // generate root level structure for each node in shrunken level for (i=0; i depth) { // level structure with greater depth, new starting node s = node; depth = qdepth; news = true; break; } else { // new ending node e = node; minwidth = qwidth; } } } delete[] q;

235

240

245

250

} if (e != node) { // generate root level structure from end node rootLevelStructure(e, mask, rls, irls, rlslen, qdepth, qwidth); } // store distances from end node for each node // nodes that are not part of this component will have // zero for (i=0; i<depth; i++) { for (j=irls[i]; j
255

260

} 265

void OOSgraph::number(int W1, int W2, TIND s, TIND* rls, TIND rlslen, int *status, TIND *q, int *priority, int &labeled) { TIND i, j; TIND node;

63

TIND qlen; // lenght of queue 270

275

// initialize node status and priority for (i=0; i
280

q[0] = s; qlen = 1; int maxp; TIND maxpi;

285

290

295

while (qlen > 0) { // find node with max priority maxp = priority[q[0]]; maxpi = 0; for (i=1; i maxp) { maxp = priority[node]; maxpi = i; } } node = q[maxpi];

300

// delete node from queue q[maxpi] = q[qlen-1]; qlen--; TIND nb, anb;

305

310

// preactive node, search neighbors if (status[node] == PREACTIVE) { for (i=iadj[node]; i
315

} } 320 // assign new node number, status is now postactive

64

labeled++; status[node] = labeled; 325

// search preactive neighbors for (i=iadj[node]; i
330

// search neighbors of active neighbor for (j=iadj[nb]; j
335

340

} 345

} } }

350 OOSindexArray* OOSgraph::sloan() { return sloan(1, 2); } 355

360

365

OOSindexArray* OOSgraph::sloan(int W1, int W2) { TIND s; // start node of pseudo-diamater TIND i; TIND labeled; // number of nodes already labeled TIND *rls = new TIND[n+1]; // root level structure items TIND *irls = new TIND[n+1]; // index do root level structure TIND *q = new TIND[n]; // priority queue int *status = new int[n]; // status of node (or new label) int *priority = new int[n]; // node priority TIND componentLen; // number of nodes in this graph component OOSindexArray *perm; // result permutation vector // initialize status vector for (i=0; i
370 labeled = 0; // while there is unlabelled component of graph while (labeled < n) { pseudoDiameter(status, s, rls, irls, componentLen); 375

65

number(W1, W2, s, rls, componentLen, status, q, priority, labeled); }

380

perm = new OOSindexArray(n); for (i=0; iset(status[i]-1, i);

385

delete[] delete[] delete[] delete[] delete[]

priority; status; rls; irls; q;

return perm; 390

395

} OOSgraph::ÕOSgraph() { delete[] adj; delete[] iadj; } }

Vypis 3: Soubor OOSgraph.cpp ´ 2

7

12

/********************************************************************/ /* */ /* OOSol2 2008 */ /* Dept. of Applied Mathematics, */ /* Technical University in Ostrava, Czech Rep. */ /* */ /* Filename: OOSgraph.h */ /* Created by Martin Stachon, 2008/05 */ /* */ /* List of updates: */ /* Date Author Description */ /* */ /********************************************************************/ #ifndef OOSGRAPH_H_ #define OOSGRAPH_H_

17 #include "OOSobjects.h" #include "OOSmatrices.h" namespace oosol { 22 class OOSgraph { /** adjacency list */ TIND* adj; 27 /** length of adjacency list */

66

TIND ladj; /** index of adjacency list */ TIND* iadj;

32

/** number of nodes in graph */ TIND n; 37

/** returns the degree of node i */ TIND nodeDegree(TIND i) { return iadj[i+1]-iadj[i]; }

42

/** Sloan algorithm node status */ //POSTACTIVE > 0 // label already assigned #define ACTIVE (0) // adjacement to a postactive node #define PREACTIVE (-1) // adjacement to an active node #define INACTIVE (-2) // none of above

47 void pseudoDiameter(TIND* mask, TIND &s, TIND* rls, TIND *irls, TIND & rlslen); void sortNodeByDegree(TIND* array, TIND start, TIND end); 52

void rootLevelStructure(TIND root, TIND* mask, TIND* rls, TIND* irls, TIND &rlslen, TIND &depth, TIND &width ); void number(int W1, int W2, TIND s, TIND* rls, TIND rlslen, int *status, TIND *q, int *priority, int &labeled); public:

57 /** Constructor - constructs graph from full rectangular matrix */ OOSgraph(OOSmatrix &mat); /** Constructor - constructs graph from rectangular matrix with nonzeros given by vectors of indices I, J */ OOSgraph(TIND nrow, OOSindexArray& I, OOSindexArray& J, OOSvectorD& V);

62

ÕOSgraph(); OOSindexArray* sloan(); 67 OOSindexArray* sloan(int W1, int W2);

}; 72 } #endif /*OOSGRAPH_H_*/

Vypis 4: Soubor OOSgraph.h ´

Objektově orientovaná implementace skyline systému

Recommend Documents