Biuletyn nr 22

Biuletyn KDM
1 | 2 | 3 | 4 | 5
6 | 7 | 8 | 9 | 10
11 | 12 | 13 | 14
15 | 16 | 17 | 18
19 | 20 | 21 | 22
23 | 24 | 25 | 26
27 | 28 | 29 | 30
31 | 32
Lista biuletynów

Biuletyn nr 22 (21 czerwca 2007)

Spis treści

Narzędzia: Vi - konfiguracja i przydatne opcje

Autor: Maciej Szpindler

Pracując w trybie zdalnym, przykładowo przy użyciu SSH, często niezbędnym narzędziem jest edytor tekstu. Nie zawsze w tej sytuacji wygodne jest korzystanie z graficznych edytorów. Z pomocą przychodzi jeden z najpopularniejszych i obecnych na wiekszości platform edytor - vi.

Wielu osobom włos jeży się na glowie na dzwięk nazwy vi. Przecież to program z epoki kamienia łupanego!?

Faktycznie vi to bardzo stary edytor, ale posiada ogromne możliwości. Z początku może wydawać się trudny, ale po zapoznaniu się z podstawami jego działania korzystanie z edytora staje się całkiem wygodne.

Jak działa vi

Vi to edytor modalny. Modalność oznacza, że podczas pracy możemy przełączać się miedzy różnymi trybami; w uproszczeniu między trybem wstawiania i trybem poleceń.

Uruchamiając edytor poleceniem vi lub vi nazwa_pliku znajdziemy się domyślnie w trybie poleceń.

W edytorze można poruszać się za pomocą kursorów.

Aby wejść w tryb edycji należy nacisnąć klawisz i insert lub a append.

W tryb poleceń przełączamy się naciskając klasisz ESC. W trybie poleceń może wydawać komendy poprzedzając je dwukropkiem : i zatwierdzając klawiszem Enter.

Vim czyli VI iMproved

Nieco młodszym klonem edytora vi jest tzw. vim, całkowicie kombatybilny z pierwowzorem, wzbogacony o nową funkcjonalność. W tym krótkim wprowadzeniu zakładamy, że używamy vim-a. W wielu nowych dystrubucjach Linuxa vim jest domyślną wersją edytora vi.

Podstawowe komendy

Pomoc

Aby uzyskać pomoc w trybie poleceń należy napisać :help komenda lub :h komenda.

Nawigacja

W trybie poleceń (bez :):

^ $ - pierwszy / ostatni znak danej linii

gg G - pierwsza / ostatnia linia w pliku

n% - linia odpowiadająca n% pliku

:q! Enter - zamknij bez zapisania zmian

:wq Enter - zapisz i zamknij plik

Edycja

i a - wstawianie przed / za znakiem kursora

I A - wstawianie na początku / końcu linii

o O - nowa linia pod / nad bieżącą linią

X X - usunięcie znaku pod / przed znakiem kursora

dd D - usunięcie bieżącej linii / linii do końca pliku

:rd Enter - usunięcie r-tej linii

yy - skopiowanie bieżącej linii do rejestru

p PP - wstawienie zawartości rejestru pod / nad bieżącą linią

u - undo

Tryb Visual

v V - rozpocznij / zakończ podświetlanie znaków

Wyszukiwanie

/s ?s - szukaj w przód / w tył ciągu znaków s

n - powtórz wyszukiwanie w przód

N - powtórz wyszukiwanie w tył

# * - szukaj w przód / w tył ciągu znaków pod znakiem kursora

Inne

: new Enter - otwórz nowe okno edytora

: new plik Enter - otwórz plik w nowym oknie

:r f - wstaw zawartość pliku f za znakiem kursora

:r c - wstaw wynik komendy c za znakiem kursora

:r w - zapisz zakres r (np. numery linii oddzielone ,) do bieżącego pliku

:r w f - zapisz zakres r do pliku f

:sh Enter :!c Enter - przejście do shella / wykonanie komendy c w shellu

Konfiguracja vima

Vim-a możemy dostosować do swoich potrzeb modyfikując (lub tworząc) plik startowy "vimrc", powinien znajdować się o w katalogu domowym i nazywać .vimrc.

Plik zawiera komendy, które zostaną wywołane automatycznie przy starcie edytora. Każda linia zawiera komendę, którą można wykonać w trybie poleceń.

Przykładowe opcje

