MN04: Różnice pomiędzy wersjami

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania
Przykry (dyskusja | edycje)
mNie podano opisu zmian
Przykry (dyskusja | edycje)
mNie podano opisu zmian
Linia 1: Linia 1:
==BLAS, LAPACK i ATLAS===


W zadaniach dotyczących macierzy gęstych, warto skorzystać z klasycznego tandemu
''Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi? poprawki''
bibliotek: BLAS (''Basic Linear Algebra Subprograms '') \cite{BLAS-home-page}
oraz LAPACK (''Linear Algebra PACKage '') \cite{LAPACK-home-page}. Dla macierzy
rozrzedzonych (zawierających wiele elementów zerowych) warto zastosować bardziej
wyspecjalizowane biblioteki. Aby wykorzystać do maksimum moc oferowaną przez te
biblioteki (w połączeniu z naszym komputerem) warto skorzystać z optymalizowanej
wersji BLASów i LAPACKa, czyli z ATLASa,  \cite{ATLAS-home-page}. Istnieje inna
wersja optymalizowanych BLASów, tzw. Goto BLAS. Niektóre procedury ATLASa są
istotnie szybsze od standardowych BLASów, dając, przykładowo dla zadania
mnożenia dwóch macierzy (na komputerze z procesorem Pentium 4 i
dostatecznie dużych macierzy), '''ponaddziesięciokrotne przyspieszenie''' na
zmiennych typu <code>float</code> i <code>double</code> i około pięciokrotne  na zmiennych
typu <code>complex </code> i <code>double complex</code>.


Aby osiągnąć największe przyspieszenie, bibliotekę ATLAS należy skompilować
==Ćwiczenia: równania nieliniowe skalarne==
samemu na własnej (''nieobciążonej'' w trakcie kompilacji!) architekturze. W plikach
<code>Makefile</code> ATLASa brak jednak opcji instalacji bibliotek w standardowych
miejscach --- trzeba zrobić to samemu.


BLAS i LAPACK są często wykorzystywane przez inne biblioteki numeryczne,
{{cwiczenie|Metoda Newtona może być zbieżna globalnie||
na nich opierają się również funkcje macierzowe w Octave i MATLABie.
Optymalizowane biblioteki BLAS i LAPACK dla swoich architektur oferują
producenci procesorów Intel (biblioteka MKL) oraz AMD (biblioteka ACML)


BLAS \cite{BLAS-home-page} jest kolekcją procedur służących manipulacji podstawowymi obiektami
Wykaż, że jeśli <math>\displaystyle f</math> jest rosnąca i wypukła na <math>\displaystyle [a,b]</math> oraz dla <math>\displaystyle x^*\in [a,b]\displaystyle f(x^*) = 0</math>, to metoda Newtona startująca z <math>\displaystyle x_0>x^*</math> jest zbieżna.
algebry liniowej: skalarami, wektorami i macierzami. Obsługiwane są obiekty
}}
rzeczywiste (w pojedynczej albo podwójnej precyzji) oraz zespolone (także w
dwóch precyzjach). W BLAS wyróżnia się 3 poziomy abstrakcji algorytmów:
* BLAS Level 1 -- działania typu wektor-wektor, np. operacja AXPY, czyli
uogólnione
dodawanie wektorów
<center><math>\displaystyle  
y \leftarrow \alpha x + y,
</math></center>


albo obliczanie normy wektora. BLAS Level 1 ma za zadanie porządkować kod;
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"> 
* BLAS Level 2 -- działania typu macierz--wektor, np. mnożenie macierzy
Powtarzając rachunki z dowodu lokalnej zbieżności, możemy łatwo stwierdzić, że kolejne
przez wektor
iteracje dają <math>\displaystyle x_k</math> które także spełniają warunek <math>\displaystyle x_k - x^* > 0</math>. Ponadto
<center><math>\displaystyle  
nietrudno sprawdzić, że <math>\displaystyle x_{k+1} < x_k</math>, a więc mamy ciąg malejący i
y \leftarrow \alpha A x + y
ograniczony. Jego granicą musi być <math>\displaystyle x^*</math> (dlaczego?).
</math></center>
</div></div>


