Which subroutine does complex Hermitian tridiagonal matrix?

Open discussion regarding features, bugs, issues, vendors, etc.

Which subroutine does complex Hermitian tridiagonal matrix?

Postby carolasu » Wed Aug 15, 2012 12:42 am

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!
Last edited by carolasu on Wed Aug 15, 2012 1:09 am, edited 1 time in total.
carolasu
 
Posts: 9
Joined: Wed Mar 21, 2012 12:50 am

Re: Which subroutine does complex symmetric tridiagonal matr

Postby Julien Langou » Wed Aug 15, 2012 12:46 am

Do you mean complex symmetric tridiagonal or complex Hermitian tridiagonal?
Julien Langou
 
Posts: 735
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Which subroutine does complex symmetric tridiagonal matr

Postby carolasu » Wed Aug 15, 2012 12:49 am

Thanks. I mean complex Hermitian tridiagonal.
carolasu
 
Posts: 9
Joined: Wed Mar 21, 2012 12:50 am

Re: Which subroutine does complex symmetric tridiagonal matr

Postby Julien Langou » Wed Aug 15, 2012 1:14 am

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.
Last edited by Julien Langou on Wed Aug 15, 2012 3:05 am, edited 1 time in total.
Julien Langou
 
Posts: 735
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Which subroutine does complex Hermitian tridiagonal matr

Postby carolasu » Wed Aug 15, 2012 1:56 am

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!
Last edited by carolasu on Wed Aug 15, 2012 3:51 am, edited 1 time in total.
carolasu
 
Posts: 9
Joined: Wed Mar 21, 2012 12:50 am

Re: Which subroutine does complex Hermitian tridiagonal matr

Postby mattcimd » Mon Nov 12, 2012 3:26 am

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!!!
mattcimd
 
Posts: 1
Joined: Sun Nov 11, 2012 5:49 pm


Return to User Discussion

Who is online

Users browsing this forum: Bing [Bot] and 2 guests

cron