A hierarchically blocked Jacobi SVD algorithm for single and multiple graphics processing units

We present a hierarchically blocked one-sided Jacobi algorithm for the singular value decomposition (SVD), targeting both single and multiple graphics processing units (GPUs). The blocking structure reflects the levels of GPU's memory hierarchy. The algorithm may outperform MAGMA's dgesvd, while retaining high relative accuracy. To this end, we developed a family of parallel pivot strategies on GPU's shared address space, but applicable also to inter-GPU communication. Unlike common hybrid approaches, our algorithm in a single GPU setting needs a CPU for the controlling purposes only, while utilizing GPU's resources to the fullest extent permitted by the hardware. When required by the problem size, the algorithm, in principle, scales to an arbitrary number of GPU nodes. The scalability is demonstrated by more than twofold speedup for sufficiently large matrices on a Tesla S2050 system with four GPUs vs. a single Fermi card.Comment: Accepted for publication in SIAM Journal on Scientific Computin

Paralelni algoritmi Jacobijeva tipa za singularnu i generaliziranu singularnu dekompoziciju

In this thesis, a hierarchically blocked one-sided Jacobi algorithm for the singular value decomposition (SVD) is presented. The algorithm targets both single and multiple graphics processing units (GPUs). The blocking structure reflects the levels of the GPU’s memory hierarchy. To this end, a family of parallel pivot strategies on the GPU’s shared address space has been developed, but the strategies are applicable to inter-node communication as well, with GPU nodes, CPU nodes, or, in general, any NUMA nodes. Unlike common hybrid approaches, the presented algorithm in a single-GPU setting needs a CPU for the controlling purposes only, while utilizing the GPU’s resources to the fullest extent permitted by the hardware. When required by the problem size, the algorithm, in principle, scales to an arbitrary number of GPU nodes. The scalability is demonstrated by more than twofold speedup for sufficiently large matrices on a four-GPU system vs. a single GPU. The subsequent part of the thesis describes how to modify the two-sided Hari–Zimmermann algorithm for computation of the generalized eigendecomposition of a symmetric matrix pair (A; B), where B is positive definite, to an implicit algorithm that computes the generalized singular value decomposition (GSVD) of a pair (F; G). In addition, blocking and parallelization techniques for accelerating both the CPU and the GPU computation are presented, with the GPU approach following the Jacobi SVD algorithm from the first part of the thesis. For triangular matrix pairs of a moderate size, numerical tests show that the double precision sequential pointwise algorithm is several times faster than the established DTGSJA algorithm in LAPACK, while the accuracy is slightly better, especially for the small generalized singular values. Cache-aware blocking increases the performance even further. As with the one-sided Jacobi-type (G)SVD algorithms in general, the presented algorithm is almost perfectly parallelizable and scalable on the shared memory machines, where the speedup almost solely depends on the number of cores used. A distributed memory variant, intended for huge matrices that do not fit into a single NUMA node, as well as a GPU variant, are also sketched. The thesis concludes with the affirmative answer to a question whether the onesided Jacobi-type algorithms can be an efficient and scalable choice for computing the (G)SVD of dense matrices on the massively parallel CPU and GPU architectures. Unless otherwise noted by the inline citations or implied by the context, this thesis is an overview of the original research results, most of which has already been published in [55, 58]. The author’s contributions are the one-sided Jacobi-type GPU algorithms for the ordinary and the generalized SVD, of which the latter has not yet been published, as well as the parallelization technique and some implementation details of the one-sided Hari–Zimmermann CPU algorithm for the GSVD. The rest is joint work with Sanja and Saša Singer.Singularna dekompozicija, katkad zvana prema engleskom originalu i dekompozicija singularnih vrijednosti, ili kraće SVD, jedna je od najkorisnijih matričnih dekompozicija, kako za teorijske, tako i za praktične svrhe. Svaka matrica $G \in \mathbb{C}^{m \times n}$ (zbog jednostavnijeg zapisa, uobičajeno se smatra da je $m \geq n$; u protivnom, traži se SVD matrice $G^\ast$) može se rastaviti u produkt tri matrice $$G = U \Sigma V^\ast,$$ gdje su $U \in \mathbb{C}^{m \times m}$ i $V \in \mathbb{C}^{n \times n}$ unitarne, a $\Sigma \in \mathbb{R}^{m \times n}$ je 'dijagonalna' s nenegativnim dijagonalnim elementima. Osim ovog oblika dekompozicije, koristi se i skraćeni oblik $$G = U'\Sigma'V^\ast,$$ pri čemu je $U' \in \mathbb{C}^{m \times n}$ matrica s ortonormiranim stupcima, a $\Sigma' = diag(\sigma_1, \dots, \sigma_n), \sigma_i \geq 0$ za $i = 0, \dots, n$, je sada stvarno dijagonalna. Izvan matematike, u 'stvarnom' životu, SVD se koristi u procesiranju slika (rekonstrukciji, sažimanju, izoštravanju) i signala, s primjenama u medicini (CT, tj. kompjuterizirana tomografija; MR, tj. magnetna rezonancija), geoznanostima, znanosti o materijalima, kristalografiji, sigurnosti (prepoznavanje lica), izvlačenja informacija iz velike količine podataka (na primjer, LSI, tj. latent semantic indexing), ali i drugdje. Većina primjena koristi svojstvo da se iz SVD-a lako čita najbolja aproksimacija dane matrice matricom fiksnog (niskog) ranga. Čini se da je lakše reći gdje se SVD ne koristi, nego gdje se koristi, stoga se SVD često naziva i "švicarskim nožićem matričnih dekompozicija"$^1$. Prvi počeci razvoja SVD-a sežu u 19. stoljeće, kad su poznati matematičari Eugenio Beltrami, Camille Jordan, James Joseph Sylvester, Erhard Schmidt i Herman Weyl pokazali njezinu egzistenciju i osnovna svojstva (za detalje pogledati [74]). Pioniri u numeričkom računanju SVD-a su Ervand George Kogbetliantz, te Gene Golub i William Kahan, koji su razvili algoritam za računanje (bidijagonalni QR), koji je dvadeset i pet godina vladao scenom numeričkog računanja SVD-a. U to vrijeme, sveučilište Stanford (gdje je Gene Golub radio) bilo je 'glavno sjedište' za razvoj primjena SVD-a. Početkom devedesetih godina, 'sjedište SVD-a' preseljeno je u Europu, nakon objave članka [21] o relativnoj točnosti računanja svojstvenih vrijednosti simetričnih pozitivno definitnih matrica korištenjem Jacobijeve metode. Naime, problem računanja svojstvene dekompozicije pozitivno definitne matrice i problem računanja SVD-a usko su vezani. Ako je poznata dekompozicija singularnih vrijednosti matrice $G$ punog stupčanog ranga, $G \in \mathbb{C}^{m \times n} = U \Sigma V^\ast$, pri čemu je $G$ faktor matrice $A$, $A = G \ast G$, onda je $A$ simetrična i pozitivno definitna i vrijedi $$A = G \ast G = V \Sigma^T U^\ast U \Sigma V^\ast = V diag(\sigma_1^2, \dots, \sigma_m^2)V^\ast .$$ Matrica $V$ je matrica svojstvenih vektora, a svojstvene vrijednosti su kvadrati singularnih vrijednosti. Stoga se algoritmi za računanje svojstvenih vrijednosti, kod kojih se transformacija vrši dvostranim (i slijeva i zdesna) djelovanjem na matricu $A$, mogu napisati implicitno, tako da se transformacija vrši ili zdesna na faktor $G$ ili slijeva na faktor $G^\ast$. U svojoj doktorskoj disertaciji Drmač [24] je napravio daljnju analizu, ne samo singularne dekompozicije računate Jacobijevim algoritmom, nego i generalizirane singularne dekompozicije (GSVD). Temeljem tih istraživanja, SVD baziran na Jacobijevim rotacijama ušao je i u numeričku biblioteku LAPACK. U međuvremenu, gotovo sva računala postala su višejezgrena, a moderni klasteri računala za znanstveno računanje sastoje se od nekoliko tisuća do nekoliko stotina tisuća višejezgrenih procesora$^2$, pa standardni sekvencijalni algoritmi nipošto više nisu primjereni za numeričko računanje. Stoga se ubrzano razvijaju paralelni algoritmi koji poštuju i hijerarhijsku memorijsku strukturu odgovarajućih računala, težeći iskoristiti brzu cache memoriju za procesiranje potproblema u blokovima, na koje je moguće primijeniti BLAS-3 operacije. Ideja blokiranja je u primjeni što više (tipično, kubično u dimenziji matrice) numeričkih operacija nad podacima u brzoj memoriji. Nadalje, pojavom grafičkih procesnih jedinica namijenjenih znanstvenom računanju, kao i drugih visokoparalelnih numeričkih akceleratora (npr. Intel Xeon Phi), otvorio se novi segment istraživanja, koji poštuje njihov masivni paralelizam, s pojedinačno slabašnom snagom svake dretve u odnosu na središnji procesor. Generaliziranu singularnu dekompoziciju (GSVD) uveli su Van Loan [77], te Paige i Saunders [62]. Definicija GSVD-a nešto je manje poznata. Ako su zadane matrice $F \in \mathbb{C}^{m \times n}$ i $G \in \mathbb{C}^{p \times n}$, za koje vrijedi $$K = {F \brack G} , k = rank(K),$$ tad postoje unitarne matrice $U \in \mathbb{C}^{m \times m}, V \in \mathbb{C}^{p \times p}$, i matrica $X \in \mathbb{C}^{k \times n}$, takve da je $$F = U \Sigma_F X, \qquad G = V \Sigma_G X, \qquad \Sigma_F \in \mathbb{R}^{m \times k}, \qquad \Sigma_G \in \mathbb{R}^{p \times k}.$$ Elementi matrica $\Sigma_F$ i $\Sigma_G$ su nula, osim dijagonalnih elemenata, koji su realni i nenegativni. Nadalje, $\Sigma_F$ i $\Sigma_G$ zadovoljavaju $$\Sigma_F^T\Sigma_F + \Sigma_G^T\Sigma_G = I.$$ Omjeri $(\Sigma_F)_{ii} / (\Sigma_G)_{ii}$ su generalizirane singularne vrijednosti para $(F, G)$. Ako je $G$ punog stupčanog ranga, tada je $rank(K) = n$ i generalizirane singularne vrijednosti su konačni brojevi. Ako je par $(F, G)$ realan, onda su realne sve matrice u dekompoziciji. Odavde nadalje, zbog jednostavnoti pretpostavlja se da je par realan. Može se pokazati da, ako je $k = n$, tada se relacija između GSVD-a i reducirane forme CS (kosinus-sinus) dekompozicije (vidjeti, na primjer, [26]) može iskoristiti za njezino računanje (pogledati, na primjer članke Stewarta [72, 73] i Suttona [75]). Slično kao i SVD, generalizirana singularna dekompozicija ima primjene u mnogim područjima, kao što je usporedna analiza podataka vezanih uz genome [1], nepotpuna singularna metoda rubnih elemeneata [47], ionosferna tomografija [9], ali i mnogo drugih. GSVD para matrica $(F, G)$ blisko je vezana s hermitskim generaliziranim svojstvenim problemom za par $(A, B) := (F^\ast F, G^\ast G)$, tako da se metode za istovremenu dijagonalizaciju para $(A, B)$ mogu modificirati za računanje GSVD-a para $(F, G)$. U ovoj radnji razvijen je brzi i efikasan algoritam za računanje generalizirane singularne dekompozicije realnog para $(F, G)$. Metoda razvijena u radnji bazirana je na algoritmu za računanje generalizirane svojstvene dekompozicije, $$Ax = \lambda Bx; \quad x \neq 0; \qquad (1)$$ gdje su $A$ i $B$ simetrične matrice, a par je definitan, tj. postoji realna konstanta $\mu$ takva da je matrica $A-\mu B$ pozitivno definitna. Članke s metodom objavili su 1960. Falk i Langemeyer [31, 32] u slabo poznatom priručniku. Kad je paralelna verzija metode testirana, pokazalo se da pati zbog problema rastuće skale stupaca matrice tijekom procesa ortogonalizacije. Treba još primijetiti da pozitivna definitnost matrice $B$ odmah znači da je definitan i par $(A, B)$. Gotovo desetljeće nakon Falka i Langemeyera, Katharina Zimmermann je u svojoj doktorskoj disertaciji [81] grubo skicirala metodu za rješavanje generaliziranog svojstvenog problema (1) ako je B pozitivno definitna. Gose [34] je predložio optimalnu ne-cikličku pivotnu strategiju i dokazao globalnu konvergenciju originalne metode. Hari je u svojoj disertaciji [37], potaknut Zimmermanninom skicom metode, izveo algoritam i pokazao njegovu globalnu i kvadratičnu konvergenciju uz cikličke pivotne strategije. Kvadratičnu konvergenciju originalne Falk–Langemeyerove metode dokazao je 1988. Slapničar u svojem magisteriju, četiri godine nakon dokaza konvergencije Hari–Zimmermann metode. Hari je u [37] pokazao ključnu vezu između Hari–Zimmermannine i Falk–Langemeyerove varijante algoritma. Ako je matrica $B$ obostrano skalirana dijagonalnom matricom $D$, tako da su joj dijagonalni elementi jednaki 1 prije svakog koraka poništavanja u Falk–Langemeyerovoj metodi, dobiva se Hari–Zimmermannina metoda. Dakle, nova metoda imala je ključno svojstvo normiranosti stupaca barem jedne matrice, što se pokazalo iznimno bitnim za uspjeh algoritma (izbjegavanje skaliranja matrica tijekom procesa ortogonalizacije). Treba reći da se GSVD može računati i na druge načine. Drmač je u [26] izveo algoritam za računanje GSVD-a para $(F, G)$, kad je $G$ punog stupčanog ranga. Algoritam transformira problem na samo jednu matricu, a nakon toga primjenjuje jednostrani Jacobijev SVD algoritam. Taj algoritam računa generalizirane singularne vrijednosti s malom relativnom greškom. Algoritam svođenja na jednu matricu sastoji se od tri koraka: skaliranje stupaca matrica $F$ i $G$, QR faktorizacije sa stupčanim pivotiranjem već skalirane matrice $G$, i konačno, rješavanjem trokutastog linearnog sustava s $k$ desnih strana. Posljednja dva koraka su sekvencijalna i vrlo ih je teško paralelizirati. Sama ideja korištenja implicitne (tj. jednostrane) Falk–Langemeyerove metode za GSVD para $(F, G)$, s $G$ punog stupčanog ranga, sreće se u disertaciji Annette Deichmöller [17], međutim, tamo se ne spominju usporedbe te metode s drugim metodama. S druge strane, algoritam za računanje GSVD-a u biblioteci LAPACK (potprogram xGGSVD), je modificirani Kogbetliantzov algoritam (vidjeti Paige [61]) s obveznim pretprocesiranjem (vidjeti Bai i Demmel [5]). Algoritam pretprocesiranja [6] transformira zadani matrični par $(F_0, G_0)$ u par $(F, G)$, takav da su $F$ i $G$ gornjetrokutaste, a $G$ je i nesingularna. Ako se unaprijed zna da je $G$ punog stupčanog ranga, i implicitna Falk–Langemeyerova i implicitna Hari–Zimmermannina metoda će raditi i bez pretprocesiranja. Ako su $F$ i $G$ vitke (engl. "tall and skinny"), QR factorizacija obje matrice će ubrzati ortogonalizaciju. Ako $G$ nije punog ranga, onda treba koristiti isto pretprocesiranje kao u LAPACK-u, budući da puni stupčani rang matrice $G$ garantira pozitivnu definitnost matrice $B := G^T G$. U ovoj radnji razvijen je i hijerarhijski, blokirani jednostrani algoritam za računanje SVD-a. Opisani algoritam može raditi na višeprocesorskom računalu, računalnim klasterima, jednoj ili više grafičkih procesnih jedinica. Princip rada algoritma na svim arhitekturama je sličan. Posebno je opisan algoritam koji radi na grafičkim procesnim jedinicama. Struktura blokiranja reflektira razine memorijske strukture grafičke procesne jedninice. Da bi se to postiglo, razvijene su familije paralelnih pivotnih strategija za dijeljenu (engl. shared) memoriju grafičkih procesnih jedinica. Uz dodatak rasporeda po procesima, strategije se mogu koristiti i kao strategije za komuniciranje među računalnim čvorovima (bili oni grafičke procesne jedinice, jezgre procesora ili tzv. NUMA čvorovi). Razvijeni algoritam nije hibridni, tj. centralnu procesnu jedinicu koristi samo za kontrolne svrhe, a cjelokupno računanje odvija se na grafičkoj procesnoj jedinici. Kad je zbog veličine problema potrebno, algoritam se može rasprostrijeti (skalirati) na proizvoljan broj grafičkih procesnih jedinica. Na dovoljno velikim matricama, skalabilnost je pokazana ubrzanjem od preko dva puta na četiri grafičke procesne jedinice, obzirom na jednu. U drugom dijelu radnje opisuje se jedan način modifikacije dvostranog Hari–Zimmermanninog algoritma za računanje generalizirane svojstvene dekompozicije matričnog para $(A, B)$, gdje su obje matrice simetrične, a $B$ je pozitivno definitna. Implicitni algoritam računa GSVD para $(F, G)$, pri čemu je $(A, B) := (F^T F, G^T G)$. Nadalje, pokazuje se kako treba blokirati algoritam, te kako ga paralelizirati, i u slučaju standardnih, i u slučaju grafičkih procesora. Za trokutaste matrične parove srednje velikih dimenzija (približno 5 000), pokazano je da je već sekvencijalni, neblokirani algoritam u dvostrukoj točnosti, predložen u radnji, nekoliko desetaka puta brži no što je to LAPACK potprogram DTGSJA i pritom ima nešto bolju točnost, posebno za male generalizirane singularne vrijednosti. Blokiranje algoritma koje odgovara cacheima znatno ubrzava algoritam. Pokazuje se da je i ovaj algoritam, slično kao jednostrani Jacobijev algoritam za SVD, gotovo idealno paralelizabilan i skalabilan na računalima s dijeljenom memorijom, te da njegovo ubrzanje gotovo isključivo ovisi o broju korištenih jezgara. U vrijeme testiranja, pokazalo se da je paralelizirani i blokirani Hari–Zimmermannin algoritam preko sto puta brži od LAPACK potprograma DTGESJA s višedretvenim BLAS potprogramima. Varijanta algoritma za razdijeljenu (engl. distributed) memoriju namijenjena je ogromnim matricama koje ne stanu u jedan NUMA čvor. Također, skicirana je i GPU varijanta algoritma, koja je vrlo slična jednostranom Jacobijevom algoritmu za SVD. Disertacija završava zaključkom da su ovi algoritmi Jacobijevog tipa efikasni i skalabilni i izvrstan su izbor za računanje (G)SVD-a punih matrica na masivno paralelnim standardnim arhitekturama i na grafičkim procesnim jedinicama. Ova doktorska disertacija bazirana je na originalnim znanstvenim radovima [55, 58], te proširena nekim novim rezultatima. Autorov doprinos u ovoj disertaciji su novi paralelni algoritmi za (G)SVD za grafičke procesne jedinice, tehnike paralelizacije, te detalji implementacije jednostranog Hari–Zimmermannina algoritma. Ostatak je zajednički rad sa Sanjom Singer i Sašom Singerom. $^1$ Diane O’Leary, 2006. $^2$ https://www.top500.or

