GSL
Re: GSL
Ni temu namenjena. To funkcijo ima blas (v vseh najrazlicnejsih implementacijah, ki imajo standardizirani in dobro poznan nabor funkcij, ki pokrijejo vse potrebne operacije). gsl prinese svoj blas: linkas ga preko -lgslcblas, lahko pa seveda uporabis katerokoli implementacijo.
Re: GSL
Ok, za GSL nisem našel rešitve, sem pa poskusil dve različni sinsusni fft (ena iz NR druga nekje na netu pobrana) pa še matlabovo in sem primerjal med seboj. No, vse tri so dale različne rezultate! Bi pa poudaril, da je tista iz NR brezveze, ker zna računati zgolj za \(n = 2^k\) elementov.
Re: GSL
Za tisto vprasanje prej:
http://www.gnu.org/software/gsl/manual/ ... brary.html
http://www.gnu.org/software/gsl/manual/ ... S-Examples
http://en.wikipedia.org/wiki/BLAS
Matlab ni recimo nic drugega kot zakamufliran blas in lapack (pa se par dodatkov, ampak osnovna funkcionalnost z matrikami je tukaj).
fft algoritem je po definiciji dizajniran za n=2^k elementov, saj je v tem vir njegove pohitritve: deluje na principu, da razdeli problem na dva sorodna podproblema, ki ju lahko resujes hkrati, in to rekurzivno naredis k-krat, in je zato casovna zahtevnost n*k=n*log(n). Za razliko od naivne transformacije, ki je n^2 (in zato neuporabna). Ce hoces delat na n-ju, ki ni potenca stevila 2, imas dve moznosti. 1) zafilas z niclami do naslednje 2^k in transformiras, pri cemer moras potem pazit da rezultat ni enak kot bi bil sicer - za nekatere rabe to ni ustrezno. 2) posplosen/nadgrajen fft algoritem lahko razcepi stevilo na prafaktorje, in zna rekurzivno razbit tudi za faktorje, razlicne od 2. Za te faktorje je potem pohitritev manjsa. In za velike faktorje gre potem ponavadi vsaka implementacija nazaj na ta tezko racunanje.
Omejitev na 2^k je v bistvu prednost - vodilo, da programer pac svoje probleme diskretizira na ta nacin, da to prednost izkoristi. 2^k je tudi zaradi drugih programerskih razlogov dobra ideja, in tudi za algoritme, ki niso omejeni na 2^k, bo delalo hitreje, in enostavneje bo programirat, ce prilagodis izbiro vzorca na 2^k. Tako da tudi, ce imas fft knjiznico, ki je sposobna prezvecit razlicne dolzine, se je najbolje drzat potenc stevila dva - je zdravo in udobno.
Ena in edina optimalna izbira knjiznice je fftw (verzija fftw3), ki je optimizirana za hitrost, in ce me spomin ne vara, ima pohitritev implementirano do stevil s faktorji do 17 (ceprav so najhitrejsi se vedno 2^k). Mocno priporocam. GSL ni za fft misljen. Navadi se, da imas specializirane in optimizirane knjiznice za dolocene tipe problemov - blas za osnovno mnozenje matrik in vektorjev in tega, lapack/arpack/atlas za bolj napredne matricne zadeve (resevanje sistemov, razcepi in obracanje matrik - QR, LU, SVD, lastne vrednosti,...), fftw za fourierovo transformacijo, GSL pa potem pokrije bolj numericne algoritme (specialne funkcije, diferencialne enacbe, minimizacija, iskanje nicel,..). Ni nekega razloga, da bi gsl izumljal toplo vodo za stvari, ki so 100% optimizirane ze 30 let.
NR je seveda referencna implementacija... NR je ucbenik za ucenje algoritmov, ne pa knjiznica, ki bi jo dejansko kdo uporabljal v praksi (tudi nevarno jo je uporabljat, in kvari programerski stil, se posebej ce vzames tisto varianto za C, ki je prepisana iz fortrana in ima indeksiranje polomljeno).
Vprasanje je tudi, ce si prav klical fft. Pazi na razporeditev elementov na vhodu, in na izhodu (kateri indeksi so katere frekvence). Kako to mislis da je rezultat razlicen? Je razlika na decimalkah, ali je oblika spektra napacna? Vprasanje je tudi pri interpretaciji vhodnih podatkov. Recimo pri fftw bos videl, da imas za kosinusno in sinusno transformacijo vec moznosti kako je periodicno razsirjeno (ali je simetrijska tocka na prvi tocki, ali pol tocke v levo...).
Sicer iz fft transformiranke ni problem dobit tega... saj ves da je realni del komponent kosinus, imaginarni pa sinus. Ce podatke pred transformacijo sodo razsiris, bodo pa imaginarne komponente itak 0. Podobno, ce jih liho razsiris in dobis sinusno transformacijo.
http://www.gnu.org/software/gsl/manual/ ... brary.html
http://www.gnu.org/software/gsl/manual/ ... S-Examples
http://en.wikipedia.org/wiki/BLAS
Matlab ni recimo nic drugega kot zakamufliran blas in lapack (pa se par dodatkov, ampak osnovna funkcionalnost z matrikami je tukaj).
fft algoritem je po definiciji dizajniran za n=2^k elementov, saj je v tem vir njegove pohitritve: deluje na principu, da razdeli problem na dva sorodna podproblema, ki ju lahko resujes hkrati, in to rekurzivno naredis k-krat, in je zato casovna zahtevnost n*k=n*log(n). Za razliko od naivne transformacije, ki je n^2 (in zato neuporabna). Ce hoces delat na n-ju, ki ni potenca stevila 2, imas dve moznosti. 1) zafilas z niclami do naslednje 2^k in transformiras, pri cemer moras potem pazit da rezultat ni enak kot bi bil sicer - za nekatere rabe to ni ustrezno. 2) posplosen/nadgrajen fft algoritem lahko razcepi stevilo na prafaktorje, in zna rekurzivno razbit tudi za faktorje, razlicne od 2. Za te faktorje je potem pohitritev manjsa. In za velike faktorje gre potem ponavadi vsaka implementacija nazaj na ta tezko racunanje.
Omejitev na 2^k je v bistvu prednost - vodilo, da programer pac svoje probleme diskretizira na ta nacin, da to prednost izkoristi. 2^k je tudi zaradi drugih programerskih razlogov dobra ideja, in tudi za algoritme, ki niso omejeni na 2^k, bo delalo hitreje, in enostavneje bo programirat, ce prilagodis izbiro vzorca na 2^k. Tako da tudi, ce imas fft knjiznico, ki je sposobna prezvecit razlicne dolzine, se je najbolje drzat potenc stevila dva - je zdravo in udobno.
Ena in edina optimalna izbira knjiznice je fftw (verzija fftw3), ki je optimizirana za hitrost, in ce me spomin ne vara, ima pohitritev implementirano do stevil s faktorji do 17 (ceprav so najhitrejsi se vedno 2^k). Mocno priporocam. GSL ni za fft misljen. Navadi se, da imas specializirane in optimizirane knjiznice za dolocene tipe problemov - blas za osnovno mnozenje matrik in vektorjev in tega, lapack/arpack/atlas za bolj napredne matricne zadeve (resevanje sistemov, razcepi in obracanje matrik - QR, LU, SVD, lastne vrednosti,...), fftw za fourierovo transformacijo, GSL pa potem pokrije bolj numericne algoritme (specialne funkcije, diferencialne enacbe, minimizacija, iskanje nicel,..). Ni nekega razloga, da bi gsl izumljal toplo vodo za stvari, ki so 100% optimizirane ze 30 let.
NR je seveda referencna implementacija... NR je ucbenik za ucenje algoritmov, ne pa knjiznica, ki bi jo dejansko kdo uporabljal v praksi (tudi nevarno jo je uporabljat, in kvari programerski stil, se posebej ce vzames tisto varianto za C, ki je prepisana iz fortrana in ima indeksiranje polomljeno).
Vprasanje je tudi, ce si prav klical fft. Pazi na razporeditev elementov na vhodu, in na izhodu (kateri indeksi so katere frekvence). Kako to mislis da je rezultat razlicen? Je razlika na decimalkah, ali je oblika spektra napacna? Vprasanje je tudi pri interpretaciji vhodnih podatkov. Recimo pri fftw bos videl, da imas za kosinusno in sinusno transformacijo vec moznosti kako je periodicno razsirjeno (ali je simetrijska tocka na prvi tocki, ali pol tocke v levo...).
Sicer iz fft transformiranke ni problem dobit tega... saj ves da je realni del komponent kosinus, imaginarni pa sinus. Ce podatke pred transformacijo sodo razsiris, bodo pa imaginarne komponente itak 0. Podobno, ce jih liho razsiris in dobis sinusno transformacijo.
Re: GSL
Hvala, sicer sem kar sam napisal zanke za množenje matrik z vektorjem in podobno. Tu najbrž itak ni prostora za nadaljno optimizacijo.Aniviller napisal/-a:Za tisto vprasanje prej:
http://www.gnu.org/software/gsl/manual/ ... brary.html
http://www.gnu.org/software/gsl/manual/ ... S-Examples
http://en.wikipedia.org/wiki/BLAS
Maš pa menda rutine za QR, Cholesky, ki pa so prilagojene specialnim matrikam (redke in slično).Aniviller napisal/-a:Ena in edina optimalna izbira knjiznice je fftw (verzija fftw3), ki je optimizirana za hitrost, in ce me spomin ne vara, ima pohitritev implementirano do stevil s faktorji do 17 (ceprav so najhitrejsi se vedno 2^k). Mocno priporocam. GSL ni za fft misljen. Navadi se, da imas specializirane in optimizirane knjiznice za dolocene tipe problemov - blas za osnovno mnozenje matrik in vektorjev in tega, lapack/arpack/atlas za bolj napredne matricne zadeve (resevanje sistemov, razcepi in obracanje matrik - QR, LU, SVD, lastne vrednosti,...), fftw za fourierovo transformacijo, GSL pa potem pokrije bolj numericne algoritme (specialne funkcije, diferencialne enacbe, minimizacija, iskanje nicel,..). Ni nekega razloga, da bi gsl izumljal toplo vodo za stvari, ki so 100% optimizirane ze 30 let.
Aniviller napisal/-a:NR je seveda referencna implementacija... NR je ucbenik za ucenje algoritmov, ne pa knjiznica, ki bi jo dejansko kdo uporabljal v praksi (tudi nevarno jo je uporabljat, in kvari programerski stil, se posebej ce vzames tisto varianto za C, ki je prepisana iz fortrana in ima indeksiranje polomljeno).
To je najbrž tista varianta, ki štarta indeksiranje z 1. Uporabljam sicer verzijo C++, a ne tiste, ki ima classe in podobno.Pri sinusni foruirerovi transformaciji, potegnjeni s spleta, je razlika predvsem na decimalkah. Pri NR pa predznakih in celih številih.Aniviller napisal/-a:Vprasanje je tudi, ce si prav klical fft. Pazi na razporeditev elementov na vhodu, in na izhodu (kateri indeksi so katere frekvence). Kako to mislis da je rezultat razlicen? Je razlika na decimalkah, ali je oblika spektra napacna? Vprasanje je tudi pri interpretaciji vhodnih podatkov. Recimo pri fftw bos videl, da imas za kosinusno in sinusno transformacijo vec moznosti kako je periodicno razsirjeno (ali je simetrijska tocka na prvi tocki, ali pol tocke v levo...).
Dokler se nisem ukvarjal s fourierovo transformacijo pri numeriki, sem bil prepričan, da je tako. Ampak imaginarni del fft ni enako dst.Aniviller napisal/-a:Sicer iz fft transformiranke ni problem dobit tega... saj ves da je realni del komponent kosinus, imaginarni pa sinus. Ce podatke pred transformacijo sodo razsiris, bodo pa imaginarne komponente itak 0. Podobno, ce jih liho razsiris in dobis sinusno transformacijo.
Re: GSL
Tukaj je VELIKO prostora za optimizacijo. Ko gre enkrat za optimizacijo tako trivialnih procesov kot je mnozenje in sestevanje, se pozna vsaka neumnost, in je vitalnega pomena znanje o delovanju procesorja. Vrstni red kako gres po matriki recimo vpliva na to, kako ucinkovita je raba Cache-ja in kolikokrat procesor ustavi in caka na dostavo podatkov iz RAM-a. Ze samo to lahko pridela nekaj deset procentov, ce pazis. Vprasanje je recimo tudi ali najprej jih par zmnozis, in spravis v lokalne spremenljivke in na koncu sestejes, ali delas sproti,... Da ne govorimo o uporabi posebnih procesorskih instrukcij za namen sestevanja (imas instrukcije, ki naredijo ax+b v eni potezi - mnozenje in sestevanje), pa vzporedno racunanje (moderni intel in amd imajo vsi sse (sse2,3,...) instrukcije, ki dovoljujejo mnozenje / sestevanje po 4 pare stevil hkrati, kar sicer prevajalnik mogoce zna sam pogruntat, ampak v praksi je cloveski poseg boljsi)... potem imas pa seveda se algoritme, ki mnozijo matrike tako, da naredijo manj mnozenj kot bi bilo naivno treba: glej recimo Strassenov algoritem http://en.wikipedia.org/wiki/Strassen_algorithm ki spravi mnozenje 2x2 matrik iz 8 na 7 mnozenj (in spet dela rekurzivno, kar pomeni da je spet blazno boljse imet matrike dimenzij 2^k).Hvala, sicer sem kar sam napisal zanke za množenje matrik z vektorjem in podobno. Tu najbrž itak ni prostora za nadaljno optimizacijo.
Stvar je tako kriticna (posebej za ljudi, ki imajo opravka z gromozanskimi matrikami in programi tecejo po vec dni), da pri prevajanju knjiznic, ki implementirajo blas funkcije, zenejo poskusanje in kalibracijo na tocno dolocen racunalnik (kombinacija ram, procesor, maticna plosca - vse hitrosti in odzivni casi), poskusa prevest na precej razlicnih nacinov in izmeri case, in uporabi najhitrejsega. ATLAS recimo pocne tocno to - poskusa zagotovit, da dobis blas, ki na danem racunalniku dela kar se da hitro. Ce dobis 15% pohitritev se ze kar pozna na hitrosti, ce imas res kriticne racune.
Ja, nekaj tega imajo ja. Sicer je pa tudi za sparse matrike par knjiznic ki so specializirane za to. Imas sparse BLAS, pa arpack recimo je za sparse lastne vrednosti, in tako naprej... GSL je novejsa zadeva kot ostale knjiznice, ampak nekako se vedno izpade kot nekoliko bolj amaterski - dober zacetek, ampak ni tako razsirjen v profesionalni rabi. Je pa zelo nepogresljiv na podrocju minimizacij, iskanja nicel, resevanja diferencialnih enacb,... ker to je nekako tisto kar je tezje najt v drugih knjiznicah. Nekako lahko reces, da ima GSL precejsen del NR-ja, ampak narejen na pravi nacin.Maš pa menda rutine za QR, Cholesky, ki pa so prilagojene specialnim matrikam (redke in slično).
No na decimalkah je lahko... odvisno od implementacije, od natancnosti (float/double), pa mogoce je tudi kak premik za pol tocke ali kaj takega. NR moras pa jasno znat gledat. Ze ce si filal array od 0 namesto od 1 naprej, ali pa nisi dobro gledal kako imajo realne in imaginarne komponente noter zlozene, je lahko marsikaj narobe. NR uporabljaj kot zelo dober ucbenik kako delujejo algoritmi, ampak ne programiraj z njim. Samo probleme ti naredi.Pri sinusni foruirerovi transformaciji, potegnjeni s spleta, je razlika predvsem na decimalkah. Pri NR pa predznakih in celih številih.
Ja, pred tem korakom moras se sodo/liho razsirit - zato, ker med dst/dct baznimi funkcijami so tudi tiste, ki imajo samo pol sinusa/pol kosinusa (za faktor 2 so razlicne frekvence). Tako da ce nimas specialisticne dst/dct implementacije, pac porabis 2x vec spomina ampak to niti ni hudo. V fftw pa imas vse to notri. Pa malo poglej, tukaj so tiste podrobnosti kjer je vazno kako je vse to zrcaljeno - imas DCT-I, DCT-II, DCT-III in DCT-IV.Dokler se nisem ukvarjal s fourierovo transformacijo pri numeriki, sem bil prepričan, da je tako. Ampak imaginarni del fft ni enako dst.
Re: GSL
Aja se to, ce si na C++, pa ce hoces enostavnost implementacije in te hitrost ne boli toliko - ce hoces bolj na matlabovski nacin pisat - si lahko pogledas kako ima matrike urejene knjiznica boost, ki je OGROMNA razsiritvena knjiznica za C++ (prava ropotarnica funkcij). Ta ti naceloma dovoli, da z matrikami delas kar s pomocjo operatorjev, v stilu
C=A*B;
Narejeno je s Templati, tako da izberes bazni tip matrike. Recimo matrix<double> C(n,m);
Ker je to C++, ne pa C, in ker je narejeno na genericen nacin (in s tem ni ravno optimizirano za hitrost ampak za prakticnost), lahko pricakujes, da ne bo ravno optimalna hitrost, ampak dvomim da imas take vrste problem, kjer bi bilo to ravno kriticnega pomena za vsak procent hitrosti.
http://www.boost.org/doc/libs/1_52_0/li ... matrix.htm
http://www.boost.org/doc/libs/1_54_0/li ... /index.htm
C=A*B;
Narejeno je s Templati, tako da izberes bazni tip matrike. Recimo matrix<double> C(n,m);
Ker je to C++, ne pa C, in ker je narejeno na genericen nacin (in s tem ni ravno optimizirano za hitrost ampak za prakticnost), lahko pricakujes, da ne bo ravno optimalna hitrost, ampak dvomim da imas take vrste problem, kjer bi bilo to ravno kriticnega pomena za vsak procent hitrosti.
http://www.boost.org/doc/libs/1_52_0/li ... matrix.htm
http://www.boost.org/doc/libs/1_54_0/li ... /index.htm
Re: GSL
Fajn je ekonomistom, uporabljajo Excel pa Matlab, za kaj več se jim pa ni potrebno beliti glave. Nisem si mislil, da se dejansko gre pri teh računanjih v take detajle. Ponavadi se pohitritev naredi s paralelizacijo z večimi jedri ozirom kar z uporabo kake farme grafičnih kartic.Aniviller napisal/-a:Tukaj je VELIKO prostora za optimizacijo. Ko gre enkrat za optimizacijo tako trivialnih procesov kot je mnozenje in sestevanje, se pozna vsaka neumnost, in je vitalnega pomena znanje o delovanju procesorja. Vrstni red kako gres po matriki recimo vpliva na to, kako ucinkovita je raba Cache-ja in kolikokrat procesor ustavi in caka na dostavo podatkov iz RAM-a. Ze samo to lahko pridela nekaj deset procentov, ce pazis. Vprasanje je recimo tudi ali najprej jih par zmnozis, in spravis v lokalne spremenljivke in na koncu sestejes, ali delas sproti,... Da ne govorimo o uporabi posebnih procesorskih instrukcij za namen sestevanja (imas instrukcije, ki naredijo ax+b v eni potezi - mnozenje in sestevanje), pa vzporedno racunanje (moderni intel in amd imajo vsi sse (sse2,3,...) instrukcije, ki dovoljujejo mnozenje / sestevanje po 4 pare stevil hkrati, kar sicer prevajalnik mogoce zna sam pogruntat, ampak v praksi je cloveski poseg boljsi)... potem imas pa seveda se algoritme, ki mnozijo matrike tako, da naredijo manj mnozenj kot bi bilo naivno treba: glej recimo Strassenov algoritem http://en.wikipedia.org/wiki/Strassen_algorithm ki spravi mnozenje 2x2 matrik iz 8 na 7 mnozenj (in spet dela rekurzivno, kar pomeni da je spet blazno boljse imet matrike dimenzij 2^k).Hvala, sicer sem kar sam napisal zanke za množenje matrik z vektorjem in podobno. Tu najbrž itak ni prostora za nadaljno optimizacijo.
Stvar je tako kriticna (posebej za ljudi, ki imajo opravka z gromozanskimi matrikami in programi tecejo po vec dni), da pri prevajanju knjiznic, ki implementirajo blas funkcije, zenejo poskusanje in kalibracijo na tocno dolocen racunalnik (kombinacija ram, procesor, maticna plosca - vse hitrosti in odzivni casi), poskusa prevest na precej razlicnih nacinov in izmeri case, in uporabi najhitrejsega. ATLAS recimo pocne tocno to - poskusa zagotovit, da dobis blas, ki na danem racunalniku dela kar se da hitro. Ce dobis 15% pohitritev se ze kar pozna na hitrosti, ce imas res kriticne racune.
Re: GSL
Racunat moras s tem, da so paketi za linearno algebro obstajali desetletja pred vecjedrnimi procesorji, univerzalnimi graficnimi cipi in celo pred SSE instrukcijami. Atlas, ki sicer pridela lahko priblizno faktor 2 v primerjavi z referencnim BLAS-om, ne izkorisca vecjedrnosti na ravni knjiznice (kar ni nujno slabo, saj lahko da hoces raje vec vzporednih matrik racunat in na tem nivoju izkoristit vec jeder, kar je ponavadi bolj ucinkovito). Vecjedrnost in dodatne optimizacije izkoriscajo paketi kot sta Atlas39 in GOTO blas, ki so zato se bistveno hitrejsi. Seveda je tudi tukaj treba pazit kaj delas, naivna uporaba vecjedrnosti lahko zaradi delitve pomnilnika med jedri pridela vcasih celo upocasnitev! Vse te implementacije so kompatibilne med seboj in knjiznico menjas brez sprememb v kodi. To je dobra prednost - izberes med razlicnimi implementacijami istih funkcij.
Za C++ imas sicer tudi Eigen, ki ni BLAS kompatibilen, je pa zelo hiter.
Za C++ imas sicer tudi Eigen, ki ni BLAS kompatibilen, je pa zelo hiter.