Compute some eigenvalues and eigenvectors of the
generalized eigenvalue problem. This can be used for
the case
. The values for
are real or complex. When the values are complex the
eigenvalues may be complex and are not expected to occur in complex conjugate
pairs.
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 ZMV(…). 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 |
|
ARPACK_xt_A |
|
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 USE ARPACK_INT statement
in the calling program to define this derived type. If EXTYPE is not included
as an argument to ARPACK_COMPLEX it
should be ignored in the user-function, F.
The function F must be written according to the abstract interface for ZMV. If F is not use-associated nor contained in the calling program, declare it with PROCEDURE(ZMV) 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
number of eigenvalues computed accurately is optionally available as the
component EXTYPE%NACC of the
base class EXTYPE.
Place — Defines the output
content of VALUES.
(Input)
PLACE
specifies the placement within the spectrum for the required eigenvalues.
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: PLACE = 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 |
|
ARPACK_Generalized |
|
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 |
|
ARPACK_Regular_Inverse |
|
ARPACK_Shift_Invert |
|
ARPACK_Complex_Part_Shift_Invert |
|
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 USE ARPACK_INT statement
in the calling program to define this derived type. If EXTYPE is not included
as an argument to ARPACK_COMPLEX 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).
Generic: ARPACK_COMPLEX (N, F,ZVALUES [,…])
Specific: The specific interface name is Z_ARPACK_COMPLEX.
A Fortran 90 compiler version is not available for this routine.
Routine ARPACK_COMPLEX calls ARPACK subroutines to compute partial eigenvalue-eigenvector decompositions for complex matrices. The ARPACK routines are dzaupd and dzeupd (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.
For purposes of checking results the complex
residualshould be small in norm
relative to the norm of
. For the
generalized eigenvalue problem
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
based on solving
for
,
where
.
The eigenvalue is finite and valid if the denominator is non-zero.
The Raleigh Quotient for eigenvalues of generalized problems is used only 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.
This example is a variation of the first example for ARPACK_SYMMETRIC. We approximate eigenvalues and eigenfunctions of the Laplacian operator
,
on the unit square, . But now the parameter
is complex. Thus the eigenvalues and eigenfunctions are
complex.
Link to example source (arpack_complex_ex1.f90)
MODULE ARPACK_COMPLEX_EX1_INT
USE ARPACK_INT
IMPLICIT NONE
TYPE, EXTENDS(ARPACKBASE), PUBLIC :: ARPACKBASE_EXT
REAL(DKIND) :: H =0._DKIND
REAL(DKIND) :: HSQ=0._DKIND
COMPLEX(DKIND) :: RHO=(0._DKIND,0._DKIND)
COMPLEX(DKIND) :: DL=(0._DKIND,0._DKIND)
COMPLEX(DKIND) :: DD=(0._DKIND,0._DKIND)
COMPLEX(DKIND) :: DU=(0._DKIND,0._DKIND)
INTEGER :: NX=0
INTEGER :: NV=0
END TYPE ARPACKBASE_EXT
CONTAINS
FUNCTION FZ1(X, TASK, EXTYPE) RESULT(Y)
USE UMACH_INT
CLASS (ARPACKBASE), INTENT(INOUT) :: EXTYPE
COMPLEX (DKIND), INTENT(INOUT) :: X(:)
INTEGER, INTENT(IN) :: TASK
COMPLEX (DKIND) Y(size(X))
COMPLEX (DKIND) DT(3)
REAL(DKIND) :: ONE=1._DKIND
INTEGER J, NOUT
CALL UMACH(2, NOUT)
SELECT TYPE(EXTYPE)
TYPE IS(ARPACKBASE_EXT)
ASSOCIATE(NX => EXTYPE % NX,&
H => EXTYPE % H,&
HSQ => EXTYPE % HSQ,&
RHO => EXTYPE % RHO,&
DL => EXTYPE % DL,&
DD => EXTYPE % DD,&
DU => EXTYPE % DU,&
NV => EXTYPE % NV)
SELECT CASE(TASK)
CASE(ARPACK_A_x)
! Computes y <-- A*x, where A is the N**2 by N**2 block
! tridiagonal matrix deriving from (Laplacian u) + rho*(du/dx).
DT=(/DL,DD,DU/)
Y(1:NX)=T(NX,X(1:NX),DT) - X(NX+1:2*NX)/HSQ
DO J=NX+1,NX**2-NX,NX
Y(J:J+NX-1)=T(NX,X(J:J+NX-1),DT) &
- (X(J-NX:J-1)+ X(J+NX:J+2*NX-1))/HSQ
END DO
Y((NX-1)*NX+1:NX**2)= - X((NX-1)*NX-NX+1:(NX-1)*NX) &
/ HSQ + T(NX,X((NX-1)*NX+1:NX**2),DT)
! Total the number of matrix-vector products.
NV=NV+1
CASE(ARPACK_Prepare)
! Define 1/H**2, etc. so they are available in the evaluator.
H = ONE/REAL(NX+1,DKIND)
HSQ = H**2
DD = (4.0D+0, 0.0D+0) / HSQ
DL = -ONE/HSQ - (5.0D-1, 0.0D+0) *RHO/H
DU = -ONE/HSQ + (5.0D-1, 0.0D+0) *RHO/H
NV = 0
CASE DEFAULT
WRITE(nout,*) TASK, ': INVALID TASK REQUESTED'
STOP 'IMSL_ERROR_WRONG_OPERATION'
END SELECT
END ASSOCIATE
END SELECT
END FUNCTION
FUNCTION T(NX, X, DT)RESULT(V)
INTEGER, INTENT(IN) :: NX
COMPLEX(DKIND), INTENT(IN) :: X(:), DT(3)
COMPLEX(DKIND) :: V(NX)
INTEGER J
ASSOCIATE(DL => DT(1),&
DD => DT(2),&
DU => DT(3))
V(1) = DD*X(1) + DU*X(2)
DO J = 2,NX-1
V(J) = DL*X(J-1) + DD*X(J) + DU*X(J+1)
END DO
V(NX) = DL*X(NX-1) + DD*X(NX)
END ASSOCIATE
END FUNCTION
END MODULE
! Compute the largest magnitude eigenvalues of a discrete Laplacian,
! based on second order divided differences.
! The matrix used is obtained from the standard central difference
! discretization of the convection-diffusion operator
! (Laplacian u) + rho*(du / dx)
! on the unit squre 0,1x0,1 with zero Dirichlet boundary
! conditions.
USE ARPACK_COMPLEX_EX1_INT
USE UMACH_INT
USE WRCRN_INT
INTEGER, PARAMETER :: NEV=6
INTEGER :: J, N, NOUT
COMPLEX(DKIND) :: VALUES(NEV)
COMPLEX(DKIND), ALLOCATABLE :: RES(:), EF(:,:)
COMPLEX(DKIND), ALLOCATABLE :: VECTORS(:,:)
REAL(DKIND) NORM
LOGICAL SMALL, SOLVED
TYPE(ARPACKBASE_EXT) EX
ASSOCIATE(NX => EX % NX, &
NV => EX % NV, &
RHO => EX % RHO,&
NACC => EX % NACC)
CALL UMACH(2, NOUT)
NX=10
RHO=(100._DKIND,1._DKIND)
! Define size of matrix problem.
N=NX**2
! Note that VECTORS(:,:) does not need to be allocated
! in the calling program. That happens within the
! routine ARPACK_COMPLEX(). It is OK to do this but
! the sizes (N,NCV) are determined in ARPACK_COMPLEX.
CALL ARPACK_COMPLEX(N, FZ1, VALUES, EXTYPE=EX, VECTORS=VECTORS)
WRITE(NOUT, *) 'Number of eigenvalues requested, and accurate'
WRITE(NOUT, *) '---------------------------------------------'
WRITE(NOUT, '(5X, I4, 5X, I4)') NEV, NACC
WRITE(NOUT, *) 'Number of Matrix-Vector Products Required, ZEX-1'
WRITE(NOUT, *) '------------------------------------------------'
WRITE(NOUT, '(5X, I4)') NV
CALL WRCRN ('Largest Magnitude Operator Eigenvalues', VALUES)
! Check residuals, A*vectors = values*vectors:
ALLOCATE(RES(N))
DO J=1,NACC
RES=FZ1(VECTORS(:,J),ARPACK_A_x,EX)-VALUES(J)*VECTORS(:,J)
NORM=maxval(abs(RES))
SMALL=(NORM <= ABS(VALUES(J))*SQRT(EPSILON(NORM)))
IF(J==1) SOLVED=SMALL
SOLVED=SOLVED .and. SMALL
END DO
IF(SOLVED) THEN
WRITE(NOUT,'(A///)') &
'All Ritz Values and Vectors have small residuals.'
ELSE
WRITE(NOUT,'(A///)') &
'Some Ritz Values and Vectors have large residuals.'
END IF
END ASSOCIATE
END
Number of eigenvalues requested, and accurate
---------------------------------------------
6 6
Number of Matrix-Vector Products Required, ZEX-1
------------------------------------------------
475
Largest Magnitude Operator Eigenvalues
1 ( 727.0,-1029.6)
2 ( 705.4, 1029.6)
3 ( 698.4,-1029.6)
4 ( 676.8, 1029.6)
5 ( 653.3,-1029.6)
6 ( 631.7, 1029.6)
All Ritz Values and Vectors have small residuals.
PHONE: 713.784.3131 FAX:713.781.9260 |