Paralelni algoritmi Jacobijeva tipa za singularnu i generaliziranu singularnu dekompoziciju

In this thesis, a hierarchically blocked one-sided Jacobi algorithm for the singular value decomposition (SVD) is presented. The algorithm targets both single and multiple graphics processing units (GPUs). The blocking structure reflects the levels of the GPU’s memory hierarchy. To this end, a family of parallel pivot strategies on the GPU’s shared address space has been developed, but the strategies are applicable to inter-node communication as well, with GPU nodes, CPU nodes, or, in general, any NUMA nodes. Unlike common hybrid approaches, the presented algorithm in a single-GPU setting needs a CPU for the controlling purposes only, while utilizing the GPU’s resources to the fullest extent permitted by the hardware. When required by the problem size, the algorithm, in principle, scales to an arbitrary number of GPU nodes. The scalability is demonstrated by more than twofold speedup for sufficiently large matrices on a four-GPU system vs. a single GPU. The subsequent part of the thesis describes how to modify the two-sided Hari–Zimmermann algorithm for computation of the generalized eigendecomposition of a symmetric matrix pair (A; B), where B is positive definite, to an implicit algorithm that computes the generalized singular value decomposition (GSVD) of a pair (F; G). In addition, blocking and parallelization techniques for accelerating both the CPU and the GPU computation are presented, with the GPU approach following the Jacobi SVD algorithm from the first part of the thesis. For triangular matrix pairs of a moderate size, numerical tests show that the double precision sequential pointwise algorithm is several times faster than the established DTGSJA algorithm in LAPACK, while the accuracy is slightly better, especially for the small generalized singular values. Cache-aware blocking increases the performance even further. As with the one-sided Jacobi-type (G)SVD algorithms in general, the presented algorithm is almost perfectly parallelizable and scalable on the shared memory machines, where the speedup almost solely depends on the number of cores used. A distributed memory variant, intended for huge matrices that do not fit into a single NUMA node, as well as a GPU variant, are also sketched. The thesis concludes with the affirmative answer to a question whether the onesided Jacobi-type algorithms can be an efficient and scalable choice for computing the (G)SVD of dense matrices on the massively parallel CPU and GPU architectures. Unless otherwise noted by the inline citations or implied by the context, this thesis is an overview of the original research results, most of which has already been published in [55, 58]. The author’s contributions are the one-sided Jacobi-type GPU algorithms for the ordinary and the generalized SVD, of which the latter has not yet been published, as well as the parallelization technique and some implementation details of the one-sided Hari–Zimmermann CPU algorithm for the GSVD. The rest is joint work with Sanja and Saša Singer.Singularna dekompozicija, katkad zvana prema engleskom originalu i dekompozicija singularnih vrijednosti, ili kraće SVD, jedna je od najkorisnijih matričnih dekompozicija, kako za teorijske, tako i za praktične svrhe. Svaka matrica $G \in \mathbb{C}^{m \times n}$ (zbog jednostavnijeg zapisa, uobičajeno se smatra da je $m \geq n$; u protivnom, traži se SVD matrice $G^\ast$) može se rastaviti u produkt tri matrice $$G = U \Sigma V^\ast,$$ gdje su $U \in \mathbb{C}^{m \times m}$ i $V \in \mathbb{C}^{n \times n}$ unitarne, a $\Sigma \in \mathbb{R}^{m \times n}$ je 'dijagonalna' s nenegativnim dijagonalnim elementima. Osim ovog oblika dekompozicije, koristi se i skraćeni oblik $$G = U'\Sigma'V^\ast,$$ pri čemu je $U' \in \mathbb{C}^{m \times n}$ matrica s ortonormiranim stupcima, a $\Sigma' = diag(\sigma_1, \dots, \sigma_n), \sigma_i \geq 0$ za $i = 0, \dots, n$, je sada stvarno dijagonalna. Izvan matematike, u 'stvarnom' životu, SVD se koristi u procesiranju slika (rekonstrukciji, sažimanju, izoštravanju) i signala, s primjenama u medicini (CT, tj. kompjuterizirana tomografija; MR, tj. magnetna rezonancija), geoznanostima, znanosti o materijalima, kristalografiji, sigurnosti (prepoznavanje lica), izvlačenja informacija iz velike količine podataka (na primjer, LSI, tj. latent semantic indexing), ali i drugdje. Većina primjena koristi svojstvo da se iz SVD-a lako čita najbolja aproksimacija dane matrice matricom fiksnog (niskog) ranga. Čini se da je lakše reći gdje se SVD ne koristi, nego gdje se koristi, stoga se SVD često naziva i "švicarskim nožićem matričnih dekompozicija"$^1$. Prvi počeci razvoja SVD-a sežu u 19. stoljeće, kad su poznati matematičari Eugenio Beltrami, Camille Jordan, James Joseph Sylvester, Erhard Schmidt i Herman Weyl pokazali njezinu egzistenciju i osnovna svojstva (za detalje pogledati [74]). Pioniri u numeričkom računanju SVD-a su Ervand George Kogbetliantz, te Gene Golub i William Kahan, koji su razvili algoritam za računanje (bidijagonalni QR), koji je dvadeset i pet godina vladao scenom numeričkog računanja SVD-a. U to vrijeme, sveučilište Stanford (gdje je Gene Golub radio) bilo je 'glavno sjedište' za razvoj primjena SVD-a. Početkom devedesetih godina, 'sjedište SVD-a' preseljeno je u Europu, nakon objave članka [21] o relativnoj točnosti računanja svojstvenih vrijednosti simetričnih pozitivno definitnih matrica korištenjem Jacobijeve metode. Naime, problem računanja svojstvene dekompozicije pozitivno definitne matrice i problem računanja SVD-a usko su vezani. Ako je poznata dekompozicija singularnih vrijednosti matrice $G$ punog stupčanog ranga, $G \in \mathbb{C}^{m \times n} = U \Sigma V^\ast$, pri čemu je $G$ faktor matrice $A$, $A = G \ast G$, onda je $A$ simetrična i pozitivno definitna i vrijedi $$A = G \ast G = V \Sigma^T U^\ast U \Sigma V^\ast = V diag(\sigma_1^2, \dots, \sigma_m^2)V^\ast .$$ Matrica $V$ je matrica svojstvenih vektora, a svojstvene vrijednosti su kvadrati singularnih vrijednosti. Stoga se algoritmi za računanje svojstvenih vrijednosti, kod kojih se transformacija vrši dvostranim (i slijeva i zdesna) djelovanjem na matricu $A$, mogu napisati implicitno, tako da se transformacija vrši ili zdesna na faktor $G$ ili slijeva na faktor $G^\ast$. U svojoj doktorskoj disertaciji Drmač [24] je napravio daljnju analizu, ne samo singularne dekompozicije računate Jacobijevim algoritmom, nego i generalizirane singularne dekompozicije (GSVD). Temeljem tih istraživanja, SVD baziran na Jacobijevim rotacijama ušao je i u numeričku biblioteku LAPACK. U međuvremenu, gotovo sva računala postala su višejezgrena, a moderni klasteri računala za znanstveno računanje sastoje se od nekoliko tisuća do nekoliko stotina tisuća višejezgrenih procesora$^2$, pa standardni sekvencijalni algoritmi nipošto više nisu primjereni za numeričko računanje. Stoga se ubrzano razvijaju paralelni algoritmi koji poštuju i hijerarhijsku memorijsku strukturu odgovarajućih računala, težeći iskoristiti brzu cache memoriju za procesiranje potproblema u blokovima, na koje je moguće primijeniti BLAS-3 operacije. Ideja blokiranja je u primjeni što više (tipično, kubično u dimenziji matrice) numeričkih operacija nad podacima u brzoj memoriji. Nadalje, pojavom grafičkih procesnih jedinica namijenjenih znanstvenom računanju, kao i drugih visokoparalelnih numeričkih akceleratora (npr. Intel Xeon Phi), otvorio se novi segment istraživanja, koji poštuje njihov masivni paralelizam, s pojedinačno slabašnom snagom svake dretve u odnosu na središnji procesor. Generaliziranu singularnu dekompoziciju (GSVD) uveli su Van Loan [77], te Paige i Saunders [62]. Definicija GSVD-a nešto je manje poznata. Ako su zadane matrice $F \in \mathbb{C}^{m \times n}$ i $G \in \mathbb{C}^{p \times n}$, za koje vrijedi $$K = {F \brack G} , k = rank(K),$$ tad postoje unitarne matrice $U \in \mathbb{C}^{m \times m}, V \in \mathbb{C}^{p \times p}$, i matrica $X \in \mathbb{C}^{k \times n}$, takve da je $$F = U \Sigma_F X, \qquad G = V \Sigma_G X, \qquad \Sigma_F \in \mathbb{R}^{m \times k}, \qquad \Sigma_G \in \mathbb{R}^{p \times k}.$$ Elementi matrica $\Sigma_F$ i $\Sigma_G$ su nula, osim dijagonalnih elemenata, koji su realni i nenegativni. Nadalje, $\Sigma_F$ i $\Sigma_G$ zadovoljavaju $$\Sigma_F^T\Sigma_F + \Sigma_G^T\Sigma_G = I.$$ Omjeri $(\Sigma_F)_{ii} / (\Sigma_G)_{ii}$ su generalizirane singularne vrijednosti para $(F, G)$. Ako je $G$ punog stupčanog ranga, tada je $rank(K) = n$ i generalizirane singularne vrijednosti su konačni brojevi. Ako je par $(F, G)$ realan, onda su realne sve matrice u dekompoziciji. Odavde nadalje, zbog jednostavnoti pretpostavlja se da je par realan. Može se pokazati da, ako je $k = n$, tada se relacija između GSVD-a i reducirane forme CS (kosinus-sinus) dekompozicije (vidjeti, na primjer, [26]) može iskoristiti za njezino računanje (pogledati, na primjer članke Stewarta [72, 73] i Suttona [75]). Slično kao i SVD, generalizirana singularna dekompozicija ima primjene u mnogim područjima, kao što je usporedna analiza podataka vezanih uz genome [1], nepotpuna singularna metoda rubnih elemeneata [47], ionosferna tomografija [9], ali i mnogo drugih. GSVD para matrica $(F, G)$ blisko je vezana s hermitskim generaliziranim svojstvenim problemom za par $(A, B) := (F^\ast F, G^\ast G)$, tako da se metode za istovremenu dijagonalizaciju para $(A, B)$ mogu modificirati za računanje GSVD-a para $(F, G)$. U ovoj radnji razvijen je brzi i efikasan algoritam za računanje generalizirane singularne dekompozicije realnog para $(F, G)$. Metoda razvijena u radnji bazirana je na algoritmu za računanje generalizirane svojstvene dekompozicije, $$Ax = \lambda Bx; \quad x \neq 0; \qquad (1)$$ gdje su $A$ i $B$ simetrične matrice, a par je definitan, tj. postoji realna konstanta $\mu$ takva da je matrica $A-\mu B$ pozitivno definitna. Članke s metodom objavili su 1960. Falk i Langemeyer [31, 32] u slabo poznatom priručniku. Kad je paralelna verzija metode testirana, pokazalo se da pati zbog problema rastuće skale stupaca matrice tijekom procesa ortogonalizacije. Treba još primijetiti da pozitivna definitnost matrice $B$ odmah znači da je definitan i par $(A, B)$. Gotovo desetljeće nakon Falka i Langemeyera, Katharina Zimmermann je u svojoj doktorskoj disertaciji [81] grubo skicirala metodu za rješavanje generaliziranog svojstvenog problema (1) ako je B pozitivno definitna. Gose [34] je predložio optimalnu ne-cikličku pivotnu strategiju i dokazao globalnu konvergenciju originalne metode. Hari je u svojoj disertaciji [37], potaknut Zimmermanninom skicom metode, izveo algoritam i pokazao njegovu globalnu i kvadratičnu konvergenciju uz cikličke pivotne strategije. Kvadratičnu konvergenciju originalne Falk–Langemeyerove metode dokazao je 1988. Slapničar u svojem magisteriju, četiri godine nakon dokaza konvergencije Hari–Zimmermann metode. Hari je u [37] pokazao ključnu vezu između Hari–Zimmermannine i Falk–Langemeyerove varijante algoritma. Ako je matrica $B$ obostrano skalirana dijagonalnom matricom $D$, tako da su joj dijagonalni elementi jednaki 1 prije svakog koraka poništavanja u Falk–Langemeyerovoj metodi, dobiva se Hari–Zimmermannina metoda. Dakle, nova metoda imala je ključno svojstvo normiranosti stupaca barem jedne matrice, što se pokazalo iznimno bitnim za uspjeh algoritma (izbjegavanje skaliranja matrica tijekom procesa ortogonalizacije). Treba reći da se GSVD može računati i na druge načine. Drmač je u [26] izveo algoritam za računanje GSVD-a para $(F, G)$, kad je $G$ punog stupčanog ranga. Algoritam transformira problem na samo jednu matricu, a nakon toga primjenjuje jednostrani Jacobijev SVD algoritam. Taj algoritam računa generalizirane singularne vrijednosti s malom relativnom greškom. Algoritam svođenja na jednu matricu sastoji se od tri koraka: skaliranje stupaca matrica $F$ i $G$, QR faktorizacije sa stupčanim pivotiranjem već skalirane matrice $G$, i konačno, rješavanjem trokutastog linearnog sustava s $k$ desnih strana. Posljednja dva koraka su sekvencijalna i vrlo ih je teško paralelizirati. Sama ideja korištenja implicitne (tj. jednostrane) Falk–Langemeyerove metode za GSVD para $(F, G)$, s $G$ punog stupčanog ranga, sreće se u disertaciji Annette Deichmöller [17], međutim, tamo se ne spominju usporedbe te metode s drugim metodama. S druge strane, algoritam za računanje GSVD-a u biblioteci LAPACK (potprogram xGGSVD), je modificirani Kogbetliantzov algoritam (vidjeti Paige [61]) s obveznim pretprocesiranjem (vidjeti Bai i Demmel [5]). Algoritam pretprocesiranja [6] transformira zadani matrični par $(F_0, G_0)$ u par $(F, G)$, takav da su $F$ i $G$ gornjetrokutaste, a $G$ je i nesingularna. Ako se unaprijed zna da je $G$ punog stupčanog ranga, i implicitna Falk–Langemeyerova i implicitna Hari–Zimmermannina metoda će raditi i bez pretprocesiranja. Ako su $F$ i $G$ vitke (engl. "tall and skinny"), QR factorizacija obje matrice će ubrzati ortogonalizaciju. Ako $G$ nije punog ranga, onda treba koristiti isto pretprocesiranje kao u LAPACK-u, budući da puni stupčani rang matrice $G$ garantira pozitivnu definitnost matrice $B := G^T G$. U ovoj radnji razvijen je i hijerarhijski, blokirani jednostrani algoritam za računanje SVD-a. Opisani algoritam može raditi na višeprocesorskom računalu, računalnim klasterima, jednoj ili više grafičkih procesnih jedinica. Princip rada algoritma na svim arhitekturama je sličan. Posebno je opisan algoritam koji radi na grafičkim procesnim jedinicama. Struktura blokiranja reflektira razine memorijske strukture grafičke procesne jedninice. Da bi se to postiglo, razvijene su familije paralelnih pivotnih strategija za dijeljenu (engl. shared) memoriju grafičkih procesnih jedinica. Uz dodatak rasporeda po procesima, strategije se mogu koristiti i kao strategije za komuniciranje među računalnim čvorovima (bili oni grafičke procesne jedinice, jezgre procesora ili tzv. NUMA čvorovi). Razvijeni algoritam nije hibridni, tj. centralnu procesnu jedinicu koristi samo za kontrolne svrhe, a cjelokupno računanje odvija se na grafičkoj procesnoj jedinici. Kad je zbog veličine problema potrebno, algoritam se može rasprostrijeti (skalirati) na proizvoljan broj grafičkih procesnih jedinica. Na dovoljno velikim matricama, skalabilnost je pokazana ubrzanjem od preko dva puta na četiri grafičke procesne jedinice, obzirom na jednu. U drugom dijelu radnje opisuje se jedan način modifikacije dvostranog Hari–Zimmermanninog algoritma za računanje generalizirane svojstvene dekompozicije matričnog para $(A, B)$, gdje su obje matrice simetrične, a $B$ je pozitivno definitna. Implicitni algoritam računa GSVD para $(F, G)$, pri čemu je $(A, B) := (F^T F, G^T G)$. Nadalje, pokazuje se kako treba blokirati algoritam, te kako ga paralelizirati, i u slučaju standardnih, i u slučaju grafičkih procesora. Za trokutaste matrične parove srednje velikih dimenzija (približno 5 000), pokazano je da je već sekvencijalni, neblokirani algoritam u dvostrukoj točnosti, predložen u radnji, nekoliko desetaka puta brži no što je to LAPACK potprogram DTGSJA i pritom ima nešto bolju točnost, posebno za male generalizirane singularne vrijednosti. Blokiranje algoritma koje odgovara cacheima znatno ubrzava algoritam. Pokazuje se da je i ovaj algoritam, slično kao jednostrani Jacobijev algoritam za SVD, gotovo idealno paralelizabilan i skalabilan na računalima s dijeljenom memorijom, te da njegovo ubrzanje gotovo isključivo ovisi o broju korištenih jezgara. U vrijeme testiranja, pokazalo se da je paralelizirani i blokirani Hari–Zimmermannin algoritam preko sto puta brži od LAPACK potprograma DTGESJA s višedretvenim BLAS potprogramima. Varijanta algoritma za razdijeljenu (engl. distributed) memoriju namijenjena je ogromnim matricama koje ne stanu u jedan NUMA čvor. Također, skicirana je i GPU varijanta algoritma, koja je vrlo slična jednostranom Jacobijevom algoritmu za SVD. Disertacija završava zaključkom da su ovi algoritmi Jacobijevog tipa efikasni i skalabilni i izvrstan su izbor za računanje (G)SVD-a punih matrica na masivno paralelnim standardnim arhitekturama i na grafičkim procesnim jedinicama. Ova doktorska disertacija bazirana je na originalnim znanstvenim radovima [55, 58], te proširena nekim novim rezultatima. Autorov doprinos u ovoj disertaciji su novi paralelni algoritmi za (G)SVD za grafičke procesne jedinice, tehnike paralelizacije, te detalji implementacije jednostranog Hari–Zimmermannina algoritma. Ostatak je zajednički rad sa Sanjom Singer i Sašom Singerom. $^1$ Diane O’Leary, 2006. $^2$ https://www.top500.or

