EIG

 

   more...

Computes the eigenvalue-eigenvector decomposition of an ordinary or generalized eigenvalue problem.

Function Return Value

Rank-1 or rank-2 real or complex array of eigenvalues. (Output)

Required Argument

A — Matrix for which the eigenexpansion is to be computed. This is a square rank-2 array or a rank-3 array with square first rank-2 sections of type single, double, complex, or double complex. (Input)

Optional Arguments, Packaged Options

B — Matrix B for the generalized problem, Ax = eBx. B must be the same type as A. (Input)

D — Array containing the real diagonal matrix factors of the generalized eigenvalues. (Output)

V — Array of real eigenvectors for both the ordinary and generalized problem. Used only for the generalized problem when matrix B is singular. (Output)

W — Array of complex eigenvectors for both the ordinary and generalized problem. Do not use optional argument V when W is present. (Output)

This function uses LIN_EIG_SELF, LIN_EIG_GEN, and LIN_GEIG_GEN to compute the decompositions. See Chapter 2, “Eigensystem Analysis”.

The option and derived type names are given in the following tables:

Option Names for EIG	Option Value
Options_for_lin_eig_self	1
Options_for_lin_eig_gen	2
Options_for_lin_geig_gen	3
Skip_error_processing	5

 

Name of Unallocated Option Array to Use for Setting Options	Use	Derived Type
?_eig_options(:)	Use when setting options for calls hereafter.	?_options
?_eig_options_once(:)	Use when setting options for next call only.	?_options

For a description on how to use these options, see “Matrix Optional Data Changes”. See LIN_EIG_SELF, LIN_EIG_GEN, and LIN_GEIG_GEN located in Chapter 2, “Eigensystem Analysis” for the specific options for these routines.

FORTRAN 90 Interface

EIG (A [, …] )

Description

Computes the eigenvalue-eigenvector decomposition of an ordinary or generalized eigenvalue problem.

For the ordinary eigenvalue problem Ax = ex, optional argument B is not used. With the generalized problem Ax = eBx, the optional argument B is used to input the matrix B. Optional output argument D is an array required only for the generalized problem and then only when the matrix B is singular.

The array of real eigenvectors is an optional output (V) for both the ordinary and the generalized problem. If any eigenvectors are complex, optional output W must be present. In that case V is not used.

Examples

Dense Matrix Example 1 (operator_ex26.f90)

 

use linear_operators

 

implicit none

 

! This is the equivalent of Example 2 (using operators) for LIN_EIG_SELF.

 

integer, parameter :: n=8

real(kind(1e0)), parameter :: one=1e0

real(kind(1e0)), dimension(n,n) :: A, d(n), v_s

 

! Generate a random self-adjoint matrix.

A = rand(A); A = A + .t.A

 

! Compute the eigenvalues and eigenvectors.

D = EIG(A,V=v_s)

 

! Check the results for small residuals.

if (norm((A .x. v_s) - (v_s .x. diag(D)))/abs(d(1)) <= &

sqrt(epsilon(one))) then

write (*,*) 'Example 2 for LIN_EIG_SELF (operators) is correct.'

end if

 

end

Dense Matrix Example 2 (operator_ex33.f90)

 

use linear_operators

 

implicit none

 

! This is the equivalent of Example 1 (using operators) for LIN_GEIG_GEN.

 

integer, parameter :: n=32

real(kind(1d0)), parameter :: one=1d0

real(kind(1d0)) A(n,n), B(n,n), bta(n), beta_t(n), err

complex(kind(1d0)) alpha(n), alpha_t(n), V(n,n)

 

! Generate random matrices for both A and B.

A = rand(A); B = rand(B)

 

! Compute the generalized eigenvalues.

alpha = EIG(A, B=B, D=bta)

 

! Compute the full decomposition once again, A*V = B*V*values,

! and check for any error messages.

alpha_t = EIG(A, B=B, D=beta_t, W = V)

 

! Use values from the first decomposition, vectors from the

! second decomposition, and check for small residuals.

err = norm((A .x. V .x. diag(bta)) - (B .x. V .x. diag(alpha)),1)/&

(norm(A,1)*norm(bta,1) + norm(B,1)*norm(alpha,1))

if (err <= sqrt(epsilon(one))) then

write (*,*) 'Example 1 for LIN_GEIG_GEN (operators) is correct.'

end if

 

end

Parallel Example (parallel_ex04.f90)

Here an alternate node is used to compute the majority of a single application, and the user does not need to make any explicit calls to MPI functions. The time-consuming parts are the evaluation of the eigenvalue-eigenvector expansion, the solving step, and the residuals. To do this, the rank-2 arrays are changed to a box data type with a unit third dimension. This uses parallel computing. The node priority order is established by the initial function call, MP_SETUP(n). The root is restricted from working on the box data type by assigning MPI_ROOT_WORKS=.false.. This example anticipates that the most efficient node, other than the root, will perform the heavy computing. Two nodes are required to execute.

 

use linear_operators

use mpi_setup_int

 

implicit none

 

! This is the equivalent of Parallel Example 4 for matrix exponential.

! The box dimension has a single rack.

integer, parameter :: n=32, k=128, nr=1

integer i

real(kind(1e0)), parameter :: one=1e0, t_max=one, delta_t=t_max/(k-1)

real(kind(1e0)) err(nr), sizes(nr), A(n,n,nr)

real(kind(1e0)) t(k), y(n,k,nr), y_prime(n,k,nr)

complex(kind(1e0)), dimension(n,nr) :: x(n,n,nr), z_0, &

Z_1(n,nr,nr), y_0, d

 

! Setup for MPI. Establish a node priority order.

! Restrict the root from significant computing.

! Illustrates using the 'best' performing node that

! is not the root for a single task.

MP_NPROCS=MP_SETUP(n)

 

MPI_ROOT_WORKS=.false.

 

! Generate a random coefficient matrix.

A = rand(A)

 

! Compute the eigenvalue-eigenvector decomposition

! of the system coefficient matrix on an alternate node.

D = EIG(A, W=X)

 

! Generate a random initial value for the ODE system.

y_0 = rand(y_0)

 

! Solve complex data system that transforms the initial

! values, X z_0=y_0.

 

z_1= X .ix. y_0 ; z_0(:,nr) = z_1(:,nr,nr)

 

! The grid of points where a solution is computed:

t = (/(i*delta_t,i=0,k-1)/)

 

 

! Compute y and y' at the values t(1:k).

! With the eigenvalue-eigenvector decomposition AX = XD, this

! is an evaluation of EXP(A t)y_0 = y(t).

y = X .x.exp(spread(d(:,nr),2,k)*spread(t,1,n))*spread(z_0(:,nr),2,k)

 

! This is y', derived by differentiating y(t).

y_prime = X .x. &

spread(d(:,nr),2,k)*exp(spread(d(:,nr),2,k)*spread(t,1,n))* &

spread(z_0(:,nr),2,k)

 

! Check results. Is y' - Ay = 0?

err = norm(y_prime-(A .x. y))

sizes=norm(y_prime)+norm(A)*norm(y)

if (ALL(err <= sqrt(epsilon(one))*sizes) .and. MP_RANK == 0) &

write (*,*) 'Parallel Example 4 is correct.'

 

! See to any error messages and quit MPI.

MP_NPROCS=MP_SETUP('Final')

 

end