Compute some eigenvalues and eigenvectors of the generalized eigenvalue problem Ax = λBx. This can be used for the case B = I. The values for A,B are real, but eigenvalues may be complex and occur in conjugate pairs.
Required Arguments
N — The dimension of the problem. (Input)
F — User-supplied FUNCTION to return matrix-vector operations or linear solutions. This user function is written corresponding to the abstract interface for the function DMV(…). The usage is F (X, TASK, EXTYPE), where
Function Return Value
F — An array of length N containing matrix-vector operations or linear equations solutions. Operations provided as code in the function F will be made depending upon the value of argument TASK.
Required Arguments
X — An array of length N containing the vector to which the operator will be applied. (Input)
TASK — An enumerated type which specifies the operation to be performed. (Input) TASK is an enumerated integer value, use-associated from the module ARPACK_INT. It will be one of the following:
Value
Description
ARPACK_Prepare
Take initial steps to prepare for the operations that follow. These steps can include defining the data for the matrices, factorizations for upcoming linear system solves, or recording the vectors used in the operations.
ARPACK_A_x
y = Ax
ARPACK_B_x
y = Bx
ARPACK_inv_A_minus_Shift_x
y=(A-σ I)-1x
ARPACK_inv_B_x
y=B-1x
ARPACK_inv_A_minus_Shift_B_x
y=(A-σ B)-1x
EXTYPE — A derived type of the extensible class ARPACKBASE, which may be used to pass additional information to/from the user-supplied function. (Input/Output) The user must include a USEARPACK_INT statement in the calling program to define this derived type. If EXTYPE is not included as an argument to ARPACK_NONSYMMETRIC it should be ignored in the user-function, F.
The function F must be written according to the abstract interface for DMV. If F is not use-associated nor contained in the calling program, declare it with PROCEDURE(DMV) F.
ZVALUES — A complex array of eigenvalues. (Output) The value NEV=size(ZVALUES) defines the number of eigenvalues to be computed. The calling program declares or allocates the array ZVALUES(1:NEV). The size value NEV should account for pairs of complex conjugates. The number of eigenvalues computed accurately is optionally available as the component EXTYPE%NACC of the base class EXTYPE.
Optional Arguments
PLACE — Defines the placement in the spectrum for required eigenvalues. (Input) PLACE can be one of the following enumerated integers as defined in ARPACK_INT:
Value
ARPACK_Largest_Magnitude
ARPACK_Smallest_Magnitude
ARPACK_Largest_Real_Parts
ARPACK_Smallest_Real_Parts
ARPACK_Largest_Imag_Parts
ARPACK_Smallest_Imag_Parts
Default: ARPACK_Largest_Magnitude.
TYPE — Defines the eigenvalue problem as either a standard or generalized eigenvalue problem. (Input) TYPE can be one of the following enumerated integers as defined in ARPACK_INT:
Value
Description
ARPACK_Standard
Ax = λx
ARPACK_Generalized
Ax = λBx
Default: TYPE = ARPACK_Standard.
CATEGORY —CATEGORY and TYPE define the operation sequence provided in the user-written function. (Input)
CATEGORY can be one of the following enumerated integers as defined in ARPACK_INT:
Value
Description
ARPACK_Regular
Ax = λx
ARPACK_Regular_Inverse
y=Ax, y=Bx, y=B-1x
ARPACK_Shift_Invert
y=Ax, y= (A-σ I)-1x, y = (A-σ B)-1x
ARPACK_Complex_Part_Shift_Invert
y=Ax, y= (A-σ I)-1x, y = (A-σ B)-1x y=Ax,y = Bx, y = Im{(A-σ B)-1x} y=Ax,y = Bx, y = Re{(A-σ B)-1x}
Default: CATEGORY = ARPACK_Regular.
EXTYPE — A derived type of the extensible class ARPACKBASE, which may be used to pass additional information to/from the user-supplied function. (Input/Output) The user must include a USEARPACK_INT statement in the calling program to define this derived type. If EXTYPE is not included as an argument to ARPACK_NONSYMMETRIC it must still be supplied as an argument to user-function, F, but is not used.
VECTORS— An allocatable array of approximate eigenvectors. (Output) It is not necessary to allocate VECTORS(:,:). If this argument is used the allocation occurs within the routine ARPACK_NONSYMMETRIC. The output sizes are VECTORS(1:N,1:NCV). The second dimension value is NCV=min(N, max(FACTOR_MAXNCV*NEV,NEV+1)), where the value FACTOR_MAXNCV is a component of the base class, ARPACKBASE. The first NEV columns of VECTORS(:,:) represent the eigenvectors (see Comments).
FORTRAN 2003 Interface
Generic: ARPACK_NONSYMMETRIC (N, F,ZVALUES [,…])
Specific: The specific interface name is D_ARPACK_NONSYMMETRIC.
Fortran 90 Interface
A Fortran 90 compiler version is not available for this routine.
Description
Routine ARPACK_NONSYMMETRIC calls ARPACK subroutines to compute partial eigenvalue-eigenvector decompositions for real matrices. The ARPACK routines are dnaupd and dneupd (see ARPACK Users’ Guide, SIAM Publications, (1998)), which use “reverse communication” to obtain the required matrix-vector operations for this approximation. Responses to these requests are made by calling the user-written function F. By including the class object EXTYPE as an argument to this function, user data and procedure pointers are available for the evaluations. A user code must extend the base class EXTYPE to include the extra data and procedure pointers.
Comments
The non-symmetric problem may have complex eigenvalues that occur in conjugate pairs, and the eigenvectors are returned in the REAL(DKIND) array VECTORS(:,:) but with a compact representation: If the eigenvalue λj has an imaginary part with a negative value, construct the complex eigenvector from the relation wj = vj+ivj+1. The real vectors vj, vj+1 are consecutive columns of the array VECTORS (:,:). The eigenvalue-eigenvector relationship is Awj = λjwj. Since A is real, λ is also an eigenvalue; thus the conjugate relationshipwill hold. For purposes of checking results the complex residual rj = Awj -λjwj should be small in norm relative to the norm of A. If that is true, there is no need to check the alternate relationship. This compact representation of the eigenvectors can be expanded to require twice the storage requirements, but that is not done here in the interest of saving large blocks of storage.
For the generalized eigenvalue problem Ax = λBx the eigenvalues are optionally computed based on the Raleigh Quotient. Because of the shifts used, only the eigenvectors may be computed. The eigenvalues are returned by solving Awj = λBwj for λ: . λj is valid if the denominator is non-zero. If λj has a non-zero imaginary part, then the complex conjugate is also an eigenvalue. The Raleigh Quotient for eigenvalues of generalized problems is used when vectors are requested and the user has requested it be used with the base class component EXTYPE%RALEIGH_QUOTIENT == .TRUE. This is the component’s default value.
Example
We solve the generalized eigenvalue problem Ax = λBx using the shift-invert category. The matrix A is tri-diagonal with the values 2 on the diagonal, -2 on the sub-diagonal, and 3 on the super-diagonal. The matrix B is tri-diagonal with the values 4 on the diagonal and 1 on the off-diagonals. We use the complex shift σ = 0.4 + 0.6i and increase the factor for the number of Ritz vectors from 2.5 to 5. Two strategies of shift-invert are illustrated, y = Re(A-σB)-1Bx and y = Im(A-σB)-1Bx. In each case NEV=6 eigenvalues are obtained, each with 3 pairs of complex conjugate values.
(Example arpack_nonsymmetric_ex1.f90)
MODULE ARPACK_NONSYMMETRIC_EX1_INT
USE ARPACK_INT
USE LSLCQ_INT
USE N1RTY_INT
IMPLICIT NONE
TYPE, EXTENDS(ARPACKBASE) :: ARPACKBASE_EXT
INTEGER :: NX=0
INTEGER :: NV=0
! This example extends the base type to
! information for solving complex tridiagonal systems.
COMPLEX(DKIND), ALLOCATABLE :: A(:), B(:), C(:)
INTEGER, ALLOCATABLE :: IR(:), IS(:)
! This controls the type of shifting. When
! the value is 1, use real part of inv(A-*M)*x.
! If value is 2, use imaginary part of same.
INTEGER :: SHIFT_STRATEGY=1
END TYPE ARPACKBASE_EXT
CONTAINS
FUNCTION FCN(X, TASK, EX) RESULT(Y)
USE UMACH_INT
CLASS (ARPACKBASE), INTENT(INOUT) :: EX
REAL (DKIND), INTENT(INOUT) :: X(:)
INTEGER, INTENT(IN) :: TASK
INTEGER, PARAMETER :: NSIZE=12
REAL(DKIND) Y(size(X))
REAL(DKIND) DL, DD, DU
COMPLEX(DKIND) Z(2*size(X))
REAL(DKIND) U(2*size(X))
INTEGER J, IERR, NOUT, IJOB
CALL UMACH(2, NOUT)
SELECT TYPE(EX)
TYPE IS(ARPACKBASE_EXT)
ASSOCIATE(N => EX % NCOLS,&
NV => EX % NV, &
SHIFT => EX % SHIFT)
SELECT CASE(TASK)
CASE(ARPACK_A_x)
DL = -2._DKIND
DD = 2._DKIND
DU = 3._DKIND
Y(1) = DD*X(1) + DU*X(2)
DO J = 2,N-1
Y(J) = DL*X(J-1) + DD*X(J) + DU*X(J+1)
END DO
Y(N) = DL*X(N-1) + DD*X(N)
NV=NV+1
CASE(ARPACK_B_x)
Y(1) = 4._DKIND*X(1) + X(2)
DO J = 2,N-1
Y(J) = X(J-1) + 4._DKIND*X(J) + X(J+1)
END DO
Y(N) = X(N-1) + 4._DKIND*X(N)
NV=NV+1
CASE(ARPACK_inv_A_minus_Shift_B_x)
! Compute Y=REAL/AIMAG(inv(A-*B)*x). This is done with a solve
! step, using the LU factorization. Note that the data
! for the factorization is stored in the user's extended
! data type.
Z=CMPLX(X,0._DKIND,DKIND)
IJOB = 2
CALL LSLCQ (EX%C, EX%A, EX%B, Z, U, &
EX%IR, EX%IS, N=N, IJOB=IJOB)
IERR= N1RTY(1)
IF (IERR==4 .OR. IERR==5) &
STOP 'IMSl_FATAL_ERROR_SOLVING'
IF(EX % SHIFT_STRATEGY == 1) THEN
Y(1:N)=REAL( Z(1:N),DKIND)
ELSE IF (EX % SHIFT_STRATEGY == 2)THEN
Y(1:N)=AIMAG(Z(1:N))
END IF
! Total number of solve steps.
NV=NV+1
CASE(ARPACK_Prepare)
! Set up storage areas for factored tridiagonal matrix.
IF (ALLOCATED(EX%A)) DEALLOCATE(EX%A)
IF (ALLOCATED(EX%B)) DEALLOCATE(EX%B)
IF (ALLOCATED(EX%C)) DEALLOCATE(EX%C)
IF (ALLOCATED(EX%IR)) DEALLOCATE(EX%IR)
IF (ALLOCATED(EX%IS)) DEALLOCATE(EX%IS)
ALLOCATE(EX%A(2*N), EX%B(2*N), EX%C(2*N), &
EX%IR(NSIZE), EX%IS(NSIZE), STAT=IERR)
IF (IERR /= 0) STOP 'IMSL_ERROR_ALLOCATING_WORKSPACE'