Real Time Log Length Measurement Using GPU Accelerated Visual Odometry

This thesis studies GPU accelerated visual odometry in measuring log length. The visual odometry would not suffer slippage nor require recalibration depending type of wood or temperature conditions compared to mechanical measurement. The requirement of the real-time performance is quite high. Image capturing in 120 Hz frequency is needed as log is moved several meters per second by harvester heads. Here GPU acceleration will be used as it can give speedup in magnitude of hundreds or more. Real-time performance is targeted by selecting fast algorithms for subtasks of measurement pipeline and considering possibilities to parallelize algorithm. In many cases performance boost is achieved, but not in expected magnitude. Physical constraints of the graphics card hardware become easily the limiting factor in parallelization. Real-time performance was achieved in this thesis but not with required accuracy. It remained for future work to find out which algorithms would give both targets.  Taman lisensiaatintutkimuksen aiheena on GPU laskennan kaytto konenäköön perustuvassa tukin pituuden mittauksessa. Konenäköön perustuva pituuden mittaus ei tarvitse uudelleen kalibrointia puulajin tai lämpötilan mukaan. Konenäköön perustuvassa mittauksessa myöskaan mittapyöra ei voi luistaa tukin pinnalla. Realiaikaisuuden vaatimus on tässä sovelluksessa korkea. Kuvat on otettu 120 Hz taajuudella, koska leikkuupää liikuttaa tukkia useita metrejä sekunnissa. GPU laskenta potentiaalisesti nopeuttaisi laskentaa tarvittavissa määrin. Realiaikaista vastetta haettiin seka algoritmien valinnalla etta harkitsemalla mahdollisuuksia rinnaisohjelmoinnin käyttämiseen. Monessa tapauksessa vasteet paranivat, vaikka grafiikkakortin ominaisuudet usein rajoittivat rinnaikkaisohjelmoinnista saatavaa hyötyä. Realiaikainen vaste saavutettiin, mutta ei tarvittavalla pituuden mittaamisen tarkkuudella. Molempien tavoitteiden saavuttaminen jai mahdollisten jatkotöiden tehtäväksi

