Which subroutine does complex Hermitian tridiagonal matrix?

[carolasu]

which subroutine computes all eigenvalues and, optionally, eigenvectors of complex symmetric tridiagonal matrix A?

I used dstev.f for real symmetric tridiagonal matrix, but now I need to use matrix in complex number, I am wondering which subroutine I can call for?

Thank you!

[Julien Langou]

Do you mean complex symmetric tridiagonal or complex Hermitian tridiagonal?

[carolasu]

Thanks. I mean complex Hermitian tridiagonal.

[Julien Langou]

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.

[carolasu]

I see. It's very clear.

(1)complex A stored as real D, complex E.

(2)Use S to scale complex E to real E.

(3)Call dstev.f with real E and real D to get V, W.

(4)Use S to scale real V to complex V.

AA and TT are not needed, but we can use those to check.

Thank you!

[mattcimd]

Does anyone know if there is a theorem, which says that you can transform a hermitian tridiagonalmatrix in a real symmetric tridiagonalmatrix by similarity transform with a unitary diagonalmatrix? And where i can find the theorem?

Thanks for help!!!