EIG
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 routines. 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