병렬화 용이한 통계계산 방법론과 현대 고성능 컴퓨팅 환경에의 적용

학위논문 (박사) -- 서울대학교 대학원 : 자연과학대학 통계학과, 2020. 8. 원중호.Technological advances in the past decade, hardware and software alike, have made access to high-performance computing (HPC) easier than ever. In this dissertation, easily-parallelizable, inversion-free, and variable-separated algorithms and their implementation in statistical computing are discussed. The first part considers statistical estimation problems under structured sparsity posed as minimization of a sum of two or three convex functions, one of which is a composition of non-smooth and linear functions. Examples include graph-guided sparse fused lasso and overlapping group lasso. Two classes of inversion-free primal-dual algorithms are considered and unified from a perspective of monotone operator theory. From this unification, a continuum of preconditioned forward-backward operator splitting algorithms amenable to parallel and distributed computing is proposed. The unification is further exploited to introduce a continuum of accelerated algorithms on which the theoretically optimal asymptotic rate of convergence is obtained. For the second part, easy-to-use distributed matrix data structures in PyTorch and Julia are presented. They enable users to write code once and run it anywhere from a laptop to a workstation with multiple graphics processing units (GPUs) or a supercomputer in a cloud. With these data structures, various parallelizable statistical applications, including nonnegative matrix factorization, positron emission tomography, multidimensional scaling, and ℓ1-regularized Cox regression, are demonstrated. The examples scale up to an 8-GPU workstation and a 720-CPU-core cluster in a cloud. As a case in point, the onset of type-2 diabetes from the UK Biobank with 400,000 subjects and about 500,000 single nucleotide polymorphisms is analyzed using the HPC ℓ1-regularized Cox regression. Fitting a half-million variate model took about 50 minutes, reconfirming known associations. To my knowledge, the feasibility of a joint genome-wide association analysis of survival outcomes at this scale is first demonstrated.지난 10년간의 하드웨어와 소프트웨어의 기술적인 발전은 고성능 컴퓨팅의 접근장벽을 그 어느 때보다 낮추었다. 이 학위논문에서는 병렬화 용이하고 역행렬 연산이 없는 변수 분리 알고리즘과 그 통계계산에서의 구현을 논의한다. 첫 부분은 볼록 함수 두 개 또는 세 개의 합으로 나타나는 구조화된 희소 통계 추정 문제에 대해 다룬다. 이 때 함수들 중 하나는 비평활 함수와 선형 함수의 합성으로 나타난다. 그 예시로는 그래프 구조를 통해 유도되는 희소 융합 Lasso 문제와 한 변수가 여러 그룹에 속할 수 있는 그룹 Lasso 문제가 있다. 이를 풀기 위해 역행렬 연산이 없는 두 종류의 원시-쌍대 (primal-dual) 알고리즘을 단조 연산자 이론 관점에서 통합하며 이를 통해 병렬화 용이한 precondition된 전방-후방 연산자 분할 알고리즘의 집합을 제안한다. 이 통합은 점근적으로 최적 수렴률을 갖는 가속 알고리즘의 집합을 구성하는 데 활용된다. 두 번째 부분에서는 PyTorch와 Julia를 통해 사용하기 쉬운 분산 행렬 자료 구조를 제시한다. 이 구조는 사용자들이 코드를 한 번 작성하면 이것을 노트북 한 대에서부터 여러 대의 그래픽 처리 장치 (GPU)를 가진 워크스테이션, 또는 클라우드 상에 있는 슈퍼컴퓨터까지 다양한 스케일에서 실행할 수 있게 해 준다. 아울러, 이 자료 구조를 비음 행렬 분해, 양전자 단층 촬영, 다차원 척 도법, ℓ1-벌점화 Cox 회귀 분석 등 다양한 병렬화 가능한 통계적 문제에 적용한다. 이 예시들은 8대의 GPU가 있는 워크스테이션과 720개의 코어가 있는 클라우드 상의 가상 클러스터에서 확장 가능했다. 한 사례로 400,000명의 대상과 500,000개의 단일 염기 다형성 정보가 있는 UK Biobank 자료에서의 제2형 당뇨병 (T2D) 발병 나이를 ℓ1-벌점화 Cox 회귀 모형을 통해 분석했다. 500,000개의 변수가 있는 모형을 적합시키는 데 50분 가량의 시간이 걸렸으며 알려진 T2D 관련 다형성들을 재확인할 수 있었다. 이러한 규모의 전유전체 결합 생존 분석은 최초로 시도된 것이다.Chapter1Prologue 1 1.1 Introduction 1 1.2 Accessible High-Performance Computing Systems 4 1.2.1 Preliminaries 4 1.2.2 Multiple CPU nodes: clusters, supercomputers, and clouds 7 1.2.3 Multi-GPU node 9 1.3 Highly Parallelizable Algorithms 12 1.3.1 MM algorithms 12 1.3.2 Proximal gradient descent 14 1.3.3 Proximal distance algorithm 16 1.3.4 Primal-dual methods 17 Chapter 2 Easily Parallelizable and Distributable Class of Algorithms for Structured Sparsity, with Optimal Acceleration 20 2.1 Introduction 20 2.2 Unification of Algorithms LV and CV (g ≡ 0) 30 2.2.1 Relation between Algorithms LV and CV 30 2.2.2 Unified algorithm class 34 2.2.3 Convergence analysis 35 2.3 Optimal acceleration 39 2.3.1 Algorithms 40 2.3.2 Convergence analysis 41 2.4 Stochastic optimal acceleration 45 2.4.1 Algorithm 45 2.4.2 Convergence analysis 47 2.5 Numerical experiments 50 2.5.1 Model problems 50 2.5.2 Convergence behavior 52 2.5.3 Scalability 62 2.6 Discussion 63 Chapter 3 Towards Unified Programming for High-Performance Statistical Computing Environments 66 3.1 Introduction 66 3.2 Related Software 69 3.2.1 Message-passing interface and distributed array interfaces 69 3.2.2 Unified array interfaces for CPU and GPU 69 3.3 Easy-to-use Software Libraries for HPC 70 3.3.1 Deep learning libraries and HPC 70 3.3.2 Case study: PyTorch versus TensorFlow 73 3.3.3 A brief introduction to PyTorch 76 3.3.4 A brief introduction to Julia 80 3.3.5 Methods and multiple dispatch 80 3.3.6 Multidimensional arrays 82 3.3.7 Matrix multiplication 83 3.3.8 Dot syntax for vectorization 86 3.4 Distributed matrix data structure 87 3.4.1 Distributed matrices in PyTorch: distmat 87 3.4.2 Distributed arrays in Julia: MPIArray 90 3.5 Examples 98 3.5.1 Nonnegative matrix factorization 100 3.5.2 Positron emission tomography 109 3.5.3 Multidimensional scaling 113 3.5.4 L1-regularized Cox regression 117 3.5.5 Genome-wide survival analysis of the UK Biobank dataset 121 3.6 Discussion 126 Chapter 4 Conclusion 131 Appendix A Monotone Operator Theory 134 Appendix B Proofs for Chapter II 139 B.1 Preconditioned forward-backward splitting 139 B.2 Optimal acceleration 147 B.3 Optimal stochastic acceleration 158 Appendix C AWS EC2 and ParallelCluster 168 C.1 Overview 168 C.2 Glossary 169 C.3 Prerequisites 172 C.4 Installation 173 C.5 Configuration 173 C.6 Creating, accessing, and destroying the cluster 178 C.7 Installation of libraries 178 C.8 Running a job 179 C.9 Miscellaneous 180 Appendix D Code for memory-efficient L1-regularized Cox proportional hazards model 182 Appendix E Details of SNPs selected in L1-regularized Cox regression 184 Bibliography 188 국문초록 212Docto

