High performance implementation of block Krylov methods

Abstract

Multiple large linear systems with a same coefficient matrix but different right-hand side vectors arise in several problems of engineering such as in solution of Navier-Stokes equations in different directions (Dinler, 2007). Moreover, seeking an optimum or a special solution of a multi-dimensional PDE with a parameter-dependent moving source term or changing force term, or variable boundary conditions requires repeated solution of this PDE with changing values of the parameter. In that case, we need to solve multiple large linear systems with a same coefficient matrix but different right-hand side vectors. To solve these large block linear systems efficiently, block Krylov subspace methods (in short block Kylov methods) are an important tool. Krylov subspace methods are projection methods that search the solution of a linear system iteratively in the Krylov subspace. Comparing with direct methods such as Gauss elimination or LU decomposition, direct methods are prohibitive to solve large linear systems however Krylov subspace methods are efficient. Moreover in several engineering problems, block Krylov methods are an important tool for solution of large linear systems with multiple right hand sides. In addition to that, block methods reduce memory-access cost and increase efficient data reuse because they allow multiple matrixvector products with less memory access. However, their robustness problem and preconditioning are difficult to handle. So, to come up with a robust implementation of block Krylov methods, we initiated BIM++ (Block Iterative Methods) project, which is devoted to research and development of block Krylov methods to make them easier to use. BIM++ package is a high performance implementation of block Krylov methods written in C/C++ under GPL license, for solution of both symmetric and nonsymmetrical large linear systems with multiple right hand sides. In addition to this BIM++ includes important non-block Krylov methods such as CG (conjugate gradient), BICGSTAB (bi-conjugate gradient stabilized), GMRES (generalized minimal residual) and GMRES(m) (GMRES restarted). While developing BIM++, we stick to “keep-it-easy rule” and put performance first. Easy-to-use and easy-to-install are also important. While developing BIM++, our goals are 1) creating an environment for easy implementation and easy development of advanced numerical linear algebra algorithms, 2) creating a high performance and user friendly package, 3) enable to utilize C++ technologies, 4) enable to programming Krylov methods shorter and more readable, 5) enable to use block Krylov methods handy, 6) creating reusable, portable and easily upgradeable package, 7) enable to use BIM++ in Windows PCs, 8) creating a clear documentation and a webpage. To achieve these goals, we borrowed many good ideas from LAPACK++ and IML++. We built BIM++ over LAPACK. We improved and used CPPLapack interface (Onishi, 2003) for FORTRAN LAPACK to make the package more user friendly. Also we tried to use generic programming to achieve performance. In this technical note, we demonstrate test results of BIM++ package for different matrices and different number of right-hand side vectors. Our results are: • Block Krylov methods may be efficient for linear systems with multiple right-hand sides. • Block Krylov methods decrease the memory access cost as a result convergence can be achieved in less time. • Performance of block Krylov methods depends on performance of matrix-matrix multiplication performance. • Using ATLAS library, which is hardware dependent optimal version of LAPACK, under BIM++ increases performance. • Superiority of block Krylov methods over non-block Krylov methods increases in case there are more right-hand side vectors. Moreover, Convergence may be improved with suitable extra right-hand vectors. • Comparing with non-block Krylov methods, convergence behavior of block Krylov methods are different so different and special preconditioning is necessary.

Özet

Kısmi türevli bir denklemde parametreye bağlı değişen sınır koşulları, hareketli bir kaynak terimi ya da değişken kuvvet terimi ile özel veya optimize bir sayısal çözümün aranması durumunda denklemin bu parametrenin değişen değerlerine göre çözülmesi gerekir. Bu da sağ-taraf vektörü değişen bir lineer sistemin birçok defa aynı katsayı matrisi ile çözümünü gerektirir. Böyle bir problem üçboyutta ise problemin sayısal çözümünde ortaya çıkan katsayı matrisi büyüktür ve bu durumda uygun iteratif metot kullanılması önemlidir. Bu gibi problemlerde katsayı matrisi değişmediği halde, m tane sağ taraf ile lineer sistemin m defa çözülmesi gerekir. Bunun yerine AX = B blok lineer sistemi oluşturulur, burada A katsayı matrisi, X bilinmeyenler matrisi ve B = [$b_{1,}b_{2,...,}b_m$ ] sağ taraf vektörlerinden oluşan matristir. Bu blok sistem daha verimli ve hızlı bir şekilde blok Krylov metotları ile bir defada çözülebilir. Geliştirmekte olduğumuz Blok İteratif Metotlar paketi (BİM++), simetrik ve simetrik-olmayan blok lineer sistemleri çözen blok Krylov metotlarının yüksek başarımlı uygulamasıdır. Bu çalışmada, blok Krylov metotları ve özellikle blok GMRES (generalized minimal residual) metodu kısaca tanıtıldıktan sonra bu metotlar için geliştirmekte olduğumuz BİM++ paketinin performansı gösterildi. Dahası bazı örnek matrisler üzerinde, çok sağ-taraflı sistemlerin çözümünde blok metotların blok-olmayan metotlara göre üstünlük gösterebildiği ve sağ taraf sayısının artması ile bu blok-olmayan metotlara olan üstünlüğün arttığı gösterildi.