块Arnoldi约化

设 A 是一个阶数为 n 的矩阵，且 b > 0 为块大小。我们称

A V_{[m]} = V_{[m]} H_{[m]} + F_m E^{\ast}_m \tag{7.29}

为长度为 m 的块Arnoldi约化，当 V_{[m]}^\ast AV_{[m]}=H_{[m]} 是一个阶数为 m\cdot b 的带状上Hessenberg矩阵， V_{[m]}^\ast V_{[m]}= I_{m b}，且 V_{[m]}^\ast F_m=0。 对于我们来说，带状上Hessenberg矩阵是一个具有 b 条次对角线的上三角矩阵。 V_{[m]} = [V_1,V_2,\ldots,V_m] 的列构成了块Krylov子空间

{\mathcal K}^m(A,V_1) \equiv \mathrm{span}\{V_1, A V_1,\ldots, A^{m-1} V_1 \}.

如果 m > \bar{m} \equiv \mathrm{ceiling}(n/b)，那么 F_m = {0} 且 H_{[m]} 是 A 到带状上Hessenberg形式的正交约化。 我们暂时假设 F_m 是满秩的，并且进一步假设 H_{m+1,m} 对角线上的元素是正的。 因此，隐式Q定理的直接扩展 [198, pp. 367-368] 给出了 F_m 由初始块 V_1 唯一指定。 注意，如果 A=A^\ast，那么 H_{[m]} 是一个块三对角矩阵。 算法 7.11 列出了计算块Arnoldi约化的算法。

算法 7.11：扩展块Arnoldi约化
(1) 使用迭代CGS方法计算QR分解 V_{m+1} H_{m+1, m} = F_m
(2) V_{[m+1]} = \left[\begin{array}{ll}V_{[m]} & V_{m+1}\end{array}\right]
(3) W = AV_{m+1}
(4) H_{m+1, m+1} = V_{m+1}^* W
(5) H_{[m+1]} = \left[\begin{array}{cc} H_{[m]} & V_{[m]}^* W \\ H_{m-1,m}E^*_m & H_{m+1,m+1} \end{array}\right]
(6) F_{m+1} = W - V_{[m+1]} \begin{bmatrix} V^*_{[m]}W \\ H_{m+1, m+1}\end{bmatrix}

我们现在对算法 7.11 的一些步骤进行评论。

(1)

这里通过迭代CGS算法并使用可能的校正步骤来计算QR分解。 详见 [96]，以及用于确定是否需要校正步骤的简单测试。 这种方案的一个好处是它允许使用Level 2 BLAS [133] 矩阵-向量乘法子程序 xGEMV（另见§10.2）。 此外，这种方案还提供了一种简单的方法来填充秩不足的 F_m。例如， 如果在生成 V_{m+1} 的第 j 列时需要第三次正交化步骤，那么 F_m 的相应列在线性上依赖于 V_{m+1} 的前 j-1 列。 H_{m+1,m} 的第 j 个对角元素设置为零，并且一个随机单位向量相对于 V_{[m]} 和 V_{m+1} 的前 j-1 列进行正交化。

(3)

这一步允许对一组向量应用 A。当访问 A 的成本很高时，这可能至关重要。显然，目标是摊销对多个向量应用 A 的成本。

(4)

这允许使用Level 3 BLAS [134] 矩阵-矩阵乘法子程序 xGEMM 来计算 V_m^\ast W。

(6)

这是块经典Gram-Schmidt（bCGR）的一步。 Level 3 BLAS [134] 矩阵-矩阵乘法子程序 xGEMM 用于计算所需的秩-b 更新，以便计算 F_{m+1}。 为了确保 F_{m+1} 与 V_{m+1} 的正交性，除了 b=1 的情况外，执行第二步 bCGR。在后一种情况下，使用DGKS [96] 中的简单测试来确定是否需要第二正交化步骤。详见 [295]。

计算 V_{[m+1]} 的方案等同于Ruhe [375]（对于对称矩阵）提出的方案，并且是由隐式重启块Lanczos代码 [24] 使用的方案。 尽管 [24] 中的方法干净地处理了秩不足的 F_m 问题，但实现并没有利用如（3）中那样应用 A 的能力。相反，如 [375] 所提议的，A 应用于 F_m 的每一列，然后计算相应的 V_{m+1} 和 H_{[m+1]} 列。我们的实现进一步重组了Ruhe的方法，使得系数矩阵 V_{[m]}^\ast W 的计算与 F_m 的QR分解分开。这样做的好处是步骤（3）-（5）通过块大小减少了I/O成本，并增加了每次内存引用的浮点运算量。