Zapis algorytmów z użyciem BLAS Level 2 umożliwia potencjalnie przyspieszenie
{{cwiczenie|Fraktale||
programu, m.in. ze względu na to, że zoptymalizowane procedury BLAS Level 2 mogą np.
wykorzystywać instrukcje wektorowe nowoczesnych procesorów, zob.
Rozdział&nbsp;[[##sec:architektura|Uzupelnic: sec:architektura ]];
* BLAS Level 3 -- operacje typu macierz--macierz, np. mnożenie dwóch
macierzy:
<center><math>\displaystyle
C \leftarrow \alpha A\cdot B + C
</math></center>


W związku z tym, że operacje macierzowe wymagają wykonania <math>\displaystyle O(N^3)</math> działań
Ciekawy zbiór o charakterze fraktalnym powstaje w wyniku zastosowania metody
arytmetycznych przy
Newtona do rozwiązania równania <math>\displaystyle z^n-1=0</math> w dziedzinie zespolonej. Punkt <math>\displaystyle z_0 =
<math>\displaystyle O(N^2)</math> danych (gdzie <math>\displaystyle N</math> jest wymiarem macierzy), wykorzystanie
(x_0,y_0)</math> należy do ''basenu zbiezności metody'', jeśli metoda Newtona jest
zoptymalizowanych procedur BLAS Level 3 może znacząco przyspieszyć wykonywanie
zeń zbieżna do (jakiegokolwiek) pierwiastka w/w równania. Kolor odpowiadający
obliczeń na maszynach z pamięcią hierarchiczną (zob. Rozdział&nbsp;[[##sec:cache|Uzupelnic: sec:cache ]]).
<math>\displaystyle z_0</math> jest określany na podstawie liczby iteracji potrzebnych metodzie do
zbieżności.


Podsumowując: przyjmijmy, że dysponujemy dobrze zoptymalizowaną biblioteką BLAS.
Zupełnie miłym (i estetycznie wartościowym) doświadczeniem, jest napisanie
Wówczas dany algorytm algebry liniowej najlepiej zapisać przy użyciu procedur  BLAS Level
programu w Octave, który wyświetla baseny zbieżności metody Newtona jak poniżej.
3, naturalnie, pod warunkiem, że w ogóle daje się to zrobić; typową strategią w
takim wypadku jest tworzenie algorytmów blokowych, operujących na ''blokach''
macierzy, zamiast na jej pojedynczych elementach.


Większość użytkowników BLAS nie będzie jednak miała potrzeby pisania własnych
[[Image:MNfraktal.png|400px|Basen zbiezności metody Newtona w okolicy początku układu
algorytmów blokowych, gdyż funkcje rozwiązujące podstawowe zadania numerycznej
współrzędnych, dla równania <math>\displaystyle z^5-1=0</math>]]
algebry liniowej: rozwiązywanie układów równań (także nad- i niedookreślonych)
oraz zadania własnego, znajdują się w doskonałym pakiecie LAPACK
\cite{LAPACK-home-page}, który
intensywnie i skutecznie wykorzystuje podprogramy BLAS.


Nazwy procedur BLASów i
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
LAPACKa są cokolwiek enigmatyczne na pierwszy rzut oka, ale w istocie bardzo
Bardzo wygodnie napisać taki program w postaci skryptu Octave. Pamiętaj
łatwo je odgadywać. Każda nazwa składa się z kilku części, najczęściej jest
jednak, by uniknąć w nim pętli w rodzaju
postaci <code>PRRFF</code>, gdzie
;
<code>P</code> oznacza precyzję i może przyjmować
wartości: S,D,C,Z, odpowiadające kolejno pojedynczej  i podwójnej precyzji w
dziedzinie
rzeczywistej i pojedynczej  i podwójnej precyzji w dziedzinie zespolonej;


;
for xpixels <nowiki>=</nowiki> 1:1024<br>
<code>RR</code> oznacza rodzaj zadania, np. GE oznacza ''GEneral '', czyli zadanie ogólne
for ypixels <nowiki>=</nowiki> 1:768<br>
(praktycznie bez założeń), a SY oznacza ''SYmmetric '', czyli zadanie symetryczne;
... wyznacz kolor pixela ....<br>
end<br>
end<br>


;
Taka podwójnie zagnieżdżona pętla może skutecznie zniechęcić Cię do
:  <code>FF</code> wreszcie określa samo zadanie, np. SV oznacza
Octave'a --- obrazek będzie liczyć się wieki --- zupełnie niepotrzebnie!
''SolVe ''(w domyśle: układ równań), MV --- ''Matrix-Vector ''(w domyśle: mnożenie), 
EV --- ''EigenValues '', czyli wartości własne, itp. Są też warianty
trzyliterowe, np. TRF (''TRiangular Factorization '') i TRS  (''TRiangular
Solve ''--- w domyśle, przy użyciu wcześniej wyznaczonej faktoryzacji)
Jeśli jednak ''nie możemy zgadnąć'', jaka jest nazwa procedury BLAS/LAPACKa,
która byłaby nam potrzebna,
najlepiej przejrzeć (przeszukać) strony internetowe tych pakietów w serwisie
Netlib.


Zestawienie najczęściej wykorzystywanych procedur BLAS i LAPACKa przedstawiamy
</div></div>
poniżej. Każda z tych procedur ma swój wariant "ekspercki", np. dla
rozwiązywania układów równań metodą eliminacji Gaussa można skorzystać z
osobnych procedur generujących rozkład LU oraz z innych, rozwiązujących układy
trójkątne.


{| border=1
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
|+ <span style="font-variant:small-caps">Uzupelnij tytul</span>
Kluczem do efektywnej implementacji w Octave jest przeprowadzenie
|-
iteracji Newtona na ''wszystkich pikselach jednocześnie''. W tym celu musisz
|
skorzystać z prowadzenia działań na całych macierzach.
</div></div>


Zadanie algebry liniowej  ||  Nazwa procedury BLAS/LAPACK
}}
|-
|


mnożenie wektora przez macierz  || DGEMV
{{cwiczenie|Pierwiastkowanie||
|-
| mnożenie macierzy przez macierz  ||  DGEMM
|-
|
rozwiązywanie układu równań  ||  DGESV
|-
| rozkład LU (w miejscu)  ||  DGETRF
|-
| rozwiązywanie układu z rozkładem uzyskanym  z DGETRF  ||  DGETRS
|-
|
rozwiązywanie układu z macierzą symetryczną  ||  DSYSV
|-
| rozkład LDL<math>\displaystyle ^T</math> macierzy symetrycznej (w miejscu)  ||  DSYTRF
|-
| rozwiązywanie układu z rozkładem uzyskanym z DSYTRF  ||  DSYTRS
|-
|
rozwiązywanie układu z macierzą pasmową  ||  DGBSV
|-
| rozkład LU macierzy pasmowej (w miejscu)  ||  DGBTRF
|-
| rozwiązywanie układu z rozkładem uzyskanym  z DGBTRF  ||  DGBTRS
|-
|
zagadnienie własne  ||  DGESV
|-
|  


|}
Niech <math>\displaystyle a\geq 0</math>. Aby wyznaczyć <math>\displaystyle \sqrt{a}</math>, można skorzystać z metody Newtona dla
równania <math>\displaystyle x^2 = a</math>. Zaprogramuj tę metodę i sprawdź, jak wiele dokładnych cyfr
wyniku dostajesz w kolejnych krokach. Czy to możliwe, by liczba dokładnych cyfr
znaczących z grubsza podwajała się na każdej iteracji? Wskaż takie <math>\displaystyle a</math> dla
którego to nie będzie prawdą.


====Mnożenie macierz-wektor w BLAS====
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
Na pewno kusi Cię, by zaprogramować najpierw ogólnie funkcję "metoda
Newtona", a następnie przekazać jej jako parametr naszą funkcję. To oczywiście
będzie działać, i to praktycznie tak samo szybko jak drugi spsoób, w którym po prostu
wyliczysz na karteczce wzór na iterację dla ''tej konkretnej funkcji, <math>\displaystyle F(x) =
x^2-a</math>'' i potem go zaimplementujesz. Jednak
drugie podejście jest tutaj godne polecenia, bo ma na sobie znak charakterystyczny dla numeryki: niech
wszystko, co programujesz będzie tak proste, jak tylko możliwe (ale nie
prostsze)!
</div></div>


Zacznijmy od prostej procedury BLAS Level 2, jaką jest mnożenie macierzy przez
}}
wektor. Wykorzystamy tzw. wzorcową implementację BLASów (niezoptymalizowaną)
dostarczaną z dystrybucją np. Red Hat Linux. Jest to  biblioteka funkcji
fortranowskich.


Naszym zadaniem jest wykonanie operacji
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"> 
Nasza iteracja, którą chcemy zaprogramować, to
<center><math>\displaystyle  
<center><math>\displaystyle  
y \leftarrow  \alpha A x + y,
x_{k+1} = \frac{1}{2}\left( x_k + \frac{a}{x_k}\right).
</math></center>
 
Dla <math>\displaystyle a\neq 0</math>, spełnione są założenia twierdzenia o zbieżności metody Newtona,
więc
 
<center><math>\displaystyle |x_{k+1}-x^*|\leq K |x_{k}-x^*|^2
</math></center>
</math></center>


gdzie <math>\displaystyle A</math> jest zadaną macierzą <math>\displaystyle N\times M</math>, natomiast <math>\displaystyle y</math> jest wektorem o <math>\displaystyle M</math>
co właśnie tłumaczy się jako podwajanie liczby dokładnych cyfr
współrzędnych.
znaczących wyniku.
 
Ale jeśli <math>\displaystyle a=0</math> to metoda zbieżna jest już tylko liniowo (z ilorazem 0.5 --- bo
jaki wtedy jest wzór na iterację?), więc co mniej więcej trzy kroki będziemy
dostawali kolejne zero dokładnego wyniku <math>\displaystyle x^*=0</math>.
</div></div>


To zadanie realizuje procedura BLASów o nazwie
{{cwiczenie|Odwrotność bez dzielenia||
<code>DGEMV</code>. W rzeczywistości, ta procedura wykonuje ogólniejsze zadanie
wyznaczania wektora


Aby wyznaczyć <math>\displaystyle 1/a</math> dla <math>\displaystyle a>0</math> bez dzielenia(!), można zastosować metodę Newtona
do funkcji <math>\displaystyle f(x) = 1/x - a</math>. Pokaż, że na <math>\displaystyle k</math>-tym kroku iteracji,
<center><math>\displaystyle  
<center><math>\displaystyle  
y \leftarrow  \alpha B x + \beta y,
|ax_k - 1| = |ax_{k-1} - 1|^2.
</math></center>
</math></center>


przy czym macierz <math>\displaystyle B</math> może być równa albo <math>\displaystyle A</math>, albo <math>\displaystyle A^T</math> (jednak za każdym
Dla jakich <math>\displaystyle x_0</math> metoda będzie zbieżna do <math>\displaystyle \frac{1}{a}</math>, a dla jakich nie?
razem argumentem macierzowym, jaki przekazujemy <code>DGEMV</code>, jest wyjściowa
Oceń, ile iteracji potrzeba do spełnienia warunku <math>\displaystyle \frac{|x_k-\frac{1}{a}|}{|a|} \leq
macierz <math>\displaystyle A</math>).
\epsilon</math>, gdy <math>\displaystyle \frac{|x_0-\frac{1}{a}|}{|a|} \leq \gamma < 1</math>,
}}


Jak wiemy, że jest możliwe bezpośrednie wykorzystanie
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"> 
biblioteki fortranowskiej w programie w C, jednak musimy pamiętać, iż macierze z
Proste ćwiczenie z analizy. Warunkiem zbieżności jest oczywiście (ze względu na
jakich ona skorzysta muszą być ułożone ''kolumnami'' w jednolitym bloku
równość w wyrażeniu na błąd iteracji powyżej) <math>\displaystyle |ax_0 - 1| < 1</math>, to znaczy
pamięci.
<center><math>\displaystyle 0 < x_0 < \frac{2}{a}.</math></center>


Bazując na opisie procedury <code>DGEMV</code> ze
</div></div>
strony \pageref{opis:dgemv}, w programie w C powinniśmy
napisać prototyp tej funkcji następująco:


<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
{{cwiczenie|||
Zaimplementuj metodę bisekcji. Sprawdź, jak będzie działać m.in. dla funkcji
* <math>\displaystyle f(x) = \sin(x)</math>,
* <math>\displaystyle f(x) = \sin^2(x)</math>,
* <math>\displaystyle f(x) = (x-1)^5</math> (wzór A),
* <math>\displaystyle f(x) = (x-1)*(x^4 - 4x^3 + 6x^2 - 4x + 1)</math> (wzór B),
* <math>\displaystyle f(x) = (x-2)^13</math> (wzór C),
* <math>\displaystyle f(x) = x^13 - 26*x^12 + 312*x^11 - 2288*x^10 + ... - 8192</math> (wzór D),
   
   
int dgemv_( char* TRANS,  
gdy tolerancję błędu zadasz na poziomie <math>\displaystyle 10^{-10}</math>.
int* M,  
 
int* N,  
Jak wyjaśnić te wyniki? Czy możesz już być pewna, że masz dobrą implementację?
double* ALPHA,  
 
double* A,
}}
int* LDA,
 
double* X,
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"> 
int* INCX,
 
double* BETA,
Jeśli nie pomylisz się, metoda powinna zbiegać bez problemów do zera funkcji
double* Y,  
<math>\displaystyle \sin</math>, dawać komunikat o błędzie dla <math>\displaystyle \sin^2</math> (bo nie ma przeciwnych znaków).
int* INCY );
 
</pre></div>
Zapewne zauważyłaś, że wzory A i B są matematycznie równoważne, podobnie zresztą
 
jak wzory C i D. Niestety, tylko wzór A i C dają w efekcie dobre przybliżenia
Dla własnej wygody, a także dla przyszłego wykorzystania, umieścimy ten
miejsca zerowego, podczas gdy wzory B i D prowadzą do oszacowania na miejsce
prototyp, razem z innymi przydatnymi drobiazgami (m.in. makro <code>IJ</code>
zerowe, które ''w ogóle nie zawierają'' prawdziwego miejsca zerowego.
dające wygodny dostęp do macierzy niezależny od jej wewnętrznej reprezentacji, a
 
także zmienne całkowite
Wyjaśnieniem tego fenomenu jest zjawisko redukcji cyfr przy odejmowaniu.
<code>static int BLASONE <nowiki>=</nowiki> 1, BLASMONE <nowiki>=</nowiki> -1;</code>), w pliku
 
nagłówkowym <code>blaslapack.h</code>.
[[Image:MNbisekcjawblad.png|400px|Metoda bisekcji ma kłopoty, gdy funkcja zadana jest
wzorem D.]]
[[Image:MNbisekcjawresid.png|400px|Residuum też jest duże, gdy <math>\displaystyle f</math> zadana jest
wzorem D.]]
 
Jeśli chodzi o pewność... No cóż, sprawdziłaś, że działa w przypadkach, gdy
spodziewałaś się, że będzie działać. Że tam, gdzie spodziewałaś się kłopotów,
lub komunikatu o błędzie --- tak rzeczywiście się stało. Wreszcie,
potwierdziłaś, że zachowanie metody jest zgodne z jej znanymi właściwościami.
 
Tak więc, można ''przypuszczać'', że implementacja jest poprawna.
 
</div></div>
 
{{cwiczenie|||
Wskaż ''wszystkie'' wartości <math>\displaystyle x_0</math>, dla jakich metoda Newtona będzie zbieżna
do rozwiązania <math>\displaystyle x^*</math> równania
 
<center><math>\displaystyle \arctan(x) = 0.
</math></center>


Dalej już jest łatwo: oto pełny kod programu realizującego operację mnożenia macierzy przez wektor
Wyznacz wartość <math>\displaystyle X_0</math>, z którego startując powinieneś dostać ciąg oscylujący <math>\displaystyle X_0, -X_0,\ldots</math>.
przy użyciu procedury BLAS <code>DGEMV</code>:
Sprawdź eksperymentalnie, czy tak jest rzeczywiście.


<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Wskazówka </span><div class="mw-collapsible-content" style="display:none">
#include <stdio.h>
Aby znaleźć graniczne <math>\displaystyle X_0</math>, czyli takie, dla którego dostaniesz
#include "blaslapack.h"
oscylacje <math>\displaystyle X_0, -X_0,\ldots</math>,  musisz rozwiązać równanie nieliniowe:


double* mmread(char *filename, int* N, int* M );
<center><math>\displaystyle -X_0 = X_0 - \frac{\arctan(X_0)}{\frac{1}{1+X_0^2}}.
</math></center>


int main()
To łatwe.
{
int N, M, i, j;
double *A, *x, *y;
/* wczytujemy macierz z pliku w formacie MatrixMarket */
/* macierz jest reprezentowana w formacie kolumnowym  */
A <nowiki>=</nowiki> mmread("mbeacxc.mtx", &N, &M );
x <nowiki>=</nowiki> (double *)malloc(N*sizeof(double));
y <nowiki>=</nowiki> (double *)malloc(M*sizeof(double));
for (i <nowiki>=</nowiki> 1; i <<nowiki>=</nowiki> N; i++)
x[IJ(i,1,N)] <nowiki>=</nowiki> (double)i;


/* obliczamy y <nowiki>=</nowiki> 5*A*x, korzystając z procedury BLAS L2: DGEMV */
</div></div>


{
}}
char TRANS <nowiki>=</nowiki> 'N'; double ALPHA <nowiki>=</nowiki> 5.0, BETA <nowiki>=</nowiki> 0.0;
dgemv_(&TRANS, &N, &M, &ALPHA, A,  &N,  x,  &BLASONE,
                        &BETA, y, &BLASONE );


}
<div class="mw-collapsible mw-made=collapsible mw-collapsed"><span class="mw-collapsible-toogle mw-collapsible-toogle-default style="font-variant:small-caps">Rozwiązanie </span><div class="mw-collapsible-content" style="display:none"> 
Dla <math>\displaystyle -X_0 < x_0 < X_0</math> mamy zbieżność. Dla <math>\displaystyle |x_0| = |X_0|</math> oscylacje, dla
/* wydruk wyniku */
większych wartości --- rozbieżność. Jednak w eksperymencie wszystko zależy od
for (i <nowiki>=</nowiki> 1; i <<nowiki>=</nowiki> M; i++)
tego, jak dokładnie wyznaczyliśmy <math>\displaystyle X_0</math>. Na przykład, mój solver w Octave
printf("
wyznaczył wartość <math>\displaystyle X_0</math>, dla której metoda Newtona dała następujący ciąg:
return(0);
<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
}
</pre></div>
   
   
Zwróćmy uwagę na wykorzystanie "stałej" <code>BLASONE</code>, równej 1,
Iteracja:  1,  x <nowiki>=</nowiki> -1.39175
predefiniowanej w pliku <code>blaslapack.h</code>. Nasz program kompilujemy
Iteracja:  2,   x <nowiki>=</nowiki> 1.39175
standardowo, pamiętając o dołączeniu na etapie linkowania używanych przez nas
Iteracja:  3,   x <nowiki>=</nowiki> -1.39175
bibliotek:
Iteracja:  4,  x <nowiki>=</nowiki> 1.39175
Iteracja:  5,  x <nowiki>=</nowiki> -1.39175
 
.. itd ..


<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
Iteracja:  25,  x <nowiki>=</nowiki> -1.39176
   
Iteracja:  26,  x <nowiki>=</nowiki> 1.39178
gcc -o testblas testblas.c -llapack -lblas -lgfortran -lm
Iteracja:  27,  x <nowiki>=</nowiki> -1.39183
Iteracja:  28,  x <nowiki>=</nowiki> 1.39197
Iteracja:  29,  x <nowiki>=</nowiki> -1.39235
Iteracja:  30,  x <nowiki>=</nowiki> 1.39333
Iteracja:  31,  x <nowiki>=</nowiki> -1.39594
Iteracja: 32,  x <nowiki>=</nowiki> 1.40283
Iteracja: 33,  x <nowiki>=</nowiki> -1.42117
Iteracja:  34,  x <nowiki>=</nowiki> 1.47061
Iteracja:  35,  x <nowiki>=</nowiki> -1.60867
Iteracja: 36,  x <nowiki>=</nowiki> 2.03161
Iteracja:  37,  x <nowiki>=</nowiki> -3.67722
Iteracja:  38,  x <nowiki>=</nowiki> 15.2779
Iteracja:  39,  x <nowiki>=</nowiki> -337.619
Iteracja:  40,  x <nowiki>=</nowiki> 178376
Iteracja:  41,  x <nowiki>=</nowiki> -4.99792e+10
Iteracja:  42,  x <nowiki>=</nowiki> 3.92372e+21
Iteracja:  43,  x <nowiki>=</nowiki> -2.41834e+43
Iteracja:  44,  x <nowiki>=</nowiki> 9.18658e+86
Iteracja:  45,  x <nowiki>=</nowiki> -1.32565e+174
Iteracja:  46,  x <nowiki>=</nowiki> inf
Iteracja:  47,  x <nowiki>=</nowiki> nan
</pre></div>
</pre></div>
   
   
--- dokładnie ''w tej właśnie kolejności'' (LAPACK oczywiście w tym momencie
</div></div>
dołączamy na wyrost: nasz program nie korzysta z żadnej funkcji LAPACKa, wobec
tego opcja <code>-llapack</code> zostanie zignorowana).
 
Pamiętamy oczywiście, że standardowe BLASy i LAPACK nie są zoptymalizowane w
stopniu pozwalającym na (prawie) maksymalne wykorzystanie możliwości sprzętu.
Dla osiągnięcia maksymalnej efektywności kodu, trzeba skorzystać z
optymalizowanych BLAS, które obecnie są dostępne nawet w kilku wariantach na
architektury x86.

Wersja z 16:33, 28 sie 2006

Uwaga: przekonwertowane latex2mediawiki; prawdopodobnie trzeba wprowadzi? poprawki

Ćwiczenia: równania nieliniowe skalarne

Ćwiczenie Metoda Newtona może być zbieżna globalnie

Wykaż, że jeśli f jest rosnąca i wypukła na [a,b] oraz dla x*[a,b]f(x*)=0, to metoda Newtona startująca z x0>x* jest zbieżna.

Rozwiązanie

Ćwiczenie Fraktale

{{{3}}}

Ćwiczenie Pierwiastkowanie

{{{3}}}
Rozwiązanie

Ćwiczenie Odwrotność bez dzielenia

Aby wyznaczyć 1/a dla a>0 bez dzielenia(!), można zastosować metodę Newtona do funkcji f(x)=1/xa. Pokaż, że na k-tym kroku iteracji,

|axk1|=|axk11|2.

Dla jakich x0 metoda będzie zbieżna do 1a, a dla jakich nie? Oceń, ile iteracji potrzeba do spełnienia warunku |xk1a||a|ϵ, gdy |x01a||a|γ<1,

Rozwiązanie

Ćwiczenie

Zaimplementuj metodę bisekcji. Sprawdź, jak będzie działać m.in. dla funkcji

  • f(x)=sin(x),
  • f(x)=sin2(x),
  • f(x)=(x1)5 (wzór A),
  • f(x)=(x1)*(x44x3+6x24x+1) (wzór B),
  • f(x)=(x2)13 (wzór C),
  • f(x)=x1326*x12+312*x112288*x10+...8192 (wzór D),

gdy tolerancję błędu zadasz na poziomie 1010.

Jak wyjaśnić te wyniki? Czy możesz już być pewna, że masz dobrą implementację?

Rozwiązanie

Ćwiczenie

{{{3}}}
Rozwiązanie