Solving large linear systems with multiple right-hand sides.
- Julien LANGOU -
- PhD. Thesis defence -
- Tuesday June 10th, 2003, 17h00/19h00,St Girons (09, France) -

defence for a PhD. in Applied Mathematics of the
``Institut National des Sciences Appliquées de Toulouse''
 
Jury :
Guillaume Alléon EADS-CCR (France)
Åke Björck Linköping University (Sweden) rapporteur
Iain S. Duff CERFACS (France) & RAL (UK)
Luc Giraud CERFACS (France) directeur de thèse
Gene H. Golub Stanford University (CA, USA)
Gérard Meurant CEA (France) rapporteur
Chris C. Paige Mc Gill University (Canada) rapporteur
 
Abstract

  The starting point of this thesis is a problem posed by the electromagnetism group at EADS-CCR: How to solve several linear systems with the same coefficient matrix but various right-hand sides ? For the targeted application, the matrices are complex, dense and huge (of order of a few millions). Because such matrices cannot be computed nor stored in numerical simulations involved in a design process, the use of an iterative scheme with an approximate matrix-vector product is the only alternative. The matrix-vector product is performed using the Fast Multipole Method. In this context, the goal of this thesis is to adapt Krylov solvers so that they handle efficiently multiple right-hand sides. Some preliminary works dealing with one right hand side show that GMRES was an efficient and robust solver for that application. Consequently, we mainly focus, in this thesis, on variants of GMRES.  

  The orthogonalization schemes that we implemented for GMRES are some variants of the Gram-Schmidt algorithm. In a first part, we have investigated the effect of rounding errors in the Gram-Schmidt algorithms. Our results answer questions that have been around for the last 25 years. We give theoretical explanations of what was currently observed and accepted, namely:  

  • modified Gram-Schmidt algorithm generates a well-conditioned set of vectors,
  • and Gram-Schmidt algorithm iterated twice gives an orthogonal set of vectors.
These two sentences holds when the initial matrix is ``not-too-ill-conditioned'' in a sense that is clearly defined. Furthermore, when the Gram-Schmidt algorithm is iterated with selective reorthogonalizations then, in a third part, we give a new criterion. We prove that the resulting algorithm is robust while the most commonly used criterion might have some weakness. In a fourth part, we generalize standard results for modified Gram-Schmidt from norms to singular values. This enables us to propose an a-posteriori reorthogonalization procedure for modified Gram-Schmidt algorithm based on a low rank update. These results have several direct applications, we give examples for Krylov methods for solving linear systems and for eigenvalue computations. Finally, instead of using the Euclidean scalar product, we have derived our results on Gram-Schmidt algorithm with the A-scalar product, where A is an Hermitian definite positive matrix. The relevance of such a study is illustrated in a collaboration with the Global Change team at CERFACS, the Lanczos algorithm with reorthogonalization is used for data assimilation in a climate modeling problem.

  In a second part, we have implemented variants of the GMRES algorithm for both real and complex, single and double precision arithmetics suitable for serial, shared memory and distributed memory computers. This software implementation, that complies with scientific library standard quality, is described in detail. For the sake of the simplicity, flexibility and efficiency the GMRES solvers have been implemented using the reverse communication mechanism for the matrix-vector product, the preconditioning step and the dot product computations. Several orthogonalization procedures have been implemented to reduce the cost of the dot product calculation, that is a well known bottleneck for the Krylov method efficiency in a parallel distributed environment. The implemented stopping criterion is based on a normwise backward error. The variants available are GMRES-DR, seed-GMRES and block-GMRES (that adds to standard GMRES, flexible GMRES and SQMR). An LU-matrix-vector product step is performed in GMRES-DR in order to store the approximate eigenvectors on the first Krylov vectors. Implicit restart and implicit preconditioning is done in seed-GMRES to avoid a matrix-vector product and a preconditioning step per right-hand side and per GMRES cycle. The block-GMRES version allows the user to select either no deflation or deflation based on the singular value decomposition of the Krylov vectors. Finally we extend existing theoretical result that relates the norm of the residual to the smallest singular value of the space constructed by the Krylov solver from GMRES to block-GMRES. 

  The third part is dedicated to the improvement of these standard methods for the solution of linear systems arising in electromagnetic applications. After a deep presentation of the code we use, we study the influence of the level of non-symmetry in the SQMR algorithm and illustrate the behaviour of the GMRES-DR method in our example. Then we mainly focus on the multiple right-hand sides solves. First of all, we examined in details techniques to adapt single right-hand side method in the multiple right-hand sides case: increase the quality of the preconditioner, initial guess strategy or gathering strategy. In the context of the monostatic radar cross section calculation, we prove that the number of independent right-hand sides is finite. The dimension of the space of right-hand sides given by our theory and the dimension numerically observed corresponds well each other. This property enables us to reduce considerably the number of linear systems to solve. In this context, a particular implementation of the block-GMRES method is given. Then, some specific issues about the seed and the block methods are discussed. Finally more prospective results are given. First, different strategies to extract and add spectral informations in a GMRES cycle are presented and compared. Then, we exploit the fact that the Fast Multipole Method is an inexact matrix-vector product for which the accuracy can be tuned. The less accurate the matrix-vector product is, the fastest the computation. We show how to take advantage of this by using relaxed schemes (inexact Krylov methods) or inner-outer schemes (flexible GMRES). Finally we study the relevance of the normwise backward error as a stopping criterion for the iterative solvers in the monostatic radar cross section problem.
 

