R2026.05.1 - filtrowanie obrazu przy pomocy SSE

README.md

@@ -1,12 +1,14 @@
 # AK2 - lab 2019+
-Trochę różnych programów napisanych przy pomocy z rozwiązywaniem zadań od mgr Aleksandry Postawki, mgra Tomasza Serafina, prof. Tadeusza Tomczaka (i innych) na zajęciach laboratoryjnych z Architektury komputerów 2, w formie materiałów dydaktycznych – kod szczegółowo opisany w komentarzach.
+Trochę różnych programów napisanych przy rozwiązywaniu zadań od mgr Aleksandry Postawki, mgra Tomasza Serafina, prof. Tadeusza Tomczaka (i innych) na zajęciach laboratoryjnych z Architektury komputerów 2, Wprowadzenia do wysokowydajnych komputerów oraz Organizacji i architektury komputerów, w formie materiałów dydaktycznych – kod szczegółowo opisany w komentarzach.
 Repozytorium powstało w 2019 i będzie na bieżąco aktualizowane.
 
 Przydatne linki:  
 * Intel 64 and IA-32 Architectures Software Developer’s Manual (Combined Volumes) – opis architektury procesorów x86 i spis wszystkich rozkazów asemblera:  
-https://software.intel.com/sites/default/files/managed/39/c5/325462-sdm-vol-1-2abcd-3abcd.pdf
+https://www.intel.com/content/www/us/en/content-details/671200/intel-64-and-ia-32-architectures-software-developer-s-manual-combined-volumes-1-2a-2b-2c-2d-3a-3b-3c-3d-and-4.html
 * Linux System Call Table – spis funkcji systemowych Linuxa dla 32- i 64-bitowej wersji x86 w przejrzystej formie:  
-https://chromium.googlesource.com/chromiumos/docs/+/HEAD/constants/syscalls.md
+https://www.chromium.org/chromium-os/developer-library/reference/linux-constants/syscalls/
+* Programming from the Ground Up:  
+https://download-mirror.savannah.gnu.org/releases/pgubook/ProgrammingGroundUp-1-0-booksize.pdf
 * System V Application Binary Interface AMD64 – opis 64-bitowej konwencji wywołań x86 używanej w systemach z rodziny Linux:  
 https://www.uclibc.org/docs/psABI-x86_64.pdf
 * System V Application Binary Interface i386 – opis 32-bitowej konwencji wywołań x86 używanej w systemach z rodziny Linux:  
@@ -23,3 +25,5 @@ https://en.cppreference.com/w/
     https://en.cppreference.com/w/cpp/io/c/fprintf
     * Opis funkcji scanf i jej ciągów formatujących:  
     https://en.cppreference.com/w/cpp/io/c/fscanf
+* Using the GNU Compiler Collection (GCC) – Extended Asm – składnia wstawek asemblerowych w C:  
+https://gcc.gnu.org/onlinedocs/gcc/Extended-Asm.html

fibb_tsc/timestamp32.s

@@ -6,8 +6,11 @@
 .text
 
 timestamp:
+push %ebx
+
 xor %eax, %eax
 cpuid
 rdtsc
 
+pop %ebx
 ret

fibb_tsc/timestamp64.s

@@ -6,6 +6,8 @@
 .text
 
 timestamp:
+push %rbx
+
 xor %rax, %rax
 cpuid
 rdtsc
@@ -13,4 +15,5 @@ rdtsc
 shl $32, %rdx
 or %rdx, %rax
 
+pop %rbx
 ret

lab5-v3/Makefile

@@ -0,0 +1,7 @@
+all: png
+
+png: png.c filter.s timestamp64.s
+	gcc -g -lpng -O3 -no-pie png.c filter.s timestamp64.s -o png
+
+clean:
+	rm png

lab5-v3/README.md

@@ -0,0 +1,218 @@
+### UWAGA 2026
+Mój kod napisany w 2020 jest wariacją na temat zadania prof. Tomczaka, którego pełny opis zamieściłem poniżej, a nie jego dokładną implementacją. Stąd wykorzystanie rozkazów SSE, a nie MMX.  
+Ma pozwolić zrozumieć ideę posługiwania się operacjami SIMD, nie polecam jednak go bezkrytycznie kopiować.  
+
+Przy testach na komputerze z procesorem AMD Ryzen 5 2600 X i 16 GB pamięci DDR4 3000 MHz osiągnąłem tempo około 2.5 cykli zegara na piksel.
+
+---
+
+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.

