AK2-stuff/lab5-v3/Zadanie.md

Celem bieżącego ćwiczenia jest zapoznanie się z jednostką wektorową (MMX) połączone z praktycznym wykorzystaniem wiedzy z zakresu podstaw arytmetyki stałoprzecinkowej. Ćwiczenie należy zrealizować do dnia 4 czerwca 2020 r.

Opis jednostki MMX jest dostępny w dokumentacji procesora (Basic Architecture) rozdz. 9. Szczególną uwagę proszę zwrócić na równoległy model przetwarzania, tzn. na możliwość zapisania w jednym rejestrze wielu reprezentacji i wykonywania jednym rozkazem działań na wszystkich zapisanych reprezentacjach.

Opis zadania:

Napisać funkcję filter(...) obliczającą splot (https://en.wikipedia.org/wiki/Kernel_(image_processing)#Convolution) macierzy dwuwymiarowych.
Elementami macierzy będą reprezentacje w kodzie naturalnym binarnym.

Jedna z macierzy (oznaczona w poniższym opisie przez M) będzie dość duża (o rozmiarze co najmniej kilkaset na kilkaset elementów), druga (oznaczona przez k) będzie macierzą 3x3 elementy.
Wynikiem będzie duża macierz W o rozmiarze takim, jak macierz M.

W praktycznej realizacji elementy obu macierzy M i W będą pikselami obrazków, a elementy małej macierzy k będą specjalnie dobranymi stałymi pozwalającymi na uzyskanie konkretnego efektu.
Państwa implementacje powinny używać macierzy k postaci:

-1 -1  0
-1  0  1
 0  1  1

ponieważ dla takiej macierzy zostały przygotowane dane pozwalające na sprawdzenie poprawności implementacji.

Na operację splotu można patrzeć jak na:

• kolejne "nakładanie" małej macierzy k na kolejne pozycje w macierzy dużej M,
• mnożenie elementów dużej macierzy przez pokrywające się z nimi elementy macierzy małej,
• sumowanie uzyskanych iloczynów,
• zapis otrzymanej sumy jako elementu w macierzy wynikowej, który znajduje się na pozycji wyznaczonej przez położenie środkowego elementu macierzy k.

Przykład:

Macierz M, M(i,j) oznacza element na pozycji (i,j):

M(0,0) M(0,1) M(0,2) M(0,3) ...
M(1,0) M(1,1) M(1,2) M(1,3) ...
M(2,0) M(2,1) M(2,2) M(2,3) ...
M(3,0) M(3,1) M(3,2) M(3,3) ...
...    ...    ...    ...    ...

Macierz k:

k(0,0) k(0,1) k(0,2
k(1,0) k(1,1) k(1,2)
k(2,0) k(2,1) k(2,2)

---

Iteracja 1 (znakiem ~ oznaczono położenie środkowego elementu macierzy k):

"Nałożenie" macierzy:

k(0,0) k(0,1) k(0,2)  M(0,3) ...
k(1,0) ~k(1,1) k(1,2) M(1,3) ...
k(2,0) k(2,1) k(2,2)  M(2,3) ...
M(3,0) M(3,1) M(3,2)  M(3,3) ...
...    ...    ...     ...

Mnożenie i sumowanie:

W(1,1) =
M(0,0)*k(0,0) + M(0,1)*k(0,1) + M(0,2)*k(0,2) +
M(1,0)*k(1,0) + M(1,1)*k(1,1) + M(1,2)*k(1,2) +
M(2,0)*k(2,0) + M(2,1)*k(2,1) + M(2,2)*k(2,2)

---

Iteracja 2:

"Nałożenie" macierzy:

M(0,0) k(0,0) k(0,1)  k(0,2) ...
M(1,0) k(1,0) ~k(1,1) k(1,2) ...
M(2,0) k(2,0) k(2,1)  k(2,2) ...
M(3,0) M(3,1) M(3,2)  M(3,3) ...
...    ...    ...     ...

Mnożenie i sumowanie:

W(1,2) =
M(0,1)*k(0,0) + M(0,2)*k(0,1) + M(0,3)*k(0,2) +
M(1,1)*k(1,0) + M(1,2)*k(1,1) + M(1,3)*k(1,2) +
M(2,1)*k(2,0) + M(2,2)*k(2,1) + M(2,3)*k(2,2)

---

Jedna z kolejnych iteracji:

"Nałożenie" macierzy:

M(0,0) M(0,1) M(0,2)  M(0,3) ...
M(1,0) k(0,0) k(0,1)  k(0,2) ...
M(2,0) k(1,0) ~k(1,1) k(1,2) ...
M(3,0) k(2,0) k(2,1)  k(2,2) ...
...    ...    ...     ...

Mnożenie i sumowanie:

W(2,2) =
M(1,1)*k(0,0) + M(1,2)*k(0,1) + M(1,3)*k(0,2) +
M(2,1)*k(1,0) + M(2,2)*k(1,1) + M(2,3)*k(1,2) +
M(3,1)*k(2,0) + M(3,2)*k(2,1) + M(3,3)*k(2,2)

---

Państwa zadaniem będzie zaimplementowanie funkcji realizującej tę operację z wykorzystaniem jednostki MMX dostępnej w procesorach rodziny x86.
Zaimplementowana funkcja koniecznie musi w jak największym stopniu wykorzystywać możliwości wykonywania wielu operacji jedną instrukcją.
Oznacza to, że implementacja powinna jak najwięcej działań przeprowadzać równolegle. Implementacje wykorzystujące rozkazy MMX, ale przeprowadzające obliczenia bez równoległości ("pojedynczo", szeregowo) będą niżej punktowane.

Po zaimplementowaniu funkcji należy zmierzyć jej czas trwania korzystając z instrukcji rdtsc.
Po zmierzeniu czasu trwania należy także obliczyć i podać jakiś parametr pozwalający określić wydajność funkcji niezależnie od rozmiaru macierzy.
Przykładową wartością może być np. liczba cykli zegara na jeden element macierzy wynikowej.
Bardzo szybkie implementacje mogą uzyskać wyniki na poziomie pojedynczych cykli zegara na element, lub nawet mniej w niektórych przypadkach.
Podczas pomiarów wydajności proszę pamiętać o włączeniu optymalizacji kodu.

Uwagi techniczne

Pod adresem http://zak.iiar.pwr.wroc.pl/materials/architektura/laboratorium%20AK2/MMX/png.tar.bz2 dostępny jest kod programu (plik png.c), który:

• tworzy w pamięci reprezentację macierzy M na podstawie obrazka w formacie png o nazwie podanej jako argument programu,
• wywołuje funkcję, której napisanie jest Państwa zadaniem,
• wynik zapisuje w postaci obrazka w pliku out.png.

Oprócz kodu programu, w pliku png.tar.bz2 znajdują się także przykładowe obrazy będące danymi wejściowymi oraz spodziewane odpowiedzi (w plikach o nazwach zawierających _emboss).
Do sprawdzania uzyskanych wyników można np. wykorzystać jedną z technik wymienionych pod adresem https://askubuntu.com/questions/209517/does-diff-exist-for-images, są dostępne także inne narzędzia do porównywania obrazków.
Można także przygotować swoje funkcje porównujące wartości w poszczególnych pikselach, ale jest to rozwiązanie niezalecane ze względu na czasochłonność.
Proszę także wziąć pod uwagę, że w tym zadaniu możliwe są różne sposoby implementacji dające nieco inne wyniki.
Jeśli więc zaobserwowane różnice są niewielkie (na poziomie kilku najmniej znaczących bitów), to wyniki można uznać za poprawne.

Program png.c korzysta z biblioteki libpng.
Na moim systemie użyłem wersji 1.2, ale nie powinno to mieć znaczenia.
Biblioteka libpng powinna być automatycznie zainstalowana w większości systemów linuksowych.
W razie potrzeby powinno się ją dać bardzo łatwo doinstalować metodami standardowymi dla danego systemu.
Czasami może zachodzić potrzeba doinstalowania plików nagłówkowych, które zazwyczaj są w pakiecie zawierającym -dev w nazwie, np. libpng12-dev.

Kod programu png.c powstał jako dość swobodna kopia przykładów z dokumentacji biblioteki libpng: http://www.libpng.org/pub/png/pngbook.html
Działa dla załączonych obrazków, ale jeśli chcielibyście Państwo użyć także swoich danych, najprawdopodobniej należałoby albo przebudować program, albo odpowiednio przygotować obrazki wejściowe.
Szczegółów nie podaję, ponieważ jest ich dość dużo, a nie są bezpośrednio związane z treścią ćwiczenia.
Proszę także tego kodu nie traktować jako wzorcowego przykładu użycia biblioteki libpng, ponieważ ten kod powstał wyłącznie jako tymczasowe narzędzie ułatwiające realizację ćwiczenia.

W pierwszym kroku realizacji zadania należy skonstruować odpowiednie polecenia kompilacji i linkowania tak, aby uzyskać plik wykonywalny. Zazwyczaj wystarcza dodanie opcji -lpng podczas linkowania, ale zdecydowanie zalecam sprawdzenie w dokumentacji działania opcji gcc -I, -L i -l. W moim systemie pliki biblioteki są w katalogu /usr/lib/x86_64-linux-gnu/, a plik nagłówkowy w katalogu /usr/include.

Po uruchomieniu kodu zalecam najpierw napisanie dla rozgrzewki jakiejś prostej funkcji w C, która pozwoli Państwu nauczyć się swobodnie poruszać po reprezentacji macierzy M. Często stosowanym rozwiązaniem jest po prostu użycie dwóch pętli, ale są oczywiście możliwe także inne rozwiązania.
Dobrym pomysłem jest np. próba narysowania kreski lub paska na obrazku wynikowym. Najniższa wartość elementu odpowiada czarnemu pikselowi, najwyższa wartość elementu pikselowi białemu.

W pamięci operacyjnej macierz M zapisana jest "wierszami" (https://en.wikipedia.org/wiki/Row-_and_column-major_order), czyli pierwszym elementem (pod najniższym adresem) jest element M(0,0), potem M(0,1), M(0,2) ... M(1,0), M(1,1), ...

Po opanowaniu poruszania się po macierzy zdecydowanie zalecam najpierw napisanie w C prostej, szeregowej wersji funkcji realizującej zadany problem.
Funkcja ta powinna wszystkie obliczenia wykonywać korzystając wyłącznie z arytmetyki stałoprzecinkowej, ponieważ tylko taka jest dostępna w jednostce MMX.
Szczególną uwagę proszę zwrócić na dobranie odpowiednich typów danych dla wyników pośrednich.
Należy określić, jakie będą zakresy tych wyników i zagwarantować, że podczas obliczeń nie pojawią się błędy spowodowane przepełnieniem.
Błędy takie zazwyczaj są widoczne na obrazku wynikowym jako nagłe zmiany jasności sąsiednich pikseli.

Proszę także pamiętać, że elementy macierzy wynikowej W są zapisane w dokładnie takim samym formacie, jak elementy macierzy M.
Ponieważ więc wyniki obliczeń mogą znacznie przekraczać dopuszczalny zakres wartości, należy je LINIOWO przeskalować tak, aby:

• najmniejsza wartość elementu macierzy W odpowiadała najmniejszej wartości wyniku końcowego,
• największa wartość elementu macierzy W odpowiadała największej wartości wyniku końcowego,
• a pozostałe wartości wyniku odpowiadały wartościom elementów W z zachowaniem relacji większości, tzn. jeśli a_i < a_j, to W(a_i) <= W(a_j), gdzie a_i oznacza wartość wyniku, a W(a_i) wartość macierzy W wytworzoną z wyniku a_i.

Przykład:

Jeśli np. dopuszczalny zakres wartości elementów W to <0, 100>, a dopuszczalny zakres wartości wyniku a to <0, 1000>, to liniowe skalowanie można zrealizować zwykłym dzieleniem wyniku przez 10, czyli W(a) = a / 10.

---

Skalowanie można przeprowadzać w co najmniej dwóch miejscach:

• PRZED wykonaniem działań, czyli skalować od razu argumenty,
• PO wykonaniu działań, czyli skalować wynik końcowy.

Oczywiście skalowanie można także wykonywać w dowolnym momencie obliczeń.

Skalowanie PRZED wykonaniem obliczeń pozwala na używanie mniejszych typów danych, co może np. pozwolić na szybsze obliczenia w jednostce MMX, ponieważ więcej działań będzie można wykonać jedną instrukcją.
Skalowanie PRZED obliczeniami powoduje jednak większe błędy.

Skalowanie PO obliczeniach powoduje mniejsze błędy, ale wymaga przeprowadzania obliczeń na szerszych typach danych.

Z punktu widzenia celów tego ćwiczenia nie ma znaczenia, który z tych rodzajów skalowania Państwo wybierzecie, ważne jest jedynie, aby uzyskać poprawne wyniki.

Sztuczka techniczna – ponieważ w arytmetyce stałoprzecinkowej może być trudno zapisać dokładnie wymagane współczynniki skalowania, bardzo często stosuje się ich "rozsądne" przybliżenia.
W skrajnych przypadkach można nawet zastosować przybliżenie najbliższą potęgą 2, co sprowadza skalowanie do zwykłych przesunięć bitowych.
W Państwa implementacjach proszę spróbować poszukać rozsądnego kompromisu pomiędzy szybkością a dokładnością, ze zdecydowanym naciskiem na szybkość.
Załączone dane przykładowe powstały dla bardzo niedokładnych przybliżeń.

W operacji splotu pojawia się problem związany z właściwą obsługą wartości ze skrajnych wieszy/kolumn macierzy.
W tym ćwiczeniu poprawna obsługa tych wartości jest zupełnie nieistotna, dlatego zdecydowanie zalecam po prostu zignorowanie tych skrajnych kolumn i wierszy.
W załączonych danych skrajne wartości są ustawione przypadkowo, głównie na najmniejszą wartość.

Po opracowaniu poprawnej wersji w C, można przejść do realizacji właściwej części zadania, czyli do implementacji wykorzystującej równoległość udostępnianą przez jednostkę MMX.
Przy pierwszym kontakcie z modelem SIMD zalecam najpierw przemyślenie algorytmu i przygotowanie wstępnej koncepcji, np. w postaci obrazków pokazujących, jak będą ułożone dane w pamięci, w rejestrach, i jak na tych danych będą wykonywane równoległe operacje. Zdecydowanie pomaga rozrysowanie sobie, jak w pamięci zapisane są elementy macierzy M i W.

Następnie, przygotowanie kodu realizującego to zadanie najprościej będzie chyba Państwu zrobić po prostu w odzielnym pliku asemblerowym.
Możliwe są także różne techniki udostępniane przez narzędzie gcc, ale przy pierwszym kontakcie zdecydowanie odradzam ich stosowanie.
Jeśli jednak ktoś z Państwa ma doświadczenie w tym zakresie i potrafi poprawnie zrealizować to zadanie z poziomu kodu w C, to takie rozwiązania także są dopuszczalne. Jeszcze raz jednak przestrzegam osoby, które nigdy dotychczas nie spotkały się z takimi technikami, ponieważ wymagają one dogłębnego zrozumienia tematu.

Podczas przygotowywania implementacji MMX proszę skupić się wyłącznie na próbach opracowania kodu jak najlepiej wykorzystującego możliwości sprzętu.
Zalecam więc np.:

• wykorzystanie ustalonej postaci macierzy k (zakładamy, że ta macierz się na pewno nie zmieni, dlatego można znacznie uprościć wykonywane działania),
• całkowicie zignorować te skrajne wiersze i kolumny macierzy M, które nie będą "pasowały" Państwu do koncepcji,
• skalowanie wykonać w najprostszy możliwy sposób, nawet kosztem dokładności.