LIN_GEIG_GEN
Computes the generalized eigenvalues of an n × n matrix pencil, Av = λBv. Optionally, the generalized eigenvectors are computed. If either of A or B is nonsingular, there are diagonal matrices α and β, and a complex matrix V, all computed such that AV β = BVα.
Required Arguments
A — Array of size n × n containing the matrix A. (Input [/Output])
B — Array of size n × n containing the matrix B. (Input [/Output])
ALPHA — Array of size n containing diagonal matrix factors of the generalized eigenvalues. These complex values are in order of decreasing absolute value. (Output)
BETAV — Array of size n containing diagonal matrix factors of the generalized eigenvalues. These real values are in order of decreasing value. (Output)
Optional Arguments
NROWS = n (Input)
Uses arrays A(1:n, 1:n) and B(1:n, 1:n) for the input matrix pencil.
Default: n = SIZE(A, 1)
Uses arrays A(1:n, 1:n) and B(1:n, 1:n) for the input matrix pencil.
Default: n = SIZE(A, 1)
v = V(:,:) (Output)
Returns the complex array of generalized eigenvectors for the matrix pencil.
Returns the complex array of generalized eigenvectors for the matrix pencil.
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:
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_GEIG_GEN | ||
|---|---|---|
Option Prefix = ? | Option Name | Option Value |
s_, d_, c_, z_ | lin_geig_gen_set_small | 1 |
s_, d_, c_, z_ | lin_geig_gen_overwrite_input | 2 |
s_, d_, c_, z_ | lin_geig_gen_scan_for_NaN | 3 |
s_, d_, c_, z_ | lin_geig_gen_self_adj_pos | 4 |
s_, d_, c_, z_ | lin_geig_gen_for_lin_sol_self | 5 |
s_, d_, c_, z_ | lin_geig_gen_for_lin_eig_self | 6 |
s_, d_, c_, z_ | lin_geig_gen_for_lin_sol_lsq | 7 |
s_, d_, c_, z_ | lin_geig_gen_for_lin_eig_gen | 8 |
iopt(IO) = ?_options(?_lin_geig_gen_set_small, Small)
This tolerance, multiplied by the sum of absolute value of the matrix B, is used to define a small diagonal term in the routines lin_sol_lsq and lin_sol_self. That value can be replaced using the option flags lin_geig_gen_for_lin_sol_lsq, and lin_geig_gen_for_lin_sol_self.
Default: Small = epsilon(.), the relative accuracy of arithmetic.
This tolerance, multiplied by the sum of absolute value of the matrix B, is used to define a small diagonal term in the routines lin_sol_lsq and lin_sol_self. That value can be replaced using the option flags lin_geig_gen_for_lin_sol_lsq, and lin_geig_gen_for_lin_sol_self.
Default: Small = epsilon(.), the relative accuracy of arithmetic.
iopt(IO) = ?_options(?_lin_geig_gen_overwrite_input, ?_dummy)
Does not save the input arrays A(:, :) and B(:, :).
Default: The array is saved.
Does not save the input arrays A(:, :) and B(:, :).
Default: The array is saved.
iopt(IO) = ?_options(?_lin_geig_gen_scan_for_NaN, ?_dummy)
Examines each input array entry to find the first value such that
Examines each input array entry to find the first value such that
isNaN(a(i,j)) .or. isNaN(b(i,j)) == .true.
See the isNaN() function, Chapter 10.
Default: The arrays are not scanned for NaNs.
Default: The arrays are not scanned for NaNs.
iopt(IO) = ?_options(?_lin_geig_gen_self_adj_pos, ?_dummy)
If both matrices A and B are self-adjoint and additionally B is positive-definite, then the Cholesky algorithm is used to reduce the matrix pencil to an ordinary self-adjoint eigenvalue problem.
If both matrices A and B are self-adjoint and additionally B is positive-definite, then the Cholesky algorithm is used to reduce the matrix pencil to an ordinary self-adjoint eigenvalue problem.
iopt(IO) = ?_options(?_lin_geig_gen_for_lin_sol_self, ?_dummy)
iopt(IO+1) = ?_options((k=size of options for lin_sol_self), ?_dummy)
The options for lin_sol_self follow as data in iopt().
The options for lin_sol_self follow as data in iopt().
iopt(IO) = ?_options(?_lin_geig_gen_for_lin_eig_self, ?_dummy)
iopt(IO+1) = ?_options((k=size of options for lin_eig_self), ?_dummy)
The options for lin_eig_self follow as data in iopt().
The options for lin_eig_self follow as data in iopt().
iopt(IO) = ?_options(?_lin_geig_gen_for_lin_sol_lsq, ?_dummy)
iopt(IO+1) = ?_options((k=size of options for lin_sol_lsq), ?_dummy)
The options for lin_sol_lsq follow as data in iopt().
The options for lin_sol_lsq follow as data in iopt().
iopt(IO) = ?_options(?_lin_geig_gen_for_lin_eig_gen, ?_dummy)
iopt(IO+1) = ?_options((k=size of options for lin_eig_gen), ?_dummy)
The options for lin_eig_gen follow as data in iopt().
The options for lin_eig_gen follow as data in iopt().
FORTRAN 90 Interface
Generic: CALL LIN_GEIG_GEN (A, B, ALPHA, BETAV [,…])
Specific: The specific interface names are S_LIN_GEIG_GEN, D_LIN_GEIG_GEN, C_LIN_GEIG_GEN, and Z_LIN_GEIG_GEN.
Description
Routine LIN_GEIG_GEN implements a standard algorithm that reduces a generalized eigenvalue or matrix pencil problem to an ordinary eigenvalue problem. An orthogonal decomposition is computed
BPT = HR
The orthogonal matrix H is the product of n – 1 row permutations, each followed by a Householder transformation. Column permutations, P, are chosen at each step to maximize the Euclidian length of the pivot column. The matrix R is upper triangular. Using the default tolerance = ɛ∥B∥, where ɛ is machine relative precision, each diagonal entry of R exceeds in value. Otherwise, R is singular. In that case A and B are interchanged and the orthogonal decomposition is computed one more time. If both matrices are singular the problem is declared singular and is not solved. The interchange of A and B is accounted for in the output diagonal matrices α and β. The ordinary eigenvalue problem is Cx = λx, where
C =HT APT R-1
and
RPv = x
If the matrices A and B are self-adjoint and if, in addition, B is positive-definite, then a more efficient reduction than the default algorithm can be optionally used to solve the problem: A Cholesky decomposition is obtained, RTR R = PBPT. The matrix R is upper triangular and P is a permutation matrix. This is equivalent to the ordinary self-adjoint eigenvalue problem Cx = λx, where RPv = x and
C =R-T PAPT R-1
The self-adjoint eigenvalue problem is then solved.
Fatal, Terminal, and Warning Error Messages
See the messages.gls file for error messages for LIN_GEIG_GEN. These error messages are numbered 921–936; 941–956; 961–976; 981–996.
Examples
Example 1: Computing Generalized Eigenvalues
The generalized eigenvalues of a random real matrix pencil are computed. These values are checked by obtaining the generalized eigenvectors and then showing that the residuals
AV - BV αβ-1
are small. Note that when the matrix B is nonsingular β = I, the identity matrix. When B is singular and A is nonsingular, some diagonal entries of β are essentially zero. This corresponds to “infinite eigenvalues” of the matrix pencil. This random matrix pencil example has all finite eigenvalues. Also, see operator_ex33, Chapter 10.
use lin_geig_gen_int
use rand_gen_int
implicit none
! This is Example 1 for LIN_GEIG_GEN.
integer, parameter :: n=32
real(kind(1d0)), parameter :: one=1d0
real(kind(1d0)) A(n,n), B(n,n), betav(n), beta_t(n), err, y(n*n)
complex(kind(1d0)) alpha(n), alpha_t(n), V(n,n)
! Generate random matrices for both A and B.
call rand_gen(y)
A = reshape(y,(/n,n/))
call rand_gen(y)
B = reshape(y,(/n,n/))
! Compute the generalized eigenvalues.
call lin_geig_gen(A, B, alpha, betav)
! Compute the full decomposition once again, A*V = B*V*values.
call lin_geig_gen(A, B, alpha_t, beta_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) - &
matmul(B,V)*spread(alpha/betav,DIM=1,NCOPIES=n))) / &
sum(abs(a)+abs(b))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 1 for LIN_GEIG_GEN is correct.'
end if
end
Output
Example 1 for LIN_GEIG_GEN is correct.
Example 2: Self-Adjoint, Positive-Definite Generalized Eigenvalue Problem
This example illustrates the use of optional flags for the special case where A and B are complex self-adjoint matrices, and B is positive-definite. For purposes of maximum efficiency an option is passed to routine LIN_SOL_SELF so that pivoting is not used in the computation of the Cholesky decomposition of matrix B. This example does not require that secondary option. Also, see operator_ex34, supplied with the product examples.
use lin_geig_gen_int
use lin_sol_self_int
use rand_gen_int
implicit none
! This is Example 2 for LIN_GEIG_GEN.
integer i
integer, parameter :: n=32
real(kind(1d0)), parameter :: one=1.0d0, zero=0.0d0
real(kind(1d0)) betav(n), temp_c(n,n), temp_d(n,n), err
type(d_options) :: iopti(4)=d_options(0,zero)
complex(kind(1d0)), dimension(n,n) :: A, B, C, D, V, alpha(n)
! Generate random matrices for both A and B.
do i=1, n
call rand_gen(temp_c(1:n,i))
call rand_gen(temp_d(1:n,i))
end do
c = temp_c; d = temp_c
do i=1, n
call rand_gen(temp_c(1:n,i))
call rand_gen(temp_d(1:n,i))
end do
c = cmplx(real(c),temp_c,kind(one))
d = cmplx(real(d),temp_d,kind(one))
a = conjg(transpose(c)) + c
b = matmul(conjg(transpose(d)),d)
! Set option so that the generalized eigenvalue solver uses an
! efficient method for well-posed, self-adjoint problems.
iopti(1) = d_options(z_lin_geig_gen_self_adj_pos,zero)
iopti(2) = d_options(z_lin_geig_gen_for_lin_sol_self,zero)
! Number of secondary optional data items and the options:
iopti(3) = d_options(1,zero)
iopti(4) = d_options(z_lin_sol_self_no_pivoting,zero)
call lin_geig_gen(a, b, alpha, betav, v=v, &
iopt=iopti)
! Check that residuals are small. Use the real part of alpha
! since the values are known to be real.
err = sum(abs(matmul(a,v) - matmul(b,v)* &
spread(real(alpha,kind(one))/betav,dim=1,ncopies=n))) / &
sum(abs(a)+abs(b))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 2 for LIN_GEIG_GEN is correct.'
end if
end
Output
Example 2 for LIN_GEIG_GEN is correct.
Example 3: A Test for a Regular Matrix Pencil
In the classification of Differential Algebraic Equations (DAE), a system with linear constant coefficients is given by A + Bx = f. Here A and B are n × n matrices, and f is an n-vector that is not part of this example. The DAE system is defined as solvable if and only if the quantity det (μA + B) does not vanish identically as a function of the dummy parameter μ. A sufficient condition for solvability is that the generalized eigenvalue problem Av = λBv is nonsingular. By constructing A and B so that both are singular, the routine flags nonsolvability in the DAE by returning NaN for the generalized eigenvalues. Also, see operator_ex35, supplied with the product examples.
use lin_geig_gen_int
use rand_gen_int
use error_option_packet
use isnan_int
implicit none
! This is Example 3 for LIN_GEIG_GEN.
integer, parameter :: n=6
real(kind(1d0)), parameter :: one=1.0d0, zero=0.0d0
real(kind(1d0)) a(n,n), b(n,n), betav(n), y(n*n)
type(d_options) iopti(1)
type(d_error) epack(1)
complex(kind(1d0)) alpha(n)
! Generate random matrices for both A and B.
call rand_gen(y)
a = reshape(y,(/n,n/))
call rand_gen(y)
b = reshape(y,(/n,n/))
! Make columns of A and B zero, so both are singular.
a(1:n,n) = 0; b(1:n,n) = 0
! Set internal tolerance for a small diagonal term.
iopti(1) = d_options(d_lin_geig_gen_set_small,sqrt(epsilon(one)))
! Compute the generalized eigenvalues.
call lin_geig_gen(a, b, alpha, betav, &
iopt=iopti,epack=epack)
! See if singular DAE system is detected.
! (The size of epack() is too small for the message, so
! output is blocked with NaNs.)
if (isnan(alpha)) then
write (*,*) 'Example 3 for LIN_GEIG_GEN is correct.'
end if
end
Output
Example 3 for LIN_GEIG_GEN is correct.
Example 4: Larger Data Uncertainty than Working Precision
Data values in both matrices A and B are assumed to have relative errors that can be as large as Ɛ 1/2 where Ɛ is the relative machine precision. This example illustrates the use of an optional flag that resets the tolerance used in routine lin_sol_lsq for determining a singularity of either matrix. The tolerance is reset to the new value Ɛ 1/2∥B∥ and the generalized eigenvalue problem is solved. We anticipate that B might be singular and detect this fact. Also, see ---operator_ex36, Chapter 10.
use lin_geig_gen_int
use lin_sol_lsq_int
use rand_gen_int
use isNaN_int
implicit none
! This is Example 4 for LIN_GEIG_GEN.
integer, parameter :: n=32
real(kind(1d0)), parameter :: one=1d0, zero=0d0
real(kind(1d0)) a(n,n), b(n,n), betav(n), y(n*n), err
type(d_options) iopti(4)
type(d_error) epack(1)
complex(kind(1d0)) alpha(n), v(n,n)
! Generate random matrices for both A and B.
call rand_gen(y)
a = reshape(y,(/n,n/))
call rand_gen(y)
b = reshape(y,(/n,n/))
! Set the option, a larger tolerance than default for lin_sol_lsq.
iopti(1) = d_options(d_lin_geig_gen_for_lin_sol_lsq,zero)
! Number of secondary optional data items
iopti(2) = d_options(2,zero)
iopti(3) = d_options(d_lin_sol_lsq_set_small,sqrt(epsilon(one))*&
sqrt(sum(b**2)/n))
iopti(4) = d_options(d_lin_sol_lsq_no_sing_mess,zero)
! Compute the generalized eigenvalues.
call lin_geig_gen(A, B, alpha, betav, v=v, &
iopt=iopti, epack=epack)
if(.not. isNaN(alpha)) then
! Check the residuals.
err = sum(abs(matmul(A,V)*spread(betav,dim=1,ncopies=n) - &
matmul(B,V)*spread(alpha,dim=1,ncopies=n))) / &
sum(abs(a)+abs(b))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 4 for LIN_GEIG_GEN is correct.'
end if
end if
end
Output
Example 4 for LIN_GEIG_GEN is correct.
