Parallel computational methods for solving large linear systems

Abstract

In this thesis, numerical methods, based on the algebraic domain decomposition techniques, have been developed for solving large linear systems on high performance computing (HPC) infrastructures. In particular, a class of algebraic domain decomposition preconditioning schemes, using semi-aggregated subdomains, is proposed. Moreover, an interface semi-aggregation process has been proposed for solving the Schur complement linear system. Additionally, a novel two-stage domain decomposition method has been developed for solving large linear systems on hybrid parallel systems. Algebraic domain decomposition preconditioning schemes have been designed, along with a load balancing algorithm, for distributed memory systems with accelerators, such as manycore processors or GPGPUs. The performance of the algebraic domain decomposition preconditioning schemes has been examined in a cloud computing environment. Furthermore, a time-domain decomposition method in conjunction with semi-aggregated subdomains has been proposed. In Chapter 1, an introduction of the basic mathematical elements in the fields of linear algebra and partial differential equations are given. Further, preconditioned iterative methods for solving linear systems are introduced including domain decomposition methods. Moreover, an introduction of the modern parallel systems is provided, while, the scientific computing software libraries and environments, which have been used, are presented. In Chapter 2, the hybrid parallel Multiprojection method (MPM) for solving large linear systems, using semi-aggregated subdomains, is presented. The class of Multiprojection methods are classified in the algebraic domain decomposition methods, where the local linear systems contain both fine and aggregated (coarse) unknowns. In addition, a subspace compression algorithm has been developed for reducing the memory requirements of the Multiprojection method, namely (MPMSC). This algorithm is based on a breadth first search procedure with limited depth, which is performed on the graph corresponding to the already aggregated components. Neighborhoods of aggregated components are resulted by the aforementioned algorithm and each neighborhood is re-aggregated into one component. In Chapter 3, a hybrid parallel Schur complement solver, using an interface semi-aggregation procedure, has been proposed. Additionally, a hybrid parallel two-stage algebraic domain decomposition method, designed appropriately for distributed memory systems with multicore processors, is proposed. The two-stage Multiprojection method (2s-MPM) is based on upper Schur complement techniques, applied to the local linear systems, and two semi-aggregation procedures, one for each stage. The subdomains are assigned to the workstations of the distributed memory system and the first semi-aggregation procedure is performed in order to formulate the local semi-aggregated linear systems related to the workstations. Upper Schur complement techniques have been used for separating the fine components from the aggregated components of each local semi-aggregated linear system. The second semi-aggregation procedure is performed on the fine components of each workstation, leading to the final local linear systems related to the cores of each multicore processor. Furthermore, a subspace compression technique has been combined with the two-stage Multiprojection method (2s-MPMSC), in order to reduce the memory requirements of the method. The aggregated components outside the workstation are re-aggregated by the subspace compression algorithm, which is based on a breadth first search process with limited depth. In Chapter 4, the Multiprojection preconditioning scheme has been re-designed, using preconditioned iterative methods for solving the local linear systems, in order to be appropriate for distributed memory systems with accelerators. In particular, the preconditioned conjugate gradient method in conjunction with the symmetric factored approximate sparse inverse (SFASI) matrix has been used for solving the local linear systems on an Intel Xeon Phi cluster. The SFASI matrices are computed in the preprocessing phase of the preconditioning scheme. The semi-aggregated subdomains are assigned so that each multicore processor and each manycore processor is responsible for solving one local linear system. Moreover, flexible Krylov subspace methods have been used along with the re-designed iterative MPM scheme for solving linear systems on GPU clusters. The flexible GMRES(m) method and the inexact preconditioned conjugate gradient method for unsymmetric and symmetric coefficient matrices, respectively, have been implemented in order to exploit the available computational resources of a GPU cluster. The local iterative schemes include the preconditioned BiCGSTAB (PBiCGSTAB) method enhanced with the modified generic factored approximate sparse inverse (MGenFASpI) matrix for solving the unsymmetric linear systems. Whereas, the preconditioned conjugate gradient method in conjunction with the SFASI matrix is used for solving the symmetric cases. Additionally, a load balancing algorithm has been designed and implemented for distributing the workload among the available multicore processors and GPUs, aiming to exploit simultaneously the two types of computational resources. Hence, each core of the multicore processor is responsible for the computations concerning one local linear system, while, w local linear systems are assigned to each GPU. Furthermore, the various phenomena, related to performance degradation on cloud computing infrastructures, have been examined by conducting numerical experiments on cloud environment. The "noisy neighbor" phenomenon, where multiple virtual processors that are part of the same physical processor are deployed to different users, is a major factor for degradation of the performance on the cloud. The bandwidth and the latency of the software defined virtual networks is also related to the performance. The hardware heterogeneity in the virtual computational resources, for instance the different clock frequency, cache memory and SIMD register units, often cause performance issues. In Chapter 5, a hybrid parallel time-domain decomposition numerical scheme has been designed and implemented for solving parabolic partial differential equations. The proposed method is based on discretization of the time variable with the backward differences or the Crank-Nicolson method, whereas, the spatial variables are discretized by the finite differences method. The linear equations of multiple time steps are grouped in a sparse linear system. Solving the aforementioned linear system in parallel leads to the solution of the parabolic partial differential equation in multiple time steps, concurrently. The parallel preconditioned GMRES(m) (PPGMRES(m)) method in conjunction with the Multiprojection method is used for solving the aforementioned linear system. The local linear systems are solved by a direct method. Solving multiple time steps in parallel enables the better exploitation of supercomputing infrastructures, on the contrary to the traditional methods that compute one time step per iteration. Two different model problems have been solved in order to assess the performance of the proposed scheme on supercomputing infrastructures: a) a two dimensional thin composite plate cooling model problem and b) a heat transfer model problem in three space variables. In Chapter 6, numerical results for the proposed numerical schemes are given. Numerical experiments concerning the scalability of the proposed schemes have been conducted on high performance computing systems, such as the ARIS HPC supercomputing infrastructure, provided by the GRNET. Moreover, the convergence behavior and the performance of the proposed methods have been examined for several model problems. Comparative results with standard methods, proposed by other researchers, are presented. Finally, concluding remarks are provided concerning the proposed schemes as well as propositions for future research.

