Recursive blocked algorithms for linear systems with Kronecker product structure
- PDF / 1,742,749 Bytes
- 18 Pages / 439.642 x 666.49 pts Page_size
- 82 Downloads / 199 Views
Recursive blocked algorithms for linear systems with Kronecker product structure Minhong Chen1 · Daniel Kressner2 Received: 24 May 2019 / Accepted: 20 August 2019 / © Springer Science+Business Media, LLC, part of Springer Nature 2019
Abstract Recursive blocked algorithms have proven to be highly efficient at the numerical solution of the Sylvester matrix equation and its generalizations. In this work, we show that these algorithms extend in a seamless fashion to higher-dimensional variants of generalized Sylvester matrix equations, as they arise from the discretization of PDEs with separable coefficients or the approximation of certain models in macroeconomics. By combining recursions with a mechanism for merging dimensions, an efficient algorithm is derived that outperforms existing approaches based on Sylvester solvers. Keywords Blocked algorithm · Linear system · Tensor equation · Sylvester equation
1 Introduction In computations with matrices, recursive blocked algorithms offer an elegant way to arrive at implementations that benefit from increased data locality and efficiently utilize highly tuned kernels. See [7] for a survey and [25] for a more recent testimony of this principle. These algorithms have proven particularly effective for solving Sylvester equations, that is, matrix equations of the form A1 X + XAT2 = B, Rn1 ×n1 , A2
Rn2 ×n2 ,
Rn1 ×n2
(1.1) Rn1 ×n2
∈ and B ∈ are given and X ∈ where A1 ∈ is unknown. In the Bartels-Stewart algorithm [2], the matrices A1 and A2 are first Daniel Kressner
[email protected] Minhong Chen [email protected] 1
Department of Mathematics, Zhejiang Sci-Tech University, Hangzhou, 310029, Zhejiang, People’s Republic of China
2
Institute of Mathematics, EPF Lausanne, 1015 Lausanne, Switzerland
Numerical Algorithms
reduced to block upper form by real Schur decompositions. The reduced problem is then solved by a variant of backward substitution. Both stages of the algorithms require O(n3 ) operations, with n = max{n1 , n2 }. Entirely consisting of level 2 BLAS operations, the backward substitution step performs quite poorly. To avoid this, Jonsson and K˚agstr¨om [14, 15] have proposed recursive algorithms for triangular Sylvester and related matrix equations. The recursive algorithm for solving (1.1) with upper quasi-triangular A 1 , A2 starts with partitioning the matrix of larger size. A1,11 A1,12 with A1,11 ∈ Rk×k such that k ≈ n/2 Assuming n1 ≥ n2 , let A1 = 0 A1,22 X1 B1 and partition X = , B = correspondingly. Then (1.1) becomes X2 B2 equivalent to A1,11 X1 + X1 AT2 = B1 − A1,12 X2 ,
(1.2a)
A1,22 X2 + X2 AT2
(1.2b)
= B2 .
First the Sylvester (1.2b) is solved recursively, then the right-hand side (1.2a) is updated, and finally (1.2a) is solved recursively. Apart from the solution of smallsized Sylvester equations at the lowest recursion level, the entire algorithm consists of matrix-matrix multiplications A1,12 X2 and thus attains high performance by leveraging level 3 BLAS. As emphasized in [7, 25], recursive algorithms are less sensiti
Data Loading...