Résumé

  Le point de départ de cette thèse est un problème posé par le groupe électromagnétisme de EADS-CCR : comment résoudre plusieurs sytèmes linéaires avec la même matrice mais différents seconds membres ? Pour l'application voulue, les matrices sont complexes, denses et de grande taille (de l'ordre de quelques millions). Comme de telles matrices ne peuvent être ni calculées, ni stockées dans un processus industriel, l'utilisation d'un produit matrice-vecteur approché est la seule alternative. En l'occurence, le produit matrice-vecteur est effectué en ultilisant la méthode multipôle rapide. Dans ce contexte, le but de cette thèse est d'adapter les méthodes itératives de type Krylov de telles sorte qu'elles traitent efficacement les nombreux seconds membres. Des travaux préliminaires avec un seul second membre ont montré que la méthode GMRES était particulièrement efficace et robuste pour cette application. En conséquence, nous nous concentrons, dans cette thèse, uniquement sur des variantes de GMRES.

  Les schémas d'orthogonalisation que nous avons implanté dans GMRES sont des variantes de l'algorithme de Gram-Schmidt. Dans une première partie, nous nous intéressons à l'influence des erreurs d'arrondi dans les algorithmes de Gram-Schmidt. Nos resultats répondent à des questions vieilles de vingt-cinq ans. Nous donnons l'explication théorique de ce qui était communément observé et accepté :

  • l'algorithme de Gram-Schmidt modifié génère un ensemble de vecteurs bien conditionné,
  • et l'algorithme de Gram-Schmidt itéré deux fois fabrique un ensemble de vecteurs orthonormals.
