算法模板

包含重启机制和收缩技术以计算多个特征对（eigenpairs）的完整Jacobi-Davidson方法被称为JDQR[172]，因为它可以被视为QR算法的一种迭代方法。该算法的模板在算法4.17中给出。

要应用此算法，我们需要指定初始向量{v}_0、容差\epsilon、目标值\tau以及指定应计算多少个接近\tau的特征对的数值{k}_{\max}。{m}_{\max}表示搜索子空间的最大维度。如果超过此值，将使用指定维度{m}_{\min}的子空间进行重启。

通常，当\tau选择大于\lambda_{\max}(A)时，算法会输出{k}_{\max}个最大特征值；当\tau选择小于\lambda_{\min}(A)时，则会输出{k}_{\max}个最小特征值。计算得到的特征对(\widetilde\lambda_j,\widetilde{x}_{j})，满足\Vert\widetilde{x}_{j}\Vert _2=1，满足\Vert A \widetilde{x}_{j}-\widetilde\lambda_j\widetilde{x}_{j}\Vert _2\leq j\epsilon，其中\widetilde{x}_j表示\widetilde{X}的第j列。

原则上，该算法计算最接近指定目标值\tau的{k}_{\max}个特征值。这仅在需要{k}_{\max}个最大或最小特征值时可靠。对于内部特征值集，我们将在第4.7.4节中描述更安全的技术。现在，我们将根据前几节的讨论，对算法的一些部分进行评论。

(2)

初始化阶段。搜索子空间以t=v_0初始化。

(4)-(6)

搜索子空间的新扩展向量通过修正的Gram-Schmidt方法与当前搜索子空间正交化。为了提高数值稳定性，可以使用算法4.14中给出的模板进行替换。

如果{m}=0，这是一个空循环。

(8)-(10)

我们仅计算厄米矩阵{M}\equiv V^\ast AV（阶数为{m}）的上三角部分。

(11)

对于{m}阶矩阵{M}的特征问题，可以使用LAPACK中的标准稠密厄米特征问题求解器来解决。我们选择计算标准Ritz值，这使得算法适合于计算A的最大或最小{k}_{\max}个特征值。如果希望计算谱内部某处的{k}_{\max}个特征值，建议使用谐波Ritz值；参见第4.7.4节。

矩阵{V}表示列向量为{v}_j的n阶{m}矩阵，{V^A}\equiv AV亦然；S是列向量为s_j的{m}阶矩阵，\Theta = \mathrm{diag}(\theta_1,\ldots,\theta_m)。

(13)

停止准则是在归一化特征向量近似的残差范数低于\epsilon时接受该近似。这意味着我们接受特征值计算中的误差约为\epsilon^2，特征向量中的角度误差约为\epsilon，前提是相关特征值是简单且与其他特征值充分分离的；参见（4.4）。

偶尔，A的某个所需特征向量可能未被检测到，例如，如果v_0在相应特征向量方向上没有分量。对于随机初始向量，这种情况不太可能发生。（参见算法4.13的注释(14)。）

(16)

在接受一个Ritz对后，我们继续搜索下一个特征对，剩余的Ritz向量作为初始搜索空间的基础。这些向量在(17)-(20)中计算。

(23)

一旦当前特征向量的搜索空间维度超过{m}_{\max}，我们就重启。重启过程使用与目标值\tau最接近的{m}_{\min}个Ritz向量所张成的子空间（这些向量在(25)-(27)行中计算）。

(30)-(31)

我们将已锁定的（已计算的）特征向量收集在\widetilde{X}中，矩阵\widetilde{Q}是\widetilde{X}扩展了当前特征向量近似{u}。这样做是为了获得更紧凑的表达形式；修正方程等价于（4.50）中的形式。新的修正{t}必须与\widetilde{X}的列以及{u}正交。

当然，修正方程可以通过任何合适的迭代过程近似求解，例如预处理Krylov子空间方法。由于\widetilde{Q}的出现，在使用预处理矩阵A-\theta I时需要谨慎。预处理器的引入可以遵循与单向量Jacobi-Davidson算法相同的原则；参见算法4.18的模板。确保迭代求解器的初始向量{t}_0满足正交约束\widetilde{Q}^\ast {t}_0=0。注意，如果在几次Jacobi-Davidson迭代中保持{K}不变，算法4.18每步可以节省显著的计算量。在这种情况下，{\widehat{Q}}的列可以从之前的步骤中保存。同样，算法4.18中的矩阵{\mathcal M}及其{\mathcal L}{\mathcal U}分解也可以从之前的步骤中更新。