Wygodne może być umieszczenie w pliku strartowym vim-a przykładowych linii:

syntax on

Włącza podświetlanie składni w edytorze, obsługuje wiele języków programowania.

set nu

Powoduje wyświetlanie numeru linii na początku każdego wiersza.

set tabstop=n

Powoduje zamianę tabulatorów na ciąg n spacji.

set backupdir=jakis_katalog

Ustawienie zapisywania plików z "~" na końcu do katalogu jakis_katalog

set directory=jakis_katalog

Ustawienie zapisywania plików z ".swp" do katalogu jakis_katalog

Vim na tornado

W przypadku komputera Tornado vi i vim to dwa różne edytory. Tornado to dość nietypowy komputer, jego system operacyjny to UNICOS (a dokładniej Irix w wersji na tą architekturę) i ma on swoje rózne smaczki. W szczególności dla wielu osób uciążliwe mogą być nie działające domyślnie (w tradycyjny sposób) w vim-ie klawisze Delete i Backspace. Opiszemy krótko jak to można zmienić.

Zmieniamy zawartość pliku .vimrc w katalogu domowym dodając następujące linie:

:if &term == "xterm"
:  set t_kb=^?
:  fixdel
:endif

Onaczają one, że jeżeli nasz terminal działa w trybie xterm przypisujemy klawiszowi Backspace - w vi t_kb - ciąg znaków, które spowodują upragniony efekt po jego naciśnięciu.

Uwaga: ciąg znaków ^? należy wpisać za pomocą kombinacji klawiszy CTRL+v CTRL+Backspace.

Natomiast komenda fixdel poprawia działanie klawisza Delete.

Dalej wprowadzamy analogiczne linie dla trybu vt100:

:if &term == "vt100"
:  set t_kb=^?
:  fixdel
:endif

Potrzebna jest jeszcze jedna linia:

set bs=2

Powodująca, że działanie Backspace nie będzie blokowane.

Biblioteki: Obliczanie splotu funkcji przy wykorzystaniu biblioteki FFTW

Autor: Maciek Cytowski

Problem

Obliczanie splotu dwóch funkcji jest zagadnieniem, które pojawia się dość często w naukach obliczeniowych. Przykładami zastosowań, z którymi osobiście się zetknąłem są np. rozmycie obrazu cyfrowego lub wygładzanie nieciągłej funkcji w celu użycia jej do dalszych obliczeń np. w solwerze FEM (Finite Element Method). Obydwa te zagadnienia można rozwiązać poprzez wyliczenie splotu z funkcją Gaussa, która w 1D ma postać:

G_{\sigma} (x) = \frac{1}{2  \Pi \sigma} e^{\frac{-x^{2}}{2 \sigma ^{2}}}

Ogólny wzór na wyliczenie splotu dwóch funkcji f:R \rightarrow R oraz g:R \rightarrow R jest następujący:

(f*g)(t) = \int^{\infty}_{-\infty} f(\tau) g(t-\tau) d\tau

Bardzo klarowne wizualne przedstawienie konceptu splotu funkcji można znaleźć na stronach Wolfram Math World.

W przypadku rzeczywistych obliczeń numerycznych mamy oczywiście zwykle do czynienia z danymi dyskretnymi. Funkcja f jest wówczas reprezentowana przez dyskretne próbki f_{j} mierzone w równych przerwach czasowych. Funkcja g jest również dyskretnym zbiorem próbek g_{k}, o następującej interpretacji:

  • r_{0} mówi nam jaka wielokrotność j-tej wartości wektora f jest kopiowana do j-tego miejsca wektora wynikowego (i tak samo dla wszelkich pozycji j wektora f)
  • r_{1} mówi nam jaka wielokrotność j-tej wartości wektora f jest dodatkowo kopiowana do j+1-ego miejsca wektora wynikowego (i tak samo dla wszelkich pozycji j wektora f)
  • r_{-1} mówi nam jaka wielokrotność j-tej wartości wektora f jest dodatkowo kopiowana do j-1-ego miejsca wektora wynikowego (i tak samo dla wszelkich pozycji j wektora f)


Wówczas wzór na splot funkcji f z funkcją g niezerową w zakresie o długośći M ma następującą postać: (f*g)_{j}=\sum_{k=-M/2+1}^{M/2} f_{j-k}r_{k}