Ces deux propositions reposent sur l'hypothèse que la matrice de départ est ``numériquement non singulière' en un sens qui est clairement défini. D'autre part, quand l'agorithme de Gram-Schmidt est itéré avec un critère de reorthogonalisation, dans une troisième partie, nous proposons un nouveau critère. Nous montrons que l'algorithme obtenu est robuste alors que le critère communément utilisé est mis en défaut dans certain cas. Dans une quatrième partie, nous généralisons des normes aux valeurs singulières des résultats standards sur l'algorithme de Gram-Schmidt modifié. Ceci nous permet de dériver un schéma de réorthogonalisation a-posteriori utilisant une matrice de rang faible. Ces résultats ont plusieurs applications directes, nous donnons des exemples avec les méthodes de Krylov pour résoudre des problèmes linéaires avec plusieurs seconds membres. Finalement, au lieu d'utiliser le produit scalaire euclidien, nous avons dérivé nos résultats sur Gram-Schmidt avec le A-produit scalaire, où A est hermitienne définie positive. L'importance d'une telle étude est illustrée par une collaboration avec l'équipe Global Change du CERFACS, l'algorithme de Lanczos avec réorthogonalisation est utilisé pour faire de l'assimilation de donnée dans un problème de modélisation du climat.

  Dans une deuxième partie, nous avons implémenté des variantes de la méthode GMRES pour les arithmétiques réelle et complexes, simple et double précision. Cette implémentation convient pour des ordinateurs classiques, à mémoire partagée où distribuée. Le code résultant satisfait aux critères de qualité des librairies standards et son implémentation est largement détaillée. Pour des besoins de simplicité, flexibilité et efficacité, les solveurs utilisent un mécanisme de reverse communication pour les produits matrice-vecteur, les étapes de préconditionement et les produits scalaires. Différents schémas d'orthogonalisation sont implémentés pour réduire le coût de calcul des produits scalaires, un point particulièrement important pour l'efficacité des méthodes de Krylov dans un environnement parallèle distribué. Le critère d'arrêt implémenté est basé sur l'erreur inverse normalisée. Les variantes disponibles sont GMRES-DR, seed-GMRES et block-GMRES, ces codes s'ajoutent aux variantes déjà existantes (GMRES, flexible GMRES et SQMR). Un produit matrice-vecteur avec une décomposition LU est utilisé dans GMRES-DR de tel sorte que le stockage des approximations des vecteurs propres se fassent sur les premiers vecteurs de l'espace de Krylov. Un restart implicite et une étape de préconditionnement implicite a été implémenté dans seed-GMRES, ceci permet d'éviter un produit matrice-vecteur et une étape de préconditionement par second membre et par cycle de GMRES. La version de block-GMRES permet à l'utilisateur de sélectionner différents modes de déflation. Finalement des résultats reliant la norme du résidu de GMRES à la plus petite valeur singulière de l'espace construit par la méthode de Krylov ont été généralisés à la méthode block-GMRES.

  La troisième partie est dédiée à l'amélioration de ces techniques standards dans le cadre des problèmes électromagnétiques. Après une présentation approfondie du code, nous étudions l'influence sur la convergence du niveau de non-symmétrie dans l'algorithme SQMR. Nous étudions aussi le comportement de GMRES-DR sur nos problèmes. Ceci correspond à deux méthodes avec un seul second membre, le reste de cette partie est dédiée au cas avec plusieurs second membres. Tout d'abord, nous examinons en détail les techniques qui consistent à adapter les méthodes pour un second membre unique à plusieurs seconds membres : améliorer la qualité du préconditioneur, avoir une stratégie de solution initiale ou grouper les opérations. Dans le contexte du calcul de surface équivalente radar monostatique, nous avons montré que l'espace des seconds membres était de dimension finie. La dimension donnée par notre théorie est proche de ce que nous observons en pratique. Cette propriété nous permet de réduire considérablement le nombre de systèmes linéaires à résoudre. Dans ce contexte, une version de la méthode block-GMRES est donnée. Ensuite, nous abordons certains problèmes spécifiques des méthodes seed-GMRES et block-GMRES ; nous proposons des solutions. Pour finir, des résultats plus prospectifs sont donnés. Tout d'abord, plusieurs stratégies pour extraire et ajouter de l'information spectrale d'un cylcle de GMRES à l'autre sont proposées et comparées. Ensuite nous utilisons le fait que le méthode multipôle rapide est un produit matrice-vecteur inexacte dont la précision est règlable. Moins précis est le produit matrice-vecteur, plus rapide il est. Nous montrons comment tirer partie de cette propriété en utilisant un schéma relaché (méthode de Krylov inexacte) ou des itérations emboîtées (flexible GMRES). Finallement, le critère d'arrêt basé sur l'erreur inverse normalisée dans le cadre du calcul d'une surface équivalente radar est remis en question.

 

 
algweb@cerfacs.fr
Last Update: May 21, 2003