我已经实现了 Matrix Computations 4th Edition 一书中描述的 Aasen 算法版本。那里的版本没有旋转。这本书对如何添加旋转的描述对我来说有点粗略。有人可以用实际代码说明如何添加透视吗？

我应该补充一点，我的实现支持 Hermitian 矩阵，这与书中的不同。

from sympy import *

def aasen(A):
    """Aasen without pivoting, but with support for Hermitian matrices."""
    n = A.rows
    alpha = zeros(n,1)
    beta = zeros(n,1)
    L = eye(n)
    h = zeros(n,1)
    v = zeros(n,1)
    for j in range(n):
        if j == 0:
            alpha[0] = A[0,0]
            v[1:n,0] = A[1:n,0]
        else:
            h[0] = beta[0].conjugate() * L[j,1].conjugate()
            for k in range(1,j):
                h[k] = beta[k-1] * L[j, k-1].conjugate() + alpha[k] * L[j,k].conjugate() + beta[k].conjugate() * L[j,k+1].conjugate()
            h[j] = A[j,j] - L[j,0:j].dot(h[0:j])
            alpha[j] = h[j] - beta[j-1] * L[j,j-1].conjugate()
            v[j+1:n,0] = A[j+1:n,j] - L[j+1:n,0:j+1] * Matrix(h[0:j+1])
        if j <= n-2:
            beta[j] = v[j+1]
        if j <= n - 3:
            L[j+2:n,j+1] = Matrix(v[j+2:n]) / v[j+1]
    T = zeros(n,n)
    for i in range(n):
        for j in range(n):
            if i == j+1:
                T[i,j] = beta[j]
            elif i == j:
                T[i,j] = alpha[j]
            elif j == i+1:
                T[i,j] = beta[i].conjugate()
    return L, T