Parallel algorithms for computational fluid dynamics on unstructured meshes

La simulació numèrica directa (DNS) de fluxos complexes és actualment una utopia per la majoria d'aplicacions industrials ja que els requeriments computacionals son massa elevats. Donat un flux, la diferència entre els recursos computacionals necessaris i els disponibles és cobreix mitjançant la modelització/simplificació d'alguns termes de les equacions originals que regeixen el seu comportament. El creixement continuat dels recursos computacionals disponibles, principalment en forma de super-ordinadors, contribueix a reduir la part del flux que és necessari aproximar. De totes maneres, obtenir la eficiència esperada dels nous super-ordinadors no és una tasca senzilla i, per aquest motiu, part de la recerca en el camp de la Mecànica de Fluids Computacional es centra en aquest objectiu. En aquest sentit, algunes contribucions s'han presentat en el marc d'aquesta tesis. El primer objectiu va ser el desenvolupament d'un codi de CFD de propòsit general i paral·lel, basat en la metodologia de volums finits en malles no estructurades, per resoldre problemes de multi-física. Aquest codi, anomenat TermoFluids (TF), té un disseny orientat a objectes i pensat per ser usat de forma altament eficient en els super-ordinadors actuals. Amb el temps, ha esdevingut pel grup una eina fonamental en projectes tant de recerca bàsica com d'interès industrial. En el context d'aquesta tesis, el treball s'ha focalitzat en el desenvolupament de dos de les llibreries més bàsiques de TermoFluids: i) La Basics Objects Library (BOL), que es una plataforma de software sobre la qual estan programades la resta de llibreries del codi, i que conté els mètodes algebraics i geomètrics fonamentals per la implementació paral·lela dels algoritmes de discretització, ii) la Linear Solvers Library (LSL), que conté un gran nombre de mètodes per resoldre els sistemes d'equacions lineals derivats de les discretitzacions. El primer capítol d'aquesta tesi conté les principals idees subjacents al disseny i la implementació de la BOL i la LSL, juntament amb alguns exemples i algunes aplicacions industrials. En els capítols posteriors hi ha una explicació detallada de solvers específics per algunes aplicacions concretes. En el segon capítol, es presenta un solver paral·lel i directe per la resolució de l'equació de Poisson per casos en els quals una de les direccions del domini té condicions d'homogeneïtat. En la simulació de fluxos incompressibles, l'equació de Poisson es resol almenys una vegada en cada pas de temps, convertint-se en una de les parts més costoses i difícils de paral·lelitzar del codi. El mètode que proposem és una combinació d'una descomposició directa de Schur (DDS) i una diagonalització de Fourier. La darrera descompon el sistema original en un conjunt de sub-sistemes 2D independents que es resolen mitjançant l'algorisme DDS. Atès que no s'imposen restriccions a les direccions no periòdiques del domini, aquest mètode és aplicable a la resolució de problemes discretitzats mitjançat l'extrusió de malles 2D no estructurades. L'escalabilitat d'aquest mètode ha estat provada amb èxit amb un màxim de 8192 CPU per malles de fins a ~10⁹ volums de control. En el darrer capitol capítol, es presenta un mètode de resolució per l'equació de Transport de Boltzmann (BTE). La estratègia emprada es basa en el mètode d'Ordenades Discretes i pot ser aplicat en discretitzacions no estructurades. El flux per a cada ordenada angular es resol amb un mètode de substitució equivalent a la resolució d'un sistema lineal triangular. La naturalesa seqüencial d'aquest procés fa de la paral·lelització de l'algoritme el principal repte. Diversos algorismes de substitució han estat analitzats, esdevenint una de les heurístiques proposades la millor opció en totes les situacions analitzades, amb excel·lents resultats. Els testos d'eficiència paral·lela s'han realitzat usant fins a 2560 CPU.Direct Numerical Simulation (DNS) of complex flows is currently an utopia for most of industrial applications because computational requirements are too high. For a given flow, the gap between the required and the available computing resources is covered by modeling/simplifying of some terms of the original equations. On the other hand, the continuous growth of the computing power of modern supercomputers contributes to reduce this gap, reducing hence the unresolved physics that need to be attempted with approximated models. This growth, widely relies on parallel computing technologies. However, getting the expected performance from new complex computing systems is becoming more and more difficult, and therefore part of the CFD research is focused on this goal. Regarding to it, some contributions are presented in this thesis. The first objective was to contribute to the development of a general purpose multi-physics CFD code. referred to as TermoFluids (TF). TF is programmed following the object oriented paradigm and designed to run in modern parallel computing systems. It is also intensively involved in many different projects ranging from basic research to industry applications. Besides, one of the strengths of TF is its good parallel performance demonstrated in several supercomputers. In the context of this thesis, the work was focused on the development of two of the most basic libraries that compose TF: I) the Basic Objects Library (BOL), which is a parallel unstructured CFD application programming interface, on the top of which the rest of libraries that compose TF are written, ii) the Linear Solvers Library (LSL) containing many different algorithms to solve the linear systems arising from the discretization of the equations. The first chapter of this thesis contains the main ideas underlying the design and the implementation of the BOL and LSL libraries, together with some examples and some industrial applications. A detailed description of some application-specific linear solvers included in the LSL is carried out in the following chapters. In the second chapter, a parallel direct Poisson solver restricted to problems with one uniform periodic direction is presented. The Poisson equation is solved, at least, once per time-step when modeling incompressible flows, becoming one of the most time consuming and difficult to parallelize parts of the code. The solver here proposed is a combination of a direct Schur-complement based decomposition (DSD) and a Fourier diagonalization. The latter decomposes the original system into a set of mutually independent 2D sub-systems which are solved by means of the DSD algorithm. Since no restrictions are imposed in the non-periodic directions, the overall algorithm is well-suited for solving problems discretized on extruded 2D unstructured meshes. The scalability of the solver has been successfully tested using up to 8192 CPU cores for meshes with up to 10 9 grid points. In the last chapter, a solver for the Boltzmann Transport Equation (BTE) is presented. It can be used to solve radiation phenomena interacting with flows. The solver is based on the Discrete Ordinates Method and can be applied to unstructured discretizations. The flux for each angular ordinate is swept across the computational grid, within a source iteration loop that accounts for the coupling between the different ordinates. The sequential nature of the sweep process makes the parallelization of the overall algorithm the most challenging aspect. Several parallel sweep algorithms, which represent different options of interleaving communications and calculations, are analyzed. One of the heuristics proposed consistently stands out as the best option in all the situations analyzed. With this algorithm, good scalability results have been achieved regarding both weak and strong speedup tests with up to 2560 CPUs

