Computes the generalized eigenvalues of an n × n matrix pencil, Av = lBv. Optionally, the
gen­eralized 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, Chapter 10.
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 com­puted

The orthogonal matrix H is the product of n - 1 row permutations, each fol­lowed by a Householder transformation. Column permutations, P, are chosen at each step to max­imize 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 account­ed for in the output diagonal matrices α and β. The ordinary eigenvalue problem is Cx = lx, where

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 = lx, where RPv = x and

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.

.p>.F90CH2.DOC!EX1_LIN_GEIG_GEN;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

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.

Additional Examples

.p>.F90CH2.DOC!EX2_LIN_GEIG_GEN;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.

.p>.F90CH2.DOC!EX3_LIN_GEIG_GEN;Example 3: A Test for a Regular Matrix Pencil

In the classification of Differential Algebraic Equations (DAE), a system with linear constant
coef­ficients 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 con­dition for solvability is that the generalized eigenvalue problem Av = lBv is nonsingular. By con­structing 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.

.p>.F90CH2.DOC!EX4_LIN_GEIG_GEN;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  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  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.

Visual Numerics, Inc. Visual Numerics - Developers of IMSL and PV-WAVE http://www.vni.com/ PHONE: 713.784.3131 FAX:713.781.9260