The subroutine eigen determines all eigenvalues and (if desired) all eigenvectors of a real square n * n matrix via the QR method in the version of Martin, Parlett, Peters, Reinsch and Wilkinson.
n : integer; ( n > 0 ) size of matrix, number of eigenvalues
mat : n x n matrix; input matrix
eivec : n x n matrix; (only if vec = 1) matrix, if vec = 1 that holds the eigenvectors thus : If the jth eigenvalue of the matrix is real then the jth column is the corresponding real eigenvector; if the jth eigenvalue is complex then the jth column of eivec contains the real part of the eigenvector while its imaginary part is in column j+1. (the j+1st eigenvector is the complex conjugate vector.)
valre : vector of size n; Real parts of the eigenvalues.
valim : vector of size n; Imaginary parts of the eigenvalues
cnt : Integer vector of size n; vector containing the number of iterations for each eigenvalue. (for a complex conjugate pair the second entry is negative).
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | n | |||
real(kind=double), | intent(in), | dimension(n,n) | :: | matrix | ||
real(kind=double), | intent(out), | dimension(n,n) | :: | eigvec | ||
real(kind=double), | intent(out), | dimension(n) | :: | eigval |
subroutine eigen (n, matrix, eigvec, eigval)
!+The subroutine eigen determines all eigenvalues and (if desired)
!+all eigenvectors of a real square n * n matrix via the QR method
!+in the version of Martin, Parlett, Peters, Reinsch and Wilkinson.
!+
!+###Literature
!+1. Peters, Wilkinson: Eigenvectors of real and complex
!+ matrices by LR and QR triangularisations,
!+ Num. Math. 16, p.184-204, (1970); [PETE70]; contribution
!+ II/15, p. 372 - 395 in [WILK71].
!+2. Martin, Wilkinson: Similarity reductions of a general
!+ matrix to Hessenberg form, Num. Math. 12, p. 349-368,(1968)
!+ [MART 68]; contribution II,13, p. 339 - 358 in [WILK71].
!+3. Parlett, Reinsch: Balancing a matrix for calculations of
!+ eigenvalues and eigenvectors, Num. Math. 13, p. 293-304,
!+ (1969); [PARL69]; contribution II/11, p.315 - 326 in
!+ [WILK71].
!+
!+###Input parameters
!+ n : integer; ( n > 0 )
!+ size of matrix, number of eigenvalues
!+
!+ mat : n x n matrix;
!+ input matrix
!+
!+###Output parameters
!+ eivec : n x n matrix; (only if vec = 1)
!+ matrix, if vec = 1 that holds the eigenvectors
!+ thus :
!+ If the jth eigenvalue of the matrix is real then the
!+ jth column is the corresponding real eigenvector;
!+ if the jth eigenvalue is complex then the jth column
!+ of eivec contains the real part of the eigenvector
!+ while its imaginary part is in column j+1.
!+ (the j+1st eigenvector is the complex conjugate
!+ vector.)
!+
!+ valre : vector of size n;
!+ Real parts of the eigenvalues.
!+
!+ valim : vector of size n;
!+ Imaginary parts of the eigenvalues
!+
!+ cnt : Integer vector of size n;
!+ vector containing the number of iterations for each
!+ eigenvalue. (for a complex conjugate pair the second
!+ entry is negative).
integer ,intent(in) :: n ! nlevels
real(double) ,intent(in),dimension(n,n) :: matrix
real(double) ,intent(out),dimension(n,n) :: eigvec
real(double) ,intent(out),dimension(n) :: eigval
real(double) :: mat(0:n-1,0:n-1)
real(double) :: eivec(0:n-1,0:n-1)
real(double) :: valre(0:n-1) !real parts of eigenvalues
real(double) :: valim(0:n-1) !imaginary parts of eigenvalues
integer :: rc !return code
integer :: cnt(0:n-1) !Iteration counter
integer :: high, low
real(double) :: d(0:n-1), scale(0:n-1)
integer :: perm(0:n-1)
cnt=0 ; d=0.d0
mat(0:n-1,0:n-1)=matrix(1:n,1:n)
!balance mat for nearly
call balance(n, mat, scale, low, high) !equal row and column
!reduce mat to upper
call elmhes(n, low, high, mat, perm) !reduce mat to upper
!Hessenberg form
call elmtrans(n, low, high, mat, perm, eivec)
!QR algorithm for eigenvalues and eigenvectors
call hqr2(n, low, high, mat, valre, valim, eivec, cnt,rc)
!reverse balancing to determine eigenvectors
call balback(n, low, high, scale, eivec)
if (rc.ne.0) then
print*, 'matrix = '
print*, matrix
stop 'problem in eigen!'
endif
eigval(1:n)=valre(0:n-1)
eigvec(1:n,1:n)=eivec(0:n-1,0:n-1)
end subroutine eigen