eigen Subroutine

public 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).

Arguments

TypeIntentOptionalAttributesName
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

Calls

proc~~eigen~~CallsGraph proc~eigen eigen proc~elmhes elmhes proc~eigen->proc~elmhes proc~hqr2 hqr2 proc~eigen->proc~hqr2 proc~balance balance proc~eigen->proc~balance proc~balback balback proc~eigen->proc~balback proc~elmtrans elmtrans proc~eigen->proc~elmtrans proc~rswap RSWAP proc~elmhes->proc~rswap proc~hqrvec hqrvec proc~hqr2->proc~hqrvec proc~balance->proc~rswap proc~balback->proc~rswap proc~comdiv Comdiv proc~hqrvec->proc~comdiv

Called by

proc~~eigen~~CalledByGraph proc~eigen eigen proc~colrad colrad proc~colrad->proc~eigen proc~ndmc ndmc proc~ndmc->proc~colrad proc~fida_weights_los fida_weights_los proc~fida_weights_los->proc~colrad proc~dcx dcx proc~dcx->proc~colrad proc~fida_f fida_f proc~fida_f->proc~colrad proc~pfida_f pfida_f proc~pfida_f->proc~colrad proc~read_equilibrium read_equilibrium proc~read_equilibrium->proc~colrad proc~attenuate attenuate proc~attenuate->proc~colrad proc~halo halo proc~halo->proc~colrad proc~pfida_mc pfida_mc proc~pfida_mc->proc~colrad proc~fida_weights_mc fida_weights_mc proc~fida_weights_mc->proc~colrad proc~fida_mc fida_mc proc~fida_mc->proc~colrad program~fidasim fidasim program~fidasim->proc~ndmc program~fidasim->proc~fida_weights_los program~fidasim->proc~dcx program~fidasim->proc~fida_f program~fidasim->proc~pfida_f program~fidasim->proc~read_equilibrium program~fidasim->proc~halo program~fidasim->proc~pfida_mc program~fidasim->proc~fida_weights_mc program~fidasim->proc~fida_mc proc~pnpa_mc pnpa_mc program~fidasim->proc~pnpa_mc proc~pnpa_f pnpa_f program~fidasim->proc~pnpa_f proc~npa_mc npa_mc program~fidasim->proc~npa_mc proc~npa_f npa_f program~fidasim->proc~npa_f proc~npa_weights npa_weights program~fidasim->proc~npa_weights proc~pnpa_mc->proc~attenuate proc~pnpa_f->proc~attenuate proc~npa_mc->proc~attenuate proc~npa_f->proc~attenuate proc~npa_weights->proc~attenuate

Contents

Source Code


Source Code

  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