lab5-v3/filter.s

@@ -0,0 +1,260 @@
+# Jan Potocki 2020
+
+.globl filter
+.type filter, @function
+
+.data
+# Pozwala uzyc rozkazu movdqa zamiast movdqu
+.align 16
+
+# Macierz k w 4 kopiach - do rejestru xmm zmiescimy dane dla 4 pikseli splotu 
+# 1. wiersz macierzy dopelniony do 32 bitow
+kernel128_1:
+.byte -1, -1, 0, 0
+.byte -1, -1, 0, 0
+.byte -1, -1, 0, 0
+.byte -1, -1, 0, 0
+# 2. wiersz macierzy dopelniony do 32 bitow
+kernel128_2:
+.byte -1, 0, 1, 0
+.byte -1, 0, 1, 0
+.byte -1, 0, 1, 0
+.byte -1, 0, 1, 0
+# 3. wiersz macierzy dopelniony do 32 bitow
+kernel128_3:
+.byte 0, 1, 1, 0
+.byte 0, 1, 1, 0
+.byte 0, 1, 1, 0
+.byte 0, 1, 1, 0
+
+# Stala na potrzeby normalizacji - wartosc splotu z k miesci sie w przedziale 
+# od 3*(-1*255) do 3*(1*255) czyli <-765,765>, wiec dodanie 765 pozwoli
+# przesunac ja w przedzial <0,1530>.
+# Na tym etapie obliczen liczby beda zapisane w 16-bitowej reprezentacji.
+bias128: .word 765, 765, 765, 765, 765, 765, 765, 765
+
+# Argumenty:
+# rdi - adres M
+# rsi - adres W
+# rdx - rozdzielczosc w poziomie
+# rcx - rozdzielczosc w pionie
+
+.text
+
+filter:
+push %rbp
+mov %rsp, %rbp
+push %r15
+push %r14
+push %r13
+push %r12
+
+# r15 - warunek stopu dla petli SIMD, tyle ile zmiesci sie w wierszu
+mov %rdx, %r15
+sub $17, %r15
+# r14 - warunek stopu dla petli indeksujacej wiersze
+mov %rcx, %r14
+sub $2, %r14
+# r12 - warunek stopu dla petli pomocniczej, dokonczenie wiersza
+mov %rdx, %r12
+sub $2, %r12
+# ...w tej implementacji skrajne piksele obrazu wynikowego pomijamy
+
+# Wpisanie stalych do rejestrow aby przyspieszyc obliczenia
+movdqa bias128, %xmm12
+movdqa kernel128_1, %xmm13
+movdqa kernel128_2, %xmm14
+movdqa kernel128_3, %xmm15
+
+# r9 - pomocniczy indeks adresujacy wiersze w pamieci
+xor %r9, %r9
+# r10 - indeks rozdzielczosci pionowej (wierszy)
+xor %r10, %r10
+
+# Petla indeksujaca wiersze
+lp1:
+# r11 - indeks rozdzielczosci poziomej (pikseli w wierszu)
+xor %r11, %r11
+
+# Podstawowa petla SIMD (16 pikseli na iteracje)
+lp2:
+# r8 - adres pikseli w pamieci (suma adresu wiersza i indeksu pikseli)
+mov %r9, %r8
+add %r11, %r8
+
+# Instrukcja pinsrd kopiuje 4 kolejne 8-bitowe piksele do fragmentu xmm,
+# czyli tyle ile trzeba aby obliczyc splot z jednym wierszem k. Ostatni piksel
+# zachodzi na pierwszy kolejnego splotu - to nie ma znaczenia, bo kazdy wiersz
+# macierzy k jest do 32-bitowej postaci dopelniony zerem na ostatnim bajcie.
+# Instrukcja pmaddubsw oblicza iloczyn 8-bitowych pikseli obrazu i wartosci 
+# kernela, sumuje dwa sasiednie iloczyny i w ich miejsce zapisuje wynik
+# w 16-bitowej reprezentacji U2.
+
+# Gorny wiersz macierzy k
+pinsrd $0, 0(%rdi, %r8, 1), %xmm0
+pinsrd $1, 1(%rdi, %r8, 1), %xmm0
+pinsrd $2, 2(%rdi, %r8, 1), %xmm0
+pinsrd $3, 3(%rdi, %r8, 1), %xmm0
+pmaddubsw %xmm13, %xmm0
+
+pinsrd $0, 4(%rdi, %r8, 1), %xmm1
+pinsrd $1, 5(%rdi, %r8, 1), %xmm1
+pinsrd $2, 6(%rdi, %r8, 1), %xmm1
+pinsrd $3, 7(%rdi, %r8, 1), %xmm1
+pmaddubsw %xmm13, %xmm1
+
+pinsrd $0, 8(%rdi, %r8, 1), %xmm6
+pinsrd $1, 9(%rdi, %r8, 1), %xmm6
+pinsrd $2, 10(%rdi, %r8, 1), %xmm6
+pinsrd $3, 11(%rdi, %r8, 1), %xmm6
+pmaddubsw %xmm13, %xmm6
+
+pinsrd $0, 12(%rdi, %r8, 1), %xmm7
+pinsrd $1, 13(%rdi, %r8, 1), %xmm7
+pinsrd $2, 14(%rdi, %r8, 1), %xmm7
+pinsrd $3, 15(%rdi, %r8, 1), %xmm7
+pmaddubsw %xmm13, %xmm7
+
+add %rdx, %r8
+# Przechowanie adresu srodkowego wiersza na pozniej, do zapisu wyniku splotu
+mov %r8, %r13
+
+# Srodkowy wiersz macierzy k
+pinsrd $0, 0(%rdi, %r8, 1), %xmm2
+pinsrd $1, 1(%rdi, %r8, 1), %xmm2
+pinsrd $2, 2(%rdi, %r8, 1), %xmm2
+pinsrd $3, 3(%rdi, %r8, 1), %xmm2
+pmaddubsw %xmm14, %xmm2
+
+pinsrd $0, 4(%rdi, %r8, 1), %xmm3
+pinsrd $1, 5(%rdi, %r8, 1), %xmm3
+pinsrd $2, 6(%rdi, %r8, 1), %xmm3
+pinsrd $3, 7(%rdi, %r8, 1), %xmm3
+pmaddubsw %xmm14, %xmm3
+
+pinsrd $0, 8(%rdi, %r8, 1), %xmm8
+pinsrd $1, 9(%rdi, %r8, 1), %xmm8
+pinsrd $2, 10(%rdi, %r8, 1), %xmm8
+pinsrd $3, 11(%rdi, %r8, 1), %xmm8
+pmaddubsw %xmm14, %xmm8
+
+pinsrd $0, 12(%rdi, %r8, 1), %xmm9
+pinsrd $1, 13(%rdi, %r8, 1), %xmm9
+pinsrd $2, 14(%rdi, %r8, 1), %xmm9
+pinsrd $3, 15(%rdi, %r8, 1), %xmm9
+pmaddubsw %xmm14, %xmm9
+
+add %rdx, %r8
+
+# Dolny wiersz macierzy k
+pinsrd $0, 0(%rdi, %r8, 1), %xmm4
+pinsrd $1, 1(%rdi, %r8, 1), %xmm4
+pinsrd $2, 2(%rdi, %r8, 1), %xmm4
+pinsrd $3, 3(%rdi, %r8, 1), %xmm4
+pmaddubsw %xmm15, %xmm4
+
+pinsrd $0, 4(%rdi, %r8, 1), %xmm5
+pinsrd $1, 5(%rdi, %r8, 1), %xmm5
+pinsrd $2, 6(%rdi, %r8, 1), %xmm5
+pinsrd $3, 7(%rdi, %r8, 1), %xmm5
+pmaddubsw %xmm15, %xmm5
+
+pinsrd $0, 8(%rdi, %r8, 1), %xmm10
+pinsrd $1, 9(%rdi, %r8, 1), %xmm10
+pinsrd $2, 10(%rdi, %r8, 1), %xmm10
+pinsrd $3, 11(%rdi, %r8, 1), %xmm10
+pmaddubsw %xmm15, %xmm10
+
+pinsrd $0, 12(%rdi, %r8, 1), %xmm11
+pinsrd $1, 13(%rdi, %r8, 1), %xmm11
+pinsrd $2, 14(%rdi, %r8, 1), %xmm11
+pinsrd $3, 15(%rdi, %r8, 1), %xmm11
+pmaddubsw %xmm15, %xmm11
+
+# Suma kolumn splotu
+paddw %xmm2, %xmm0
+paddw %xmm3, %xmm1
+paddw %xmm8, %xmm6
+paddw %xmm9, %xmm7
+paddw %xmm4, %xmm0
+paddw %xmm5, %xmm1
+paddw %xmm10, %xmm6
+paddw %xmm11, %xmm7
+
+# Suma wierszy splotu
+phaddw %xmm1, %xmm0
+phaddw %xmm7, %xmm6
+
+# Normalizacja wyniku splotu - dodanie zdefiniowanego obciazenia i skalowanie
+# z powrotem do 8-bitowej reprezentacji.
+# Dokladny wynik pozwoliloby uzyskac dzielenie przez 6 (1530/6=255),
+# przesuniecie o 3 bity w prawo jest wydajniejsze, a wciaz pozwala uzyskac
+# rozsadna dokladnosc (1530>>3=191)
+paddw %xmm12, %xmm0
+paddw %xmm12, %xmm6
+psrlw $3, %xmm0
+psrlw $3, %xmm6
+
+# Spakowanie 8-bitowych wynikow do jednego rejestru i zapis do pamieci
+packuswb %xmm6, %xmm0
+movdqu %xmm0, 1(%rsi, %r13, 1)
+
+# Sprawdzenie, czy nastepna iteracja petli SIMD zmiesci sie w wierszu
+add $16, %r11
+cmp %r15, %r11
+jl lp2
+
+# Sprawdzenie, czy wiersz zostal przeliczony w calosci...
+cmp %r12, %r11
+je next_row
+
+# ...jesli nie - na dobieg pomocnicza petla (1 piksel na iteracje)
+lp2a:
+mov %r9, %r8
+add %r11, %r8
+
+# Gorny wiersz macierzy k
+pinsrd $0, 0(%rdi, %r8, 1), %xmm0
+pmaddubsw %xmm13, %xmm0
+add %rdx, %r8
+mov %r8, %r13
+
+# Srodkowy wiersz macierzy k
+pinsrd $0, 0(%rdi, %r8, 1), %xmm2
+pmaddubsw %xmm14, %xmm2
+add %rdx, %r8
+
+# Dolny wiersz macierzy k
+pinsrd $0, 0(%rdi, %r8, 1), %xmm4
+pmaddubsw %xmm15, %xmm4
+
+# Suma i normalizacja splotu
+paddw %xmm2, %xmm0
+paddw %xmm4, %xmm0
+phaddw %xmm0, %xmm0
+paddw %xmm12, %xmm0
+psrlw $3, %xmm0
+
+# Spakowanie 8-bitowego wyniku w najmlodszym fragmencie rejestru i zapis
+packuswb %xmm0, %xmm0
+pextrb $0, %xmm0, 1(%rsi, %r13, 1)
+
+# Sprawdzenie, czy tym razem wiersz jest juz przeliczony w calosci
+inc %r11
+cmp %r12, %r11
+jl lp2a
+
+next_row:
+# Przesuniecie adresu wiersza w dol
+add %rdx, %r9
+# Sprawdzenie, czy caly obraz jest przeliczony
+inc %r10
+cmp %r14, %r10
+jl lp1
+
+pop %r12
+pop %r13
+pop %r14
+pop %r15
+pop %rbp
+ret

lab5-v3/png.c

@@ -0,0 +1,147 @@
+#include <png.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+
+#define ERROR                                                   \
+	fprintf (stderr, "ERROR at %s:%d.\n", __FILE__, __LINE__) ;   \
+	return -1 ;
+
+#define HEADER_SIZE (1)
+
+void filter(unsigned char * M, unsigned char * W, int width, int height);
+unsigned long long timestamp();
+
+int main (int argc, char ** argv)
+{
+	if (2 != argc)
+	{
+		printf ("\nUsage:\n\n%s file_name.png\n\n", argv[0]) ;
+
+		return 0 ;
+	}
+
+	const char * file_name = argv [1] ;
+	
+	unsigned char header [HEADER_SIZE] ;
+
+	FILE *fp = fopen (file_name, "rb");
+	if (NULL == fp)
+	{
+		fprintf (stderr, "Can not open file \"%s\".\n", file_name) ;
+		ERROR
+	}
+
+	if (fread (header, 1, HEADER_SIZE, fp) != HEADER_SIZE)
+	{
+		ERROR
+	}
+
+	if (0 != png_sig_cmp (header, 0, HEADER_SIZE))
+	{
+		ERROR
+	}
+
+	png_structp png_ptr = 
+		png_create_read_struct
+			(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL );
+	if (NULL == png_ptr)
+	{
+		ERROR
+	}
+
+	png_infop info_ptr = png_create_info_struct (png_ptr);
+	if (NULL == info_ptr)
+	{
+		png_destroy_read_struct (& png_ptr, NULL, NULL);
+
+		ERROR
+	}
+
+	if (setjmp (png_jmpbuf (png_ptr))) 
+	{
+		png_destroy_read_struct (& png_ptr, & info_ptr, NULL);
+
+		ERROR
+	}
+
+	png_init_io       (png_ptr, fp);
+	png_set_sig_bytes (png_ptr, HEADER_SIZE);
+	png_read_info     (png_ptr, info_ptr);
+
+	png_uint_32  width, height;
+	int  bit_depth, color_type;
+	
+	png_get_IHDR
+	(
+		png_ptr, info_ptr, 
+		& width, & height, & bit_depth, & color_type,
+		NULL, NULL, NULL
+	);
+
+	if (8 != bit_depth)
+	{
+		ERROR
+	}
+	if (0 != color_type)
+	{
+		ERROR
+	}
+
+	size_t size = width ;
+	size *= height ;
+
+	unsigned char * M = malloc (size) ;
+
+	png_bytep ps [height] ;
+	ps [0] = M ;
+	for (unsigned i = 1 ; i < height ; i++)
+	{
+		ps [i] = ps [i-1] + width ;
+	}
+	png_set_rows (png_ptr, info_ptr, ps);
+	png_read_image (png_ptr, ps) ;
+
+	printf 
+	(
+		"Image %s loaded:\n"
+		"\twidth      = %lu\n"
+		"\theight     = %lu\n"
+		"\tbit_depth  = %u\n"
+		"\tcolor_type = %u\n"
+		, file_name, width, height, bit_depth, color_type
+	) ;
+
+	unsigned char * W = malloc (size) ;
+
+    unsigned long long tstart = timestamp();
+ 	filter (M, W, width, height) ;
+    unsigned long long tstop = timestamp();
+    
+    printf("\nTSC cycles: %llu\n", tstop-tstart);
+    
+	png_structp write_png_ptr =
+		png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
+	if (NULL == write_png_ptr)
+	{
+		ERROR
+	}
+
+	for (unsigned i = 0 ; i < height ; i++)
+	{
+		ps [i] += W - M ;
+	}
+	png_set_rows (write_png_ptr, info_ptr, ps);
+
+	FILE *fwp = fopen ("out.png", "wb");
+	if (NULL == fwp)
+	{
+		ERROR
+	}
+
+	png_init_io   (write_png_ptr, fwp);
+	png_write_png (write_png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);
+	fclose (fwp);
+
+	return 0;
+}

lab5-v3/timestamp64.s

@@ -0,0 +1,19 @@
+# Jan Potocki 2020
+
+.global timestamp
+
+# Segment kodu
+.text
+
+timestamp:
+push %rbx
+
+xor %rax, %rax
+cpuid
+rdtsc
+
+shl $32, %rdx
+or %rdx, %rax
+
+pop %rbx
+ret