Łatwo zauważyć, że obliczenie j-ego elementu wektora wynikowego kosztuje nas M mnożeń i dodawań. Na szczęście z pomocą przychodzi nam twierdzenie o dyskretnym splocie:

Jeśli funkcja f_{j} jest okresowa z okresem N, tak że jest całkowicie zdefiniowana przez wartości f_{0},\ldots ,f_{N-1} wówczas dyskretna transformata Fouriera jej dyskretnego splotu z funkcją g zdefiniowaną wartościami g_{0},\ldots ,g_{N-1} równa jest iloczynowi dyskretnych transformat Fouriera funkcji f oraz g.

W powyższym twierdzeniu zakładamy, że funkcja g jest zapisana w porządku nazywanym w literaturze wrap-around. Wektor wartości g_{0},\ldots ,g_{N-1} zawiera te same wartości co oryginalny wektor g_{-N/2+1},\ldots ,g_{N/2}, ale wartości o indeksach ujemnych są przestawione i odbite na jego końcu. Sytuację tą przedstawiają rysunki zamieszczone poniżej.

Co więcej przy obliczeniach będziemy również dokonywali niewielkiego preprocessingu funkcji f. Aby mieć pewność, że wartości z końca wektora f nie będą miały wpływu na wartości z początku wektora wynikowego często dokonuje się tzw. dopełniania zerami (zero padding).


Algorytm wyliczenia splotu dwóch dyskretnych funkcji f i g będzie zatem następujący:

  • wylicz dyskretną transformatę Fouriera funkcji f
  • wylicz dyskretną transformatę Fouriera funkcji g
  • dokonaj mnożenia dyskretnych transformat Fouriera
  • wylicz odwrotną transformatę Fouriera wyniku

Rozwiązanie

Do rozwiązania tego zadania użyjemy bardzo popularnej biblioteki FFTW (The Fastest Fourier Transform on The West). W ICM zainstalowana jest ona na klastrze halo.

Bardzo dokładny opis składni funkcji bibliotecznych znajduje się na stronie dokumentacji.

Załóżmy, żę funkcja f jest funkcją "schodkową" i opisana jest na odcinku [0,1] następującym wzorem:


f(x) = 
\begin{cases} 
  0, & \mbox{if }x \in [0,1/8)\\
  10, & \mbox{if }x \in [1/8,1/4)\\
  -5, & \mbox{if }x \in [1/4,3/8)\\
  10, & \mbox{if }x \in [3/8,1/2)\\
  5, & \mbox{if }x \in [1/2,5/8)\\
  -10, & \mbox{if }x \in [5/8,3/4)\\
  5, & \mbox{if }x \in [3/4,7/8)\\
  0, & \mbox{if }x \in [7/8,1]
\end{cases}


Splot funkcji f z funkcją Gaussa możemy zaprogramować w języku C np. w podany poniżej sposób:


#include "stdio.h"
#include "math.h"
#include "fftw3.h"