Modelling and simulation of flexible instruments for minimally invasive surgical training in virtual reality

Improvements in quality and safety standards in surgical training, reduction in training hours and constant technological advances have challenged the traditional apprenticeship model to create a competent surgeon in a patient-safe way. As a result, pressure on training outside the operating room has increased. Interactive, computer based Virtual Reality (VR) simulators offer a safe, cost-effective, controllable and configurable training environment free from ethical and patient safety issues. Two prototype, yet fully-functional VR simulator systems for minimally invasive procedures relying on flexible instruments were developed and validated. NOViSE is the first force-feedback enabled VR simulator for Natural Orifice Transluminal Endoscopic Surgery (NOTES) training supporting a flexible endoscope. VCSim3 is a VR simulator for cardiovascular interventions using catheters and guidewires. The underlying mathematical model of flexible instruments in both simulator prototypes is based on an established theoretical framework – the Cosserat Theory of Elastic Rods. The efficient implementation of the Cosserat Rod model allows for an accurate, real-time simulation of instruments at haptic-interactive rates on an off-the-shelf computer. The behaviour of the virtual tools and its computational performance was evaluated using quantitative and qualitative measures. The instruments exhibited near sub-millimetre accuracy compared to their real counterparts. The proposed GPU implementation further accelerated their simulation performance by approximately an order of magnitude. The realism of the simulators was assessed by face, content and, in the case of NOViSE, construct validity studies. The results indicate good overall face and content validity of both simulators and of virtual instruments. NOViSE also demonstrated early signs of construct validity. VR simulation of flexible instruments in NOViSE and VCSim3 can contribute to surgical training and improve the educational experience without putting patients at risk, raising ethical issues or requiring expensive animal or cadaver facilities. Moreover, in the context of an innovative and experimental technique such as NOTES, NOViSE could potentially facilitate its development and contribute to its popularization by keeping practitioners up to date with this new minimally invasive technique.Open Acces

