LIN_EIG_GEN
Computes the eigenvalues of an n × n matrix, A. Optionally, the eigenvectors of A or AT are computed. Using the eigenvectors of A gives the decomposition AV = VE, where V is an n × n complex matrix of eigenvectors, and E is the complex diagonal matrix of eigenvalues. Other options include the reduction of A to upper triangular or Schur form, reduction to block upper triangular form with 2 × 2 or unit sized diagonal block matrices, and reduction to upper Hessenberg form.
Required Arguments
A — Array of size n × n containing the matrix. (Input [/Output])
E — Array of size n containing the eigenvalues. These complex values are in order of decreasing absolute value. The signs of imaginary parts of the eigenvalues are in no predictable order. (Output)
Optional Arguments
NROWS = n (Input)
Uses array A(1:n, 1:n) for the input matrix.
Default: n = SIZE(A, 1)
v = V(:,:) (Output)
Returns the complex array of eigenvectors for the matrix A.
v_adj = U(:,:) (Output)
Returns the complex array of eigenvectors for the matrix AT. Thus the residuals
are small.
tri = T(:,:) (Output)
Returns the complex upper-triangular matrix T associated with the reduction of the matrix A to Schur form. Optionally a unitary matrix W is returned in array V(:,:) such that the residuals Z = AW – WT are small.
iopt = iopt(:) (Input)
Derived type array with the same precision as the input matrix. Used for passing optional data to the routine. The options are as follows:
|
Packaged Options for LIN_EIG_GEN |
||
|
Option Prefix = ? |
Option Name |
Option Value |
|
s_, d_, c_, z_ |
lin_eig_gen_set_small |
1 |
|
s_, d_, c_, z_ |
lin_eig_gen_overwrite_input |
2 |
|
s_, d_, c_, z_ |
lin_eig_gen_scan_for_NaN |
3 |
|
s_, d_, c_, z_ |
lin_eig_gen_no_balance |
4 |
|
s_, d_, c_, z_ |
lin_eig_gen_set_iterations |
5 |
|
s_, d_, c_, z_ |
lin_eig_gen_in_Hess_form |
6 |
|
s_, d_, c_, z_ |
lin_eig_gen_out_Hess_form |
7 |
|
s_, d_, c_, z_ |
lin_eig_gen_out_block_form |
8 |
|
s_, d_, c_, z_ |
lin_eig_gen_out_tri_form |
9 |
|
s_, d_, c_, z_ |
lin_eig_gen_continue_with_V |
10 |
|
s_, d_, c_, z_ |
lin_eig_gen_no_sorting |
11 |
iopt(IO) = ?_options(?_lin_eig_gen_set_small, Small)
This is the tolerance used to declare off-diagonal values effectively zero compared with the size of the numbers involved in the computation of a shift.
Default: Small = epsilon(), the relative accuracy of arithmetic
iopt(IO) = ?_options(?_lin_eig_gen_overwrite_input, ?_dummy)
Does not save the input array A(:, :).
Default: The array is saved.
iopt(IO) = ?_options(?_lin_eig_gen_scan_for_NaN, ?_dummy)
Examines each input array entry to find the first value such that
isNaN(a(i,j)) == .true..
See the isNaN() function, Chapter 10.
Default: The array is not scanned for NaNs.
iopt(IO) = ?_options(?_lin_eig_no_balance, ?_dummy)
The input matrix is not preprocessed searching for isolated eigenvalues followed by rescaling. See Golub and Van Loan (1989, Chapter 7) for references. With some optional uses of the routine, this option flag is required.
Default: The matrix is first balanced.
iopt(IO) = ?_options(?_lin_eig_gen_set_iterations, ?_dummy)
Resets the maximum number of iterations permitted to isolate each diagonal block matrix.
Default: The maximum number of iterations is 52.
iopt(IO) = ?_options(?_lin_eig_gen_in_Hess_form, ?_dummy)
The input matrix is in upper Hessenberg form. This flag is used to avoid the initial reduction phase which may not be needed for some problem classes.
Default: The matrix is first reduced to Hessenberg form.
iopt(IO) = ?_options(?_lin_eig_gen_out_Hess_form, ?_dummy)
The output matrix is transformed to upper Hessenberg form, 
. If the optional argument “v=V(:,:)” is passed by the calling program unit, then the array V(:,:) contains an orthogonal matrix Q1 such that
AQ1 - Q1H1 ≅ 0
Requires the simultaneous use of option ?_lin_eig_no_balance.
Default: The matrix is reduced to diagonal form.
iopt(IO) = ?_options(?_lin_eig_gen_out_block_form, ?_dummy)
The output matrix is transformed to upper Hessenberg form, H2, which is block upper triangular. The dimensions of the blocks are either 2 × 2 or unit sized. Nonzero subdiagonal values of H2 determine the size of the blocks. If the optional argument “v=V(:,:)” is passed by the calling program unit, then the array V(:,:) contains an orthogonal matrix Q2, such that
AQ2 - Q2H2 ≅ 0
Requires the simultaneous use of option ?_lin_eig_no_balance.
Default: The matrix is reduced to diagonal form.
iopt(IO) = ?_options(?_lin_eig_gen_out_tri_form, ?_dummy)
The output matrix is transformed to upper-triangular form, T. If the optional argument “v=V(:,:)” is passed by the calling program unit, then the array V(:,:) contains a unitary matrix W such that AW – WT ≅ 0. The upper triangular matrix T is returned in the optional argument “tri=T(:,:)”. The eigenvalues of A are the diagonal entries of the matrix T. They are in no particular order. The output array E(:)is blocked with NaNs using this option. This option requires the simultaneous use of option ?_lin_eig_no_balance.
Default: The matrix is reduced to diagonal form.
iopt(IO) = ?_options(?_lin_eig_gen_continue_with_V, ?_dummy)
As a convenience or for maintaining efficiency, the calling program unit sets the optional argument “v=V(:,:)” to a matrix that has transformed a problem to the similar matrix, 
. The contents of V(:,:) are updated by the transformations used in the algorithm. Requires the simultaneous use of option ?_lin_eig_no_balance.
Default: The array V(:,:) is initialized to the identity matrix.
iopt(IO) = ?_options(?_lin_eig_gen_no_sorting, ?_dummy)
Does not sort the eigenvalues as they are isolated by solving the 2 × 2 or unit sized blocks. This will have the effect of guaranteeing that complex conjugate pairs of eigenvalues are adjacent in the array E(:).
Default: The entries of E(:) are sorted so they are non-increasing in absolute value.
FORTRAN 90 Interface
Generic: CALL LIN_EIG_GEN (A, E [,…])
Specific: The specific interface names are S_LIN_EIG_GEN, D_LIN_EIG_GEN, C_LIN_EIG_GEN, and Z_LIN_EIG_GEN.
Description
The input matrix A is first balanced. The resulting similar matrix is transformed to upper Hessenberg form using orthogonal transformations. The double-shifted QR algorithm transforms the Hessenberg matrix so that 2 × 2 or unit sized blocks remain along the main diagonal. Any off-diagonal that is classified as “small” in order to achieve this block form is set to the value zero. Next the block upper triangular matrix is transformed to upper triangular form with unitary rotations. The eigenvectors of the upper triangular matrix are computed using back substitution. Care is taken to avoid overflows during this process. At the end, eigenvectors are normalized to have Euclidean length one, with the largest component real and positive. This algorithm follows that given in Golub and Van Loan, (1989, Chapter 7), with some novel organizational details for additional options, efficiency and robustness.
Fatal, Terminal, and Warning Error Messages
See the messages.gls file for error messages for LIN_EIG_GEN. These error messages are numbered 841–858; 861–878; 881–898; 901–918.
Examples
Example 1: Computing Eigenvalues
The eigenvalues of a random real matrix are computed. These values define a complex diagonal matrix E. Their correctness is checked by obtaining the eigenvector matrix V and verifying that the residuals R = AV ‑ VE are small. Also, see operator_ex29, supplied with the product examples.
use lin_eig_gen_int
use rand_gen_int
implicit none
! This is Example 1 for LIN_EIG_GEN.
integer, parameter :: n=32
real(kind(1d0)), parameter :: one=1d0
real(kind(1d0)) A(n,n), y(n*n), err
complex(kind(1d0)) E(n), V(n,n), E_T(n)
type(d_error) :: d_epack(16) = d_error(0,0d0)
! Generate a random matrix.
call rand_gen(y)
A = reshape(y,(/n,n/))
! Compute only the eigenvalues.
call lin_eig_gen(A, E)
! Compute the decomposition, A*V = V*values,
! obtaining eigenvectors.
call lin_eig_gen(A, E_T, v=V)
! Use values from the first decomposition, vectors from the
! second decomposition, and check for small residuals.
err = sum(abs(matmul(A,V) - V*spread(E,DIM=1,NCOPIES=n))) &
/ sum(abs(E))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 1 for LIN_EIG_GEN is correct.'
end if
end
Example 1 for LIN_EIG_GEN is correct.
Example 2: Complex Polynomial Equation Roots
The roots of a complex polynomial equation,
are required. This algebraic equation is formulated as a matrix eigenvalue problem. The equivalent matrix eigenvalue problem is solved using the upper Hessenberg matrix which has the value zero except in row number 1 and along the first subdiagonal. The entries in the first row are given by a1,j = – bj, i = 1, …, n, while those on the first subdiagonal have the value one. This is a companion matrix for the polynomial. The results are checked by testing for small values of ∣f(ei)∣, i = 1, …, n, at the eigenvalues of the matrix, which are the roots of f(z). Also, see operator_ex30, supplied with the product examples.
use lin_eig_gen_int
use rand_gen_int
implicit none
! This is Example 2 for LIN_EIG_GEN.
integer i
integer, parameter :: n=12
real(kind(1d0)), parameter :: one=1.0d0, zero=0.0d0
real(kind(1d0)) err, t(2*n)
type(d_options) :: iopti(1)=d_options(0,zero)
complex(kind(1d0)) a(n,n), b(n), e(n), f(n), fg(n)
call rand_gen(t)
b = cmplx(t(1:n),t(n+1:),kind(one))
! Define the companion matrix with polynomial coefficients
! in the first row.
a = zero
do i=2, n
a(i,i-1) = one
end do
a(1,1:n) = -b
! Note that the input companion matrix is upper Hessenberg.
iopti(1) = d_options(z_lin_eig_gen_in_Hess_form,zero)
! Compute complex eigenvalues of the companion matrix.
call lin_eig_gen(a, e, iopt=iopti)
f=one; fg=one
! Use Horner's method for evaluation of the complex polynomial
! and size gauge at all roots.
do i=1, n
f = f*e + b(i)
fg = fg*abs(e) + abs(b(i))
end do
! Check for small errors at all roots.
err = sum(abs(f/fg))/n
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 2 for LIN_EIG_GEN is correct.'
end if
end
Output
Example 2 for LIN_EIG_GEN is correct.
Example 3: Solving Parametric Linear Systems with a Scalar Change
The efficient solution of a family of linear algebraic equations is required. These systems are (A + hI)x = b. Here A is an n × n real matrix, I is the identity matrix, and b is the right-hand side matrix. The scalar h is such that the coefficient matrix is nonsingular. The method is based on the Schur form for matrix A: AW = WT, where W is unitary and T is upper triangular. This provides an efficient solution method for several values of h, once the Schur form is computed. The solution steps solve, for y, the upper triangular linear system
Then, x = x(h) = Wy. This is an efficient and accurate method for such parametric systems provided the expense of computing the Schur form has a pay-off in later efficiency. Using the Schur form in this way, it is not required to compute an LU factorization of A + hI with each new value of h. Note that even if the data A, h, and b are real, subexpressions for the solution may involve complex intermediate values, with x(h) finally a real quantity. Also, see operator_ex31, supplied with the product examples.
use lin_eig_gen_int
use lin_sol_gen_int
use rand_gen_int
implicit none
! This is Example 3 for LIN_EIG_GEN.
integer i
integer, parameter :: n=32, k=2
real(kind(1e0)), parameter :: one=1.0e0, zero=0.0e0
real(kind(1e0)) a(n,n), b(n,k), x(n,k), temp(n*max(n,k)), h, err
type(s_options) :: iopti(2)
complex(kind(1e0)) w(n,n), t(n,n), e(n), z(n,k)
call rand_gen(temp)
a = reshape(temp,(/n,n/))
call rand_gen(temp)
b = reshape(temp,(/n,k/))
iopti(1) = s_options(s_lin_eig_gen_out_tri_form,zero)
iopti(2) = s_options(s_lin_eig_gen_no_balance,zero)
! Compute the Schur decomposition of the matrix.
call lin_eig_gen(a, e, v=w, tri=t, &
iopt=iopti)
! Choose a value so that A+h*I is non-singular.
h = one
! Solve for (A+h*I)x=b using the Schur decomposition.
z = matmul(conjg(transpose(w)),b)
! Solve intermediate upper-triangular system with implicit
! additive diagonal, h*I. This is the only dependence on
! h in the solution process.
do i=n,1,-1
z(i,1:k) = z(i,1:k)/(t(i,i)+h)
z(1:i-1,1:k) = z(1:i-1,1:k) + &
spread(-t(1:i-1,i),dim=2,ncopies=k)* &
spread(z(i,1:k),dim=1,ncopies=i-1)
end do
! Compute the solution. It should be the same as x, but will not be
! exact due to rounding errors. (The quantity real(z,kind(one)) is
! the real-valued answer when the Schur decomposition method is used.)
z = matmul(w,z)
! Compute the solution by solving for x directly.
do i=1, n
a(i,i) = a(i,i) + h
end do
call lin_sol_gen(a, b, x)
! Check that x and z agree approximately.
err = sum(abs(x-z))/sum(abs(x))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 3 for LIN_EIG_GEN is correct.'
end if
end
Output
Example 3 for LIN_EIG_GEN is correct.
Example 4: Accuracy Estimates of Eigenvalues Using Adjoint and Ordinary Eigenvectors
A matrix A has entries that are subject to uncertainty. This is expressed as the realization that A can be replaced by the matrix A + ηB, where the value η is “small” but still significantly larger than machine precision. The matrix B satisfies ∥B∥ ≤ ∥A∥. A variation in eigenvalues is estimated using analysis found in Golub and Van Loan, (1989, Chapter 7, p. 344). Each eigenvalue and eigenvector is expanded in a power series in η. With
and normalized eigenvectors, the bound
is satisfied. The vectors ui and vi are the ordinary and adjoint eigenvectors associated respectively with ei and its complex conjugate. This gives an upper bound on the size of the change to each ∣ei∣ due to changing the matrix data. The reciprocal
is defined as the condition number of ei. Also, see operator_ex32, Chapter 10.
use lin_eig_gen_int
use rand_gen_int
implicit none
! This is Example 4 for LIN_EIG_GEN.
integer i
integer, parameter :: n=17
real(kind(1d0)), parameter :: one=1d0
real(kind(1d0)) a(n,n), c(n,n), variation(n), y(n*n), temp(n), &
norm_of_a, eta
complex(kind(1d0)), dimension(n,n) :: e(n), d(n), u, v
! Generate a random matrix.
call rand_gen(y)
a = reshape(y,(/n,n/))
! Compute the eigenvalues, left- and right- eigenvectors.
call lin_eig_gen(a, e, v=v, v_adj=u)
! Compute condition numbers and variations of eigenvalues.
norm_of_a = sqrt(sum(a**2)/n)
do i=1, n
variation(i) = norm_of_a/abs(dot_product(u(1:n,i), &
v(1:n,i)))
end do
! Now perturb the data in the matrix by the relative factors
! eta=sqrt(epsilon) and solve for values again. Check the
! differences compared to the estimates. They should not exceed
! the bounds.
eta = sqrt(epsilon(one))
do i=1, n
call rand_gen(temp)
c(1:n,i) = a(1:n,i) + (2*temp - 1)*eta*a(1:n,i)
end do
call lin_eig_gen(c,d)
! Looking at the differences of absolute values accounts for
! switching signs on the imaginary parts.
if (count(abs(d)-abs(e) > eta*variation) == 0) then
write (*,*) 'Example 4 for LIN_EIG_GEN is correct.'
end if
end
Output
Example 4 for LIN_EIG_GEN is correct.