int main()
{
        int N;                  // ilosc dyskretnych pomiarow probki
        int K;                  // dlugosc tzw. "zero padding"
        double sigma;           // parametr sigma funkcji Gaussa
        double exponent;        // zmienna pomocnicza
        double epsilon;         // parametr obcięcia "ogona" funkcji Gaussa
        double *result;         // tablica na wynik
        double *f,*G;           // funkcja f oraz funkcja Gaussa
        int i,j;
        double C;

        epsilon=0.000001f; // wartość poniżej której wartości funkcji G będą przypisane zera

        sigma=10;

        N=2000;

        C=(1/(sqrt(2*M_PI)*sigma));

        for(exponent=C,K=0;exponent>0.0f;K++) {
                exponent=C*exp(-K*K/(2*sigma*sigma));
        }

        /* Alokowanie pamieci na funkcje Gaussa */

        G=(double*) malloc((N+K)*sizeof(double));

        /* Wypełnienie tablicy G[] wartościami - zgodnie z regułą "wrap-around" */

        for(i=0;i<=(N+K)/2;i++) {
                exponent = C*exp(-i*i/(2*sigma*sigma));
                if(exponent>epsilon) G[i]=exponent;
                else G[i]=0;
        }
        for(i=(N+K)-1;i>(N+K)/2;i--) {
                exponent = C*exp(-(i-N-K)*(i-N-K)/(2*sigma*sigma));
                if(exponent>epsilon) G[i]=exponent;
                else G[i]=0;
        }

        /* Alokowanie pamieci na funkcje f */

        f = (double*) malloc((N+K)*sizeof(double));

        /* Wypełnienie tablicy f[] wartościami w miejscach od 0 do N-1. Reszta przeznaczona na tzw. "zero padding" */

        for(i=0;i<N;i++) {
                if(i<N/8) f[i]=0;
                if(i>=N/8 && i<N/4) f[i]=10;
                if(i>=N/4 && i<3*N/8) f[i]=-5;
                if(i>=3*N/8 && i<N/2) f[i]=10;
                if(i>=N/2 && i<5*N/8) f[i]=5;
                if(i>=5*N/8 && i<3*N/4) f[i]=-10;
                if(i>=3*N/4 && i<7*N/8) f[i]=5;
                if(i>=7*N/8 && i<N) f[i]=0;
        }

        /* Reszta przeznaczona na tzw. "zero padding" */

        for(i=N;i<N+K;i++) f[i]=0;

        result=(double*) malloc((N+K+1)*sizeof(double));

        /* Obliczenia FFT */
        {
                fftw_plan p,q,r;
                fftw_complex *fft_signal,*fft_response;

                fft_signal=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*(N+K+1));
                fft_response=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*(N+K+1));

                p=fftw_plan_dft_r2c_1d(N+K,f,fft_signal,FFTW_ESTIMATE);
                r=fftw_plan_dft_r2c_1d(N+K,G,fft_response,FFTW_ESTIMATE);

                fftw_execute(p);
                fftw_execute(r);

                for(i=0;i<N+K;i++){
                        double tmp=fft_signal[i][0];
                        double scale=1/(double)(N+K);
                        fft_signal[i][0]=(tmp*fft_response[i][0]-fft_signal[i][1]*fft_response[i][1])*scale;
                        fft_signal[i][1]=(tmp*fft_response[i][1]+fft_response[i][0]*fft_signal[i][1])*scale;
                }

                q=fftw_plan_dft_c2r_1d(N+K,fft_signal,result,FFTW_ESTIMATE);

                fftw_execute(q);

                free(fft_signal);
                free(fft_response);

        }

        fftw_free(result);
        fftw_free(G);
        fftw_free(f);

}


Załączone powyżej komentarze (kolor zielony) powinny w skrócie objaśnić działanie programu i jego kolejne kroki. Ograniczę się zatem tylko do opisania części programu wykorzystującej bibliotekę FFTW. Po pierwsze na początku programu załączamy plik nagłówkowy defininiujący składnię FFTW pisząc:

#include "fftw3.h

Następnie korzystając z dostępnych funkcji oraz typów programujemy algorytm dyskretnego splotu przedstawiony w poprzedniej sekcji (kolor czerwony).

Na samym początku definiujemy i alokujemy pamięć dla tablic fftw_signal oraz fftw_response (o typie zespolonym fftw_complex, który jest tożsamy typowi double[2]). Możemy do tego celu skorzystać z funkcji alokacji FFTW: fftw_malloc. Równie dobrze moglibyśmy skorzystać z klasycznej funkcji malloc, ale jej odpowiednik z biblioteki FFTW potrafi optymalnie alokować pamięć jeśli nasza architektura udostępnia instrukcje wektorowe SIMD (takie jak SSE czy Altivec). Następnie przygotowujemy tzw. plan (możemy go na użytek tego artykułu nazywać planem działań). Aby zdefiniować plan działań musimy wiedzieć jakiego typu transformaty FFTW chcemy dokonać, jakie są dane wejściowe, jaka jest długość transformaty oraz gdzie chcemy zapisać wynik tego przekształcenia. Przykładowo linia:

p=fftw_plan_dft_r2c_1d(N+K,f,fft_signal,FFTW_ESTIMATE);

definiuje plan działania dla obliczenia jednowymiarowej transformaty FFT z danymi wejściowymi zapisanymi w rzeczywistej tablicy f. Wyniki obliczeń mają zostać zapisane w zespolonej tablicy fftw_signal. Ciąg znaków r2c w nazwie wywoływanej funkcji określa właśnie, że wykonujemy transformatę real to complex, czyli z danych rzeczywistych do zespolonych. Długość wektora podana jest jako pierwszy argument funkcji i wynosi N+K (ilość próbek plus "zero padding"). Ostatni argument to tzw. flaga, która zwykle może przyjmować dwie wartości FFTW_MEASURE oraz FFTW_ESTIMATE. Pierwsza z nich wykonuje kilka krótkich testowych uruchomień FFT aby znaleźć najlepszą konfigurację do obliczenia żądanej transformaty. Zwykle zajmuje to trochę czasu. Druga flaga, użyta przez nas FFTW_ESTIMATE, nie zezwala na wykonywanie żadnych testowych uruchomień i tylko stara się estymować optymalną konfigurację do obliczeń.