Στην παρούσα διδακτορική διατριβή αναπτύχθηκαν αριθμητικές μέθοδοι, οι οποίες βασίζονται σε τεχνικές αλγεβρικού διαχωρισμού του χωρίου, για την επίλυση γραμμικών συστημάτων μεγάλης τάξης σε υπολογιστικά συστήματα υψηλής απόδοσης. Συγκεκριμένα, έχει προταθεί μία κλάση προσυντονισμένων σχημάτων αλγεβρικού διαχωρισμού του χωρίου, που χρησιμοποιεί ημι-συσσωματωμένα υποχωρία. Επιπροσθέτως, έχει προταθεί μία διαδικασία ημι-συσσωμάτωσης των εσωτερικών συνόρων για την επίλυση του γραμμικού συστήματος του συμπληρώματος Schur. Επιπλέον, έχει αναπτυχθεί μία νέα μέθοδος διαχωρισμού του χωρίου με δύο στάδια για την επίλυση γραμμικών συστημάτων μεγάλης τάξης σε υβριδικά παράλληλα συστήματα. Έχουν σχεδιαστεί προσυντονισμένα σχήματα αλγεβρικού διαχωρισμού του χωρίου, μαζί με έναν αλγόριθμο ισοκατανομής του υπολογιστικού φόρτου, για συστήματα κατανεμημένης μνήμης με επιταχυντές, όπως επεξεργαστές πολλών πυρήνων manycore ή κάρτες γραφικών. Έχει μελετηθεί η απόδοση του προτεινόμενου σχήματος σε περιβάλλον υπολογιστικού νέφους. Επίσης, έχει προταθεί μία μέθοδος διαχωρισμού του χωρίου και του χρόνου σε συνδυασμό με υποχωρία ημι-συσσωμάτωσης. Στο κεφάλαιο 1, γίνεται εισαγωγή στις βασικές μαθηματικές έννοιες των πεδίων της γραμμικής άλγεβρας και των μερικών διαφορικών εξισώσεων. Ακόμα, γίνεται μία ανασκόπηση των προσυντονισμένων επαναληπτικών μεθόδων για την επίλυση γραμμικών συστημάτων, συμπεριλαμβανομένων των μεθόδων διαχωρισμού του χωρίου. Επιπλέον, γίνεται εισαγωγή στα σύγχρονα παράλληλα συστήματα, ενώ παρουσιάζονται οι βιβλιοθήκες λογισμικού για επιστημονικούς υπολογισμούς που χρησιμοποιήθηκαν. Στο κεφάλαιο 2, παρουσιάζεται η υβριδική παράλληλη μέθοδος πολλαπλών προβολών (MPM) για την επίλυση γραμμικών συστημάτων μεγάλης τάξης με τη χρήση ημι-συσσωματωμένων υποχωρίων. Η κλάση των μεθόδων πολλαπλών προβολών κατατάσσεται στις μεθόδους αλγεβρικού διαχωρισμού του χωρίου, όπου τα τοπικά γραμμικά συστήματα περιέχουν αγνώστους υψηλής ακρίβειας καθώς και συσσωματωμένους (χαμηλής ακρίβειας) αγνώστους. Επιπλέον, έχει σχεδιαστεί ένας αλγόριθμος συμπίεσης του υποχωρίου για την μείωση των απαιτήσεων σε μνήμη της μεθόδου πολλαπλών προβολών (MPMSC). Ο αλγόριθμος βασίζεται στη διαδικασία της κατά πλάτος αναζήτησης με περιορισμένο βάθος, η οποία εφαρμόζεται στον γράφο που αντιστοιχεί στα ήδη συσσωματωμένα στοιχεία. Σχηματίζονται γειτονιές από συσσωματωμένα στοιχεία από τον προαναφερθέντα αλγόριθμο και κάθε γειτονιά επανασυσσωματώνεται δημιουργώντας ένα νέο στοιχείο. Στο κεφάλαιο 3, έχει προταθεί ένας υβριδικός παράλληλος επιλυτής του συμπληρώματος Schur που χρησιμοποιεί μια διαδικασία ημι-συσσωμάτωσης των εσωτερικών συνόρων. Επιπλέον, προτείνεται μία υβριδική παράλληλη μέθοδος αλγεβρικού διαχωρισμού του χωρίου με δύο στάδια, η οποία σχεδιάστηκε κατάλληλα για συστήματα κατανεμημένης μνήμης με πολυπύρηνους επεξεργαστές. Η μέθοδος πολλαπλών προβολών με δύο στάδια (2s-MPM) βασίζεται σε τεχνικές του άνω συμπληρώματος Schur, οι οποίες εφαρμόζονται στα τοπικά γραμμικά συστήματα, και σε δύο διαδικασίες ημι-συσσωμάτωσης, μία για κάθε στάδιο. Αρχικά τα υποχωρία ανατίθενται στους υπολογιστικούς κόμβους του συστήματος κατανεμημένης μνήμης και η πρώτη διαδικασία ημι-συσσωμάτωσης εφαρμόζεται ώστε να σχηματιστούν τα τοπικά ημι-συσσωματωμένα γραμμικά συστήματα που σχετίζονται με του υπολογιστικούς κόμβους. Χρησιμοποιήθηκαν τεχνικές του άνω συμπληρώματος Schur ώστε να διαχωριστούν τα στοιχεία υψήλης ακρίβειας από τα συσσωματωμένα στοιχεία του κάθε ημι-συσσωματωμένου γραμμικού συστήματος. Η δεύτερη διαδικασία ημι-συσσωμάτωσης σχετίζεται με του πυρήνες του κάθε πολυπύρηνου επεξεργαστή. Επιπλέον, μία τεχνική συμπίεσης του υποχωρίου έχει συνδυαστεί με τη μέθοδο πολλαπλών προβολών δύο σταδίων (2s-MPMSC), ώστε να μειωθούν οι απαιτήσεις σε μνήμη. Τα συσσωματωμένα στοιχεία των άλλων υπολογιστικών κόμβων επανασυσσωματώνονται από τον αλγόριθμο συμπίεσης του υποχωρίου που βασίζεται σε μία διαδικασία της κατά πλάτους αναζήτησης με περιορισμένο βάθος. Στο κεφάλαιο 4, το σχήμα προσυντονισμού πολλαπλών προβολών έχει επανασχεδιαστεί χρησιμοποιώντας προσυντονισμένες επαναληπτικές μεθόδους για την επίλυση των τοπικών γραμμικών συστημάτων, ώστε να είναι κατάλληλο για συστήματα κατανεμημένης μνήμης με επιταχυντές. Συγκεκριμένα, έχει χρησιμοποιηθεί η προσυντονισμένη μέθοδος συζυγών διευθύνσεων σε συνδυασμό με τον συμμετρικό παραγοντοποιημένο προσεγγιστικό αραιό αντίστροφο πίνακα (SFASI) για την επίλυση των τοπικών γραμμικών συστημάτων σε ένα παράλληλο σύστημα με Intel Xeon Phi. Οι πίνακες SFASI υπολογίζονται στη φάση προ-επεξεργασίας του προσυντονισμένου σχήματος. Τα ημι-συσσωματωμένα υποχωρία ανατίθενται κατά τέτοιο τρόπο ώστε κάθε πολυπύρηνος επεξεργαστής και κάθε επιταχυντής να αναλάβει την επίλυση ενός τοπικού γραμμικού συστήματος. Επιπλέον, έχουν χρησιμοποιηθεί οι flexible μέθοδοι του υποχώρου Krylov σε συνδυασμό με το επανασχεδιασμένο επαναληπτικό σχήμα πολλαπλών προβολών για την επίλυση γραμμικών συστημάτων σε παράλληλα συστήματα με κάρτες γραφικών. Έχουν υλοποιηθεί η flexible GMRES(m) (FGMRES(m)) μέθοδος και η ανακριβής προσυντονισμένη μέθοδος συζυγών διευθύνσεων για μη συμμετρικούς και συμμετρικούς συντελεστές πίνακες αντίστοιχα, με σκοπό την αξιοποίηση των διαθέσιμων υπολογιστικών πόρων του παράλληλου συστήματος με κάρτες γραφικών. Τα τοπικά επαναληπτικά σχήματα αποτελούνται από την προσυντονισμένη μέθοδο BiCGSTAB (PBiCGSTAB) σε συνδυασμό με τον τροποποιημένο γενικό παραγοντοποιημένο προσεγγιστικό αντίστροφο πίνακα (MGenFASpI) ή την προσυντονισμένη μέθοδο συζυγών διευθύνσεων σε συνδυασμό με τον πίνακα SFASI, για μη συμμετρικούς και συμμετρικούς συντελεστές πίνακες, αντίστοιχα. Επιπροσθέτως, έχει σχεδιαστεί και υλοποιηθεί ένας αλγόριθμος καταμερισμού του φόρτου για τον διαμοιρασμό του υπολογιστικού φόρτου ανάμεσα στους διαθέσιμους πολυπύρηνους επεξεργαστές και στις κάρτες γραφικών, ώστε να αξιοποιηθούν ταυτόχρονα οι δύο κατηγορίες υπολογιστικών πόρων. Επομένως, κάθε πολυπύρηνος επεξεργαστής αναλαμβάνει τους υπολογισμούς ενός τοπικού γραμμικού συστήματος, ενώ, κάθε κάρτα γραφικών αναλαμβάνει w υποχωρία. Ακόμα, έχουν μελετηθεί ποικίλα φαινόμενα που σχετίζονται με την μείωση της απόδοσης της μεθόδου λόγω της χρήσης υποδομών υπολογιστικού νέφους (cloud computing). Χαρακτηριστικό παράδειγμα αποτελεί το φαινόμενο του "θορυβώδους γείτονα'' (noisy neighbor), όπου οι εικονικοί επεξεργαστές μοιράζονται άνισα τους διαθέσιμους υπολογιστικούς πόρους του υπολογιστικού νέφους. Επίσης, το εύρος ζώνης και ο χρόνος απόκρισης του εικονικού δικτύου μπορεί να αποφέρουν μείωση στην απόδοση της μεθόδου σε περιβάλλον υπολογιστικού νέφους. Η ετερογένεια του υλικού των εικονικών υπολογιστικών πόρων, παραδείγματος χάριν η διαφορετική συχνότητα του επεξεργαστή, η κρυφή μνήμη (cache) και οι εντολές "μια εντολή - πολλαπλά δεδομένα'' (SIMD), συχνά προκαλούν ζητήματα απόδοσης. Στο κεφάλαιο 5, έχει σχεδιαστεί και υλοποιηθεί ένα υβριδικό παράλληλο αριθμητικό σχήμα διαχωρισμού του χωρίου και του χρόνου για την επίλυση παραβολικών μερικών διαφορικών εξισώσεων. Η προτεινόμενη μέθοδος βασίζεται στη διακριτοποίηση της χρονικής μεταβλητής με προς τα πίσω διαφορές ή με τη μέθοδο Crank-Nicolson, ενώ, οι χωρικές μεταβλητές διακριτοποιήθηκαν από τη μέθοδο των πεπερασμένων διαφορών. Ένα μεγάλο γραμμικό σύστημα σχηματίζεται από τη διακριτοποίηση, στην οποία συμμετέχουν πολλαπλά χρονικά βήματα. Eπομένως, περισσότερα από ένα χρονικά βήματα υπολογίζονται ταυτόχρονα με την επίλυση του. Η παράλληλη προσυντονισμένη μέθοδος PPGMRES(m) σε συνδυασμό με την μέθοδο πολλαπλών προβολών χρησιμοποιείται για την επίλυση του προαναφερθέντος γραμμικού συστήματος. Τα τοπικά γραμμικά συστήματα επιλύονται από άμεση μέθοδο. Η επίλυση πολλαπλών χρονικών βημάτων παράλληλα επιτρέπει την καλύτερη αξιοποίηση των υπερυπολογιστικών υποδομών, σε αντίθεση με τις κλασσικές μεθόδους που επιλύουν ένα γραμμικό σύστημα για κάθε χρονικό βήμα. Δύο διαφορετικά προβλήματα έχουν επιλυθεί ώστε να αξιολογηθεί η απόδοση του προτεινόμενου σχήματος σε υπερυπολογιστικές δομές: α) ένα δισδιάστατο πρόβλημα λεπτής πλάκας από σύνθετα υλικά και β) ένα τρισδιάστατο πρόβλημα μεταφοράς θερμότητας. Στο κεφάλαιο 6, δίνονται αριθμητικά αποτελέσματα για τα προτεινόμενα αριθμητικά σχήματα. Έχουν πραγματοποιηθεί αριθμητικά πειράματα σχετικά με την κλιμάκωση των προτεινόμενων μεθόδων σε υπολογιστικά συστήματα υψηλής απόδοσης, όπως η υπερυπολογιστική υποδομή ARIS HPC του Εθνικού Δικτύου Υποδομών Τεχνολογίας και Έρευνας. Επιπλέον, έχει μελετηθεί η συμπεριφορά σύγκλισης και η απόδοση των προτεινόμενων μεθόδων σε διάφορα προβλήματα που προέρχονται από μερικές διαφορικές εξισώσεις. Περαιτέρω, παρουσιάζονται συγκριτικά αποτελέσματα με κλασσικές μεθόδους που έχουν προταθεί από άλλους ερευνητές. Τέλος, παρατίθενται συμπερασματικές παρατηρήσεις σχετικά με τα προτεινόμενα υπολογιστικά σχήματα καθώς και προτάσεις για περαιτέρω μελλοντική έρευνα.