LIN_GEIG_GEN

 

more...

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)

v = V(:,:) (Output)
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:

 

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.

iopt(IO) = ?_options(?_lin_geig_gen_overwrite_input,  ?_dummy)
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

isNaN(a(i,j)) .or. isNaN(b(i,j)) == .true.

See the isNaN() function, Linear Algebra Operators and Generic Functions.
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.

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

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

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

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

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, Linear Algebra Operators and Generic Functions.

 

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, Linear Algebra Operators and Generic Functions.

 

 

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.