DEVELOPMENT OF A LAGRANGIAN-LAGRANGIAN METHODOLOGY TO PREDICT BROWNOUT DUST CLOUDS

A Lagrangian-Lagrangian dust cloud simulation methodology has been developed to help better understand the complicated two-phase nature of the rotorcraft brownout problem. Brownout conditions occur when rotorcraft land or take off from ground surfaces covered with loose sediment such as sand and dust, which decreases the pilot's visibility of the ground and poses a serious safety of flight risk. The present work involved the development of a comprehensive, computationally efficient three-dimensional sediment tracking method for dilute, low Reynolds number Stokes-type flows. The flow field generated by a helicopter rotor in ground effect operations over a mobile sediment bed was modeled by using an inviscid, incompressible, Lagrangian free-vortex method, coupled to a viscous semi-empirical approximation for the boundary layer flow near the ground. A new threshold model for the onset of sediment mobility was developed by including the effects of unsteady pressure forces that are induced in vortically dominated rotor flows, which can significantly alter the threshold conditions for particle motion. Other important aspects of particle mobility and uplift in such vortically driven dust flows were also modeled, including bombardment effects when previously suspended particles impact the bed and eject new particles. Bombardment effects were shown to be a particularly significant contributor to the mobilization and eventual suspension of large quantities of smaller-sized dust particles, which tend to remain suspended. A numerically efficient Lagrangian particle tracking methodology was developed where individual particle or clusters of particles were tracked in the flow. To this end, a multi-step, second-order accurate time-marching scheme was developed to solve the numerically stiff equations that govern the dynamics of particle motion. The stability and accuracy of this scheme was examined and matched to the characteristics of free-vortex method. One-way coupling of the flow and the particle motion was assumed. Particle collisions were not considered. To help reduce numerical costs, the methodology was implemented on graphic processing units, which gave over an order of magnitude reduction in simulation time without any loss in accuracy. Validation of the methodology was performed against available measurements, including flow field measurements that have been made with laboratory-scale and full-scale rotors in ground effect operations. The predicted dust clouds were also compared against measurements of developing dust clouds produced by a helicopter during taxi-pass and approach-to-touchdown flight maneuvers. The results showed that the problem of brownout is mostly driven by the local action of the rotor wake vortices and the grouping or bundling of vortex filaments near the sediment bed. The possibilities of mitigating the intensity of brownout conditions by diffusing the blade tip vortices was also explored. While other means of brownout mitigation may be possible, enhancing the diffusion of the tip vortices was shown to drastically reduce the quantity of mobilized particles and the overall severity of the brownout dust cloud