Gdy plan działań zostanie już utworzony możemy go uruchomić wywołując:

fftw_execute(p);

I to właściwie cała wiedza o FFTW, która jest nam potrzebna do wykonania tego zadania. W podobny sposób dokonujemy transformaty tablicy G. Następnie wymnażamy te dwa wyniki w przestrzeni zespolonej (odpowiednio przeskalowując rezultat) i na końcu dokonujemy transformaty odwrotnej (w tym wypadku c2r czyli complex to real).


Kompilacja

Kompilacja naszego programu powinna na halo wyglądać następująco:

gcc -c -I/opt/fftw/3.1.2gnu/include/ main.c
gcc main.o -L/opt/fftw/3.1.2gnu/lib/ -lm -lfftw3 

utworzony plik a.out można testowo uruchomić sobie na halo np. w kolejce test.

Przykładowe wyniki

Poniżej prezentujemy przykładowe wyniki działania programu dla różnych wartości parametru sigma:

Fftw 01.png Fftw 00.png

Wykres funkcji Gaussa dla parametru \sigma=1 (zapisanej w stylu wrap-around) oraz wyniku splotu tej funkcji z funkcją f.

Fftw 11.png Fftw 10.png

Wykres funkcji Gaussa dla parametru \sigma=10 (zapisanej w stylu wrap-around) oraz wyniku splotu tej funkcji z funkcją f.

Fftw 21.png Fftw 20.png

Wykres funkcji Gaussa dla parametru \sigma=50 (zapisanej w stylu wrap-around) oraz wyniku splotu tej funkcji z funkcją f.

Fftw 31.png Fftw 30.png

Wykres funkcji Gaussa dla parametru \sigma=100 (zapisanej w stylu wrap-around) oraz wyniku splotu tej funkcji z funkcją f.


Przydatne linki

Ciekawostki: Fluent, superkomputery i Robert Kubica

Autor: Maciek Cytowski

Albert to komputer, który wykonuje obliczenia dla pojazdów drużyny BMW-Sauber (jednej ze stajni biorącej udział w tegorocznym Grand Prix Formuły 1, której kierowcą jest nasz rodak Robert Kubica). Komputer ten znajduje się obecnie na 60. miejscu listy TOP500. Głównym zadaniem komputera Albert jest testowanie nowych kształtów bolidu (m.in. przedniego spojlera). Obliczenia aerodynamiki, termodynamiki, mechaniki płynów wykonywane za pomocą oprogramowania ANSYS/Fluent są następnie sprawdzane na rzeczywistym modelu w tunelu aerodynamicznym.

Komputer Albert posiada 1024 rdzeni obliczeniowych oraz 2048 GB pamięci. Jego wydajność zmierzona za pomocą benchmarku Linpack to 9.6 TFlops. Żadna inna stajnia F1 nie posiada mocniejszej maszyny.

W roku 2007 w oparciu o wyniki obliczeń zaprezentowany został nowy kształt spojlerów.

Aerodynamika to niezwykle istotna (obok silnika, hamulców, opon itp. itd.) część bolida F1. Znane są przypadki nieudanych startów spowodowanych nie domknięciem klapki wlewu paliwa. W bolidzie Marka Webbera ze stajni Red Bull otwarta klapka wywoływała dodatkowe strumienie powietrza, które oddziaływały na górny tylny spojler i w rezultacie spowodowały jego złamanie. Podobna rzecz spotkała Polaka Roberta Kubicę w tegorocznym Grand Prix Bahrajnu, gdzie przez niedomknięty wlew zajął dopiero szóste miejsce.

Kształt i wytrzymałość części takich jak spojler przedni i tylny mają również kluczowe znaczenie dla bezpieczeństwa. W przedostatnim wyścigu Grand Prix w Kanadzie byliśmy świadkami potwornie wyglądającego wypadku Roberta Kubicy. Jego bolid uniósł się w powietrze przy prędkości 230 km/h na skutek oderwania się przedniego spojlera.