I've implemented the version of Aasen's algorithm described in the book Matrix Computations 4th Edition. The version there doesn't have pivoting. The book's description of how to add pivoting is a bit sketchy for me. Can someone illustrate with actual code how to add pivoting?
I should add that my implementation has support for Hermitian matrices, unlike the one in the book.
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 = A[0,0] v[1:n,0] = A[1:n,0] else: h = beta.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