Knihovna pro práci s maticemi komplexních čísel ve volitelné přesnosti

2006/56– 8.12.2006

Knihovna pro práci s maticemi komplexních čísel ve volitelné přesnosti Ing. Václav Pfeifer Ing. Miroslav Balík, Ph.D. Vysoké učení technické v Brně, Fakulta elektrotechniky a komunikačních technologií, Ústav telekomunikací, Purkyňova 118, 612 00 Brno, Česká republika email: [email protected], [email protected]

Tento článek seznamuje čtenáře s možnostmi a využitím knihovny c_matice určené pro práci s maticemi komplexních čísel ve volitelné přesnosti. Knihovna je navržena a naprogramována v jazyce C++.

1. Úvod V dnešní době výpočetní schopnost hardwaru neustále vzrůstá. Ať už jde o klasické domácí počítače, nebo o nové signálové procesory, či cokoliv jiného. S tímto neustálým zvyšováním výpočetního výkonu máme možnost, nebo také někdy i potřebu, řešit stávající problémy neustále mnohem rychleji. V některých případech nám ale nemusí jít přímo o rychlost výpočtu, ale třeba o provedení některých operacích v mnohem větší přesnosti, která v dřívější době nebyla možná, právě z důvodu nedostačujícího hardwarového výkonu stávajících zařízení. V dnešní době je standardní zobrazovat reálná čísla (typ double nebo float v programovacím jazyku C++) na 8 desetinných míst. Pro většinu standardních aplikací a výpočtů je to dostačující, ale jsou aplikace, v kterých se pracuje v mnohem vyšších resp. nižších řádech a pro tyhle aplikace je vhodné pracovat ne v jednoduché přesnosti, ale v tzv. volitelné přesnosti, kdy si sami můžeme zvolit na kolik desetinných míst chceme číslo zobrazit, nebo s jakou přesností bude s číslem zacházeno. Základní funkce pro práci s čísly ve volitelné přesnosti, které jsou použity při konstrukcí naších funkcí, jsou obsaženy v knihovně APFLOAT.

2. Knihovna APFLOAT APFLOAT je knihovna naprogramovaná v C++, která efektivně zpracovává všechny základní matematické operace, jako jsou sčítání, odečítání, násobení atd. ve volitelné přesnosti. Knihovna taktéž obsahuje spoustu dalších vhodných funkcí, které mohou programátorovi ulehčit práci. Knihovna APFLOAT dokáže nejen efektivně zpracovávat reálná čísla, ale umožňuje přímo pracovat také s komplexní algebrou. Výpočetní algoritmy jsou efektivně optimalizovány a dokáží bez problémů pracovat s přesností až na milion desetinných míst. Násobení je prováděno pomocí Rychlé Fourierovy Transformace, resp. FFT a upraveno aby nedocházelo k problémům se zaokrouhlováním. Knihovna je napsána tak, aby byla co nejsnadněji aplikována, ale taktéž obsahuje assemblerovskou optimalizaci v kritických částech pro co nejvyšší možný výkon. Knihovna APFLOAT je volně šiřitelný Freeware pro nekomerční použití. Dostupné z http://www.apfloat.org.

3. Struktura Třídy Matice

56-1

2006/56– 8.12.2006

Obr. 3.1: Class diagram třídy matice Obrázek 3.1. nám ukazuje základní strukturu třídy c_matice. Kromě položky hodnoty jsou všechny proměnné a funkce statickými. Proměnná hodnoty je ukazatel, který odkazuje na položky hodnot definovaných matic. Hodnoty řádku a sloupců určují rozměry zadané matice. V našem případě jsou předávány pomocí konstruktoru. Třída c_matice ještě obsahuje další funkce pro práci jak s maticemi, tak pro alokaci a práci s dynamickým polem. Následující funkce budou rozebrány podrobněji níže. V prvním bloku jsou zobrazeny všechny privátní proměnné a funkce. Znaménko minus před názvem funkce značí že se jedná právě o metodu privátní. V dalším bloku jsou konstruktory, které se obecně řadí mezi privátní, i když se s nimi tak přímo nezachází. V předposledním bloku jsou klasické veřejné metody, které jsou označeny právě znaménkem plus a v posledním bloku jsou tzv. spřátelené funkce, které mají některé výhodné vlastnosti, jako přístup do mezi privátní proměnné a funkce. V diagramu nejsou zobrazeny přetížené operátory pro základní matematické operace a porovnávací výrazy a taktéž nejsou zobrazeny funkce pro práci s řídkými maticemi, které mají v podstatě úplně stejnou deklaraci, jako funkce klasické. 56-2

2006/56– 8.12.2006 Nakonec ještě podotýkám, že společně s třídou matice spolupracuje ještě struktura LUmatice, která obsahuje 2 prvky typu c_matice. Tahle struktura je použita pouze u funkce pro LU rozklad, ale jelikož je spjata s třídou c_matice, je vhodné ji zde připomenout.

3.1. KONSTRUKTORY Implicitní konstruktor V našem případě jsme implicitní konstruktor přepsali tak, aby po deklaraci přiřadil do privátních proměnných nulové hodnoty a příslušný ukazatel na naše hodnoty na NULL. c_matice::c_matice();

Používání tohoto konstruktoru je však v našem případě zbytečné a bezúčelové, protože nejsou definované žádné veřejné metody, které by definovaly rozměry matice. I když není problém sestrojit takovou metodu, je pohodlnější a rychlejší použití konstruktoru s parametry, které budou obsahovat právě rozměry požadované matice Konstruktor pro čtvercovou matici Jde o konstruktor, který má jen 1 parametr, který nám zároveň definuje počet řádků i sloupců matice. Jedná se tedy o čtvercovou matici řádu n, který je předán parametrem. c_matice::c_matice(const int n);

Konstruktor pro obecnou matici Následující konstruktor obsahuje v deklaraci 2 proměnné, které nám definují rozměr požadované matice, počet řádku a počet sloupců. Zároveň je taktéž provedena alokace paměti pro zadaní rozměr matice. c_matice::c_matice(const int x, const int y);

Konstruktor pro jednotkovou a nulovou matici Následující konstruktor slouží pro vytvoření čtvercové matice řádu n. Podle posledního parametru je určen typ matice. J značí matici jednotkovou a N značí matici nulovou. V případě že zadáme v parametru cokoliv jiného, je vytvořena klasická čtvercová matice zvoleného řádu. c_matice::c_matice(const int x, const char typ);

Kopírovací konstruktor Kopírovací konstruktor je také jedním z implicitních konstruktorů, které se vytvářejí automaticky, když není uživatelsky definován. Jedná se o konstruktor, pomocí kterého můžeme přímo definovat celou matici, pomocí matice jiné. Výsledná matice bude mít stejné rozměry včetně všech hodnot. c_matice::c_matice(c_matice &a);

3.2. Soukromé Metody Třída matice obsahuje celkem 9 privátních metod. Všechny z uvedených metod jsou konstrukční funkce pro složitější matematické operace s maticemi, jako je například Gaussova eliminace, nebo třeba inverze čtvercové matice.

56-3

2006/56– 8.12.2006 Prvních 7 funkcí jsou metody pro práci s maticemi, které jsou vyjádřeny klasickým způsobem. To znamená, že do pole je vždy na danou souřadnici uložena konkrétní hodnota prvku. Dalším možným způsobem je zapisovat pouze nenulové prvky a pamatovat si jejich souřadnice. Pro klasickou matici je tenhle zápis krajně nepohodlný a výpočetně zdlouhavý, ale v případě, že pracujeme s matici řídkou, je tento zápis nejvýše vhodný. Poslední 2 funkce budou tedy interní metody pro práci s řídkými maticemi.

3.3. Veřejné Metody Jako veřejné metody jsem vytvořil všechny základní operace, které se dají s maticemi provádět a některé jiné, které mají uplatnění v našem elektrotechnickém oboru. Z možnosti budoucího využití, nebo modifikace naší třídy c_matice jsem přetížil všechny základní matematické operátory, jako je například operátor součtu, rozdílu, ale také porovnávací operátory a operátor přiřazení. Základní matematické operátory Jedná se o jednoduché operace, které patří k nejčastěji používaným jak při práci s maticemi, tak při práci s jednoduchými čísly. Jedná se o sčítání, odčítání a násobení. Další operátory pro základní matematické operace Pro kompletnost jsem ještě přetížil operátory, které slučují 2 typy operátorů. Jedním typem operátorů jsou základní matematické operátory jako je sčítání. Druhým operátorem je operátor přiřazení. Funkce těchto operátorů je založena na provedení základní matematické operace, ale s tím že výsledek této operace je uložen právě na levé straně výrazu. c_matice &operator+=(c_matice &a); c_matice &operator-=(c_matice &a); c_matice &operator*=(c_matice &a);

Operátor přiřazení matice Stejně jak v C++ funguje přiřazení jedné proměnné do druhé pro základní datové typy, můžeme přetížit operátor i pro naše vlastní objekty naší třídy c_matice. Do proměnné na levé straně jsou zkopírovaná data z proměnné na pravé straně. Samozřejmě že se musí jednat o stejný datový typ, v našem případě třídy c_matice. c_matice &operator=(c_matice &a);

V případě že nejsou matice stejného řádu, je matice na levé straně realokována podle rozměrů matice na straně druhé. V případě že jsou matice totožné a to včetně řádu a hodnot, je vrácena naše původní matice. Logické operátory Logické operátory jsou takové, které vyjadřují určitou matematickou skutečnost, jako je třeba rovnost, nerovnost apod. V třídě c_matice jsem přetížil konkrétně 2 nejpoužívanější operátory a to operátor porovnávání a operátor nerovnosti. bool operator==(c_matice &a); bool operator!=(c_matice &a);

56-4

2006/56– 8.12.2006 Jako u ostatních datových typů, tak i zde platí, pokud je výraz na levé straně totožný s výrazem na straně pravé, tak je vrácena hodnota TRUE. Jinak je vrácena hodnota FALSE. Pro operátor nerovnosti je tomu přesně naopak. Takže v případě rovnosti je vrácena hodnota FALSE a v případě nerovnosti hodnota TRUE. Matematicky můžeme také hodnotě TRUE přiřadit hodnotu 1 a hodnotě FALSE 0. Metoda gauss Gaussova eliminace matice je proces, kdy se snažíme dostat matici na tzv. schodovitý tvar. Právě díky tomuto procesu se Gaussova eliminace používá v oblastech výpočtu kořenů soustavy rovnic, ať už lineárních, nebo diferenciálních. Gaussova eliminace je také základ pro náš algoritmus výpočtu inverze matice a determinantu matice. c_matice &gauss();

Funkce gauss provádí Gaussovu eliminaci v tom smyslu, že pomocí výše popsaných způsobů provádí transformaci zadané matice na dolní trojúhelníkovou matici. Výsledná upravená matice je pak vrácena jako návratový argument matice. V případě, že na diagonále se vyskytuje nulový prvek, provede funkce gauss sečtení s nejbližšími řádky a to do doby, než se na diagonálním prvku nevyskytuje nula. V případě, že po sečtení se všemi řádky se na některém diagonálním prvku stále vyskytuje nula, vrátí funkce chybové hlášení. Metoda det Funkce det vrací determinant čtvercové matice. Determinant je určen za pomocí gaussovy eliminace, kdy převedeme matici na tzv. schodovitý tvar, resp. na dolní trojúhelníkovou matici resp. Na horní trojúhelníkovou matici. Z lineární algebry víme, že determinant takto upravené matice lze přímo určit tak, že vynásobíme spolu řádky na diagonále a výsledek tohoto součinu je přímo roven determinantu. apcomplex det(); apcomplex det(const int prec);

Obr. 3.1: Princip Gaussovy eliminace pro výpočet determinantů vyšších řádů Na výpočet determinantu čtvercové matice jsem vytvořil 2 přetížené funkce. Rozdíl je v tom, že první z funkcí vypíše determinant v maximální možné přesnosti a druhá funkce vypíše determinant v námi zvolené přesnosti. Návratové hodnoty jsou ale u obou funkcí stejné, rozdíl je tedy pouze v zobrazení čísla.. Jelikož lze determinant určit pouze ze čtvercové matice, ošetřil jsem funkci pro tuto situaci. V případě jiné zvolené matice než čtvercové vrátí funkce chybové hlášení. Jedná li se taktéž o singulární matici, vrátí funkce opět chybové hlášení.

Metoda inverze

56-5

2006/56– 8.12.2006 Inverze funkce je důležitým prostředkem při práci s maticemi. Pomocí inverze můžeme třeba určit podíl matic, který je definován jako součin první matice s inverzí matice druhé. Podíl matic ale nemá tak veliké použití, proto ho nebudu přímo uvádět jako funkci. Velké využití inverze je i v oblasti elektrotechniky pro analýzu obvodu, nebo v teorii elektromagnetického pole. Asi nejznámějším způsobem určení inverze je pomocí adjungované matice, kde každý prvek adjungované matice je subdeterminant matice řádů o stupeň menší. Výsledek pak dostaneme vydělením odpovídajícím determinantem. Avšak pro naši práci je tahle metoda příliš pomalá, protože musíme provádět příliš mnoho operací s maticemi a jelikož pracujeme ve volitelné řádové čárce, což zabere spoustu času. Z tohoto důvodu je vhodné inverzi určit jiným způsobem. Jako vhodný nástroj se ukázalo sepsat si k naší matici jednotkovou matici stejného řádu a pomocí stejných matematických operací převést původní matici na jednotkovou a z původní jednotkové matice je naše inverzní matice. c_matice &inverze();

Obr. 3.2: Princip určení inverzní matice(a je původní matice, b je inverzní matice) Obr. 3.2. ukazuje praktický výpočet inverze matice. Z našeho hlediska jde tedy pouze o vytvoření jednotkové matice a provedení na ní všech operací, které provádíme na naší matici v případě Gaussovy eliminace. Vstupním parametrem musí být pouze čtvercová matice, z žádné jiné nelze vytvořit inverzi. Další podmínkou je, že vstupní matice nesmí být singulární, tedy její determinant nesmí být roven nule. V případě splnění těchto podmínek můžeme určit inverzi matice. Metoda Transponuj Jak už vypovídá název slouží tato funkce k transpozici matice. Transpozice sice není tak často používaná funkce, jako jsou třeba inverze či determinant, ale pro analýzu obvodů, nebo v teorii kódování je taktéž nezbytnou operací. Principiálně jde pouze o prohození řádků a sloupců. friend matice transponuj(matice &A);

Metoda eigenvalues Následující metoda slouží k výpočtu vlastních čísel čtvercové matice libovolného řádu. Principielně je výpočet vlastních čísel velice snadný. Ale z důvodu že se v matici vyskytují parametry je algoritmizace podle obecné definice velice náročná. Z tohoto důvodu byla vytvořena spousta algoritmů, které nám určují vlastní čísla postupnou iterací dle určitého algoritmu. Mezi nejznámější algoritmy patří metoda QR a Jacobiova diagonalizace. První z metod je pro větší matice výpočetněji jednodušší, ale je o to problematičtější v algoritmizaci. Proto jsem se rozhodl pro použití Jacobiovy metody, která je co do algoritmizace jednodušší. Jelikož nemůžeme použít klasickou Gaussovu eliminaci, protože mění vlastní čísla, musíme najít takovou metodu, kterou bychom upravili matici na diagonální tvar a prvky na diagonále by byly naše vlastní čísla.

56-6

2006/56– 8.12.2006 Principem Jacobiovy metody je ortogonální transformace naší matice. Jde o to že můžeme vynásobit naši matici zprava i zleva takovými maticemi, jejichž vzájemný součin dává jedničku, tudíž se jedná o operaci která nemůže měnit vlastní čísla naší matice. Jde tedy o to určit takové matice, které nám naši požadovanou matici po vzájemném vynásobení upraví na matici diagonální.

(3.1.) Z 3.1. je vidět že musíme určit matici Q, z které pak určíme matici transponovanou a vzájemným roznásobením pak získáme požadovaný tvar naší matice. Matici Q určíme postupně pomocí posloupnosti ortogonálních matic

.

(3.2.) Jacobiova metoda má tu výhodu, že nemusíme přímo hledat ortogonální matici , ale můžeme přímo určit matici B. Výpočet provádíme v několika krocích. V prvním kroku zjišťuji největší nediagonální prvek a podle něho určím indexy p a q. V druhém kroku si určím úhel θ a z toho určím další parametry t , c , s. V posledním kroku určím podle konkrétních vzorců přímo prvky matice B. Jak již jsem se zmiňoval, tak jde o iterační metodu a tudíž naše matice B ještě zdaleka není v diagonálním tvaru s vlastními čísli na hlavní diagonále. Ale prvky mimo hlavní diagonálu se nám zmenšily. To znamená že výše uvedený algoritmus budeme nyní aplikovat na nově vytvořenou matici B. Tímto postupným upravováním zajistíme že prvky mimo hlavní diagonálu nám budou postupně konvergovat k nule. Iterace provádíme do okamžiku, kdy prvky mimo hlavní diagonálu jsou menší než nějaká námi zvolená mez. c_matice ownnumbers(apfloat prec);

3.4. Veřejné metody pro práci řídkými maticemi V našem případě nebudeme používat speciální metody pro základní matematické operace s řídkými maticemi, jelikož by to nemělo očekávaný efekt. Při práci s řídkýma maticemi budeme pro základní matematické operace využívat základní funkce naprogramované pro obecné matice, avšak pro funkce jako Gaussova eliminace, či funkce postavené na principu gaussovy eliminace, jako je například inverze, budeme využívat zvlášť naprogramovaných metod. Principielně jsou tyhle metody téměř totožné s metodami pro práci s obecnými maticemi. Rozdíl spočívá, že při operacích s řádky a sloupci jsou brány v úvahu jen nenulové prvky, protože nulové prvky se při naších operacích neprojevují. Dalším důležitým rozdílem je použití tzv. Markowitzova kritéria, které nám ještě výkonnostně urychlý požadované výpočty. c_matice &gauss_sparsematrix(); apcomplex det_sparsematrix(); apcomplex det_sparsematrix(const int prec); c_matice &inverze_sparsematrix();

Jak je vidět jsou deklarovány a definovány 4 základní metody pro práci s řídkými maticemi. Všechny z těchto metod mají stejné funkce i použití, jako tomu bylo u klasických metod v kapitole 3.2. To znamená, že můžeme tyhle metody použít při práci s obecnými maticemi, 56-7

2006/56– 8.12.2006 ovšem musíme počítat že v takovém případě nedojde k urychlení výpočtu, ale naopak se výpočet ještě prodlouží, z důvodu aplikace markowitzova algoritmu, který v případě obecné matice nepřináší žádné urychlení výpočtu. Pro připomenutí ještě uvedu, že stejně jako klasické veřejné metody využívají pomocných privátních metod, tak i pro následující metody jsou používány privátní metody, avšak upravené pro práci s řídkými maticemi.

3.5. Vstupní a výstupní funkce Vstupní a výstupní funkce nám určují formát vstupu a formát výstupu naší matice. Jelikož se předpokládá práce s maticemi vyšších řádů, vytvořil jsem funkce taktéž pro vstup i výstup do souboru. Z praktického hlediska jsem provedl přetížení všech vstupních a výstupních funkcí pro nejčastější případy, které se můžou naskytnout. Pro zadávání hodnot z klávesnice jsem taktéž přetížil operátor vstupu, což se mi jeví velice vhodné, hlavně v případě pokud bude někdo pracovat přímo s mou třídou. Taktéž jsem provedl přetížení výstupního streamového operátoru pro jednoduchý výpis na obrazovku.

3.6. Funkce Funkce LUfakt Funkce provádí tzv. LU faktorizaci, nebo také rozklad, matice. Principielně jde o to, že rozdělíme naši požadovanou matici na součin dvou trojúhelníkových matic, kde matice L je dolní trojúhelníková matice a matice U je horní trojúhelníková matice.

Obr. 3.3: LU rozklad čtvercové matice řádu 4

Jak je vidět z obr. 3.3. jsou na diagonále u matice L jedničky. Z důvodu pohodlnějšího výpočtu, hlavně v případě že provádíme výpočet ručně, dosadíme za diagonální prvky jedničky a dopočítáme zbytek. Obecně si můžeme dosadit cokoliv, ale z možnosti pohodlnějších pozdějších operací je výhodné zvolit právě jedničky. Funkce je prováděna ve dvou krocích. V prvním kroku jsou načteni na diagonálu matice L jedničky. V druhém kroku se střídavě provádí výpočet daného řádku matice U a daného sloupce matice L. Tímto způsobem dostaneme výsledné matice L a U. friend c_LUmatice LUfakt(c_LUmatice &b, c_matice &a); struct c_LUmatice{ c_matice L; c_matice U; };

Pro pohodlnější práci jsem se rozhodl vytvořit strukturu LUmatice, která v sobě obsahuje právě 2 prvky a to matici L a matici U. Oba z těchto prvků jsou datového typu c_matice. Z důvodu vzájemné asociace třídy c_matice a struktury LUmatice, jsem musel před samotnou třídou c_matice připsat deklaraci struktury LUmatice.

56-8

2006/56– 8.12.2006 LU rozklad je základním prostředkem k řešení soustavy lineárních rovnic, kdy výpočetní náročnost je stejná jako u Gaussovy eliminace tj. být výpočet inverzní matice.

. Další aplikací LU rozkladu může

4. Ukázka použití knihovny 4.1. Délka výpočtů v závislosti na řádu a přesnosti matice Pro základní informace o rychlosti knihovny jsem vytvořil testovací program, který počítal rychlosti výpočtů určitých metod. Výsledky jsem pak vynesl do grafu a zakreslil pro různé přesnosti. VÝPOČET INVERZNÍ HILBERTOVY MATICE A JEJÍ SOUČIN S MATICÍ PUVODNÍ Matematický zápis: X = Hilb(A) * inverze(Hilb(A)) ; Kde matice X je maticí diagonální - jednotkovou

Obr. 4.1: Graf závislosti rychlosti výpočtu na řádu matice Obr. 4.1. nám ukazuje rychlost výše uvedeného výpočtu v závislosti na řádu hilbertovy matice. Jak je z grafu vidět, tak přesnost zhruba o řád vyšší nijak extrémně výpočet nezpomaluje. Až zhruba o 2-3 řády výše dochází k prudkému zvýšení doby výpočtu, které je způsobeno nutností uchovávat data v mezipaměťovém bufferu na pevném disku. Jelikož se ne každý podíl dá vyjádřit na konečný počet míst, můžou se na mimodiagonálních pozicích objevit nenulové prvky. U programu MATLAB, který pracuje na maximálně 16 desetinných míst, může být tahle chyba znatelná. Se zvětšující přesností se chyba postupně zmenšuje, tudíž aproximuje k nule. Výpočty byly prováděny na následující HW konfiguraci: CHIPSET : VIA KT-600 56-9

2006/56– 8.12.2006 CPU : AMD Athlon XP 2600++ Barthon RAM : 512MB RAM - 400Mhz DDR HDD : Hitachi 7k250 250Gb - accestime 8.5ms

5. Závěr V článku jsou popsány funkce a metody pro práci s maticemi komplexních čísel ve volitelné přesnosti. Základní aritmetika pro práci s čísly ve volitelné přesnosti byla použita z knihovny APFLOAT. Možnost využití knihovny je hlavně pro analýzu elektronických obvodů, i z toho důvodu že knihovna byla původně programována pro implementaci do programu Snap, který nepodporuje výpočty ve volitelné přesnosti. Z toho důvodu obsahuje knihovna taktéž podpůrné funkce, jako například funkce pro vytvoření admitační matice z netlistu vytvořeném právě v programu snap.

6. Literatura [1] Pfeifer, V.: Matematická knihovna pro práci s maticemi komplexních čísel s volitelnou přesností pro řešení elektronických obvodů, Diplomová práce, Ústav Radioelektroniky FEKT VUT Brno, 2006 [2] Tmomila, M.: APFLOAT - Arbitrary precision arithmetic package. Dostupné na WWW:http://www.apfloat.org

56-10