This is easy then. For any complex Hermitian tridiagonal matrix, there exists a symmetric scaling (so scaling of the rows and the columns) such that the matrix becomes a real symmetric tridiagonal matrix. In another word, given A, complex Hermitian tridiagonal matrix, there exists S, unitary complex diagonal, such that T=S*A*conj(S) is real symmetric tridiagonal matrix. Note that S diagonal unitary simply means that the diagonal entries have modulus 1. So actually T = S*A*conj(S) is a similarity transformation so A and T have the same eigenvalues.
Therefore, you quickly convert your input A to S, then you use LAPACK real symmetric tridiagonal eigensolver (dstev, dstemr, dstedc, dsteqr, ... ) to get T = V * W * V^H (where V is the matrix of eigenvectors, it is real orthogonal, and W is the matrix of eigenvalues, it is real diagonal) and you do not forget to scale back the V to get your true eigenvectors: A = ( conj(S)*V ) * W * ( conj(S)*V )^H. W stays real, the eigenvectors conj(S)*V becomes complex. conj(S)*V is a unitary matrix.
(Note: this is the reason why there is no complex Hermitian tridiagonal eigensolver. The real symmetric tridiagonal eigensolvers are all that are needed.)
Cannot find any reference for the algo online so here is a quick demo in matlab.
- Code: Select all
clear;
n=10;
% A is complex Hermitian tridiagonal
% create D diagonal of A and E sub-diagonal of A
E(1:n-1,1)=randn(n-1,1)+i*randn(n-1,1);
D(1:n,1)=randn(n,1);
% create AA n-by-n from D and E (needed for check only)
AA=zeros(n); for i=1:n, AA(i,i)=D(i); end; for i=1:n-1, AA(i,i+1)=E(i); AA(i+1,i)=conj(E(i)); end;
% compute S a unitary complex diagonal matrix such that T = S*A*S' complex symmetric tridiagonal
% note that S diagonal unitary simply means that the diagonal entries have modulus 1
% this also means that S'*S = I so actually T = S*A*S' is a similarity transformation so A and T
% have the same eigenvalues
% compute at the same time the matrix T
% the matrix T is stored on D and E and replace the initial matrix A
S=zeros(n,1);
S(1) = 1;
if ( abs(E(1)) == 0 )
S(2) = 1;
else
S(2) = E(1)/abs(E(1));
E(1) = abs(E(1));
end
for i=2:n-1,
if ( abs(E(i)) == 0 )
S(i+1) = 1;
else
S(i+1) = S(i)*E(i)/abs(E(i));
E(i) = abs(E(i)); % ( note that E is real now )
end;
end;
% create TT n-by-n symmetric real from D and E
% this is needed in Matlab to call the eigensolver
% in Fortran: not needed: call directly dstemr, dstedc, dsteqr etc. directly with D and E
TT=zeros(n); for i=1:n, TT(i,i)=D(i); end; for i=1:n-1, TT(i,i+1)=E(i); TT(i+1,i)=E(i); end;
% calling the eigensolver to obtain eigenvectors V and associated eigenvalues W of T
% Note: T is real symmetric tridiagonal
[V,W] = eig(TT);
% backtransform the eigenvectors V of T to obtain the eigenvectors of A
% so perform: V <- S' * V
for i=2:n,
V(i,1:n) = conj(S(i))*V(i,1:n);
end
% check
fprintf('|| A - V W V^H || / || A || = %2.2e \n', norm( AA - V*W*V' )/norm( AA ));
fprintf('|| I - V V^H || = %2.2e \n', norm( eye(n) - V*V' ));
( In a real code, you work with array of size n for D, array of size n-1 for E and array of size n for S. No need to store n-by-n matrices. )
Cheers,
Julien.