Compute some eigenvalues and eigenvectors of the
generalized eigenvalue problem. This can be used for
the case
. The values
for
are real, but eigenvalues
may be complex and occur in 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 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 |
|
ARPACK_B_x |
|
ARPACK_inv_A_minus_Shift_x |
|
ARPACK_inv_B_x |
|
ARPACK_inv_A_minus_Shift_B_x |
|
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_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.
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 |
|
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_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).
Generic: ARPACK_NONSYMMETRIC (N, F,ZVALUES [,…])
Specific: The specific interface name is D_ARPACK_NONSYMMETRIC.
A Fortran 90 compiler version is not available for this routine.
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.
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 eigenvaluehas an imaginary part with a negative value, construct the
complex eigenvector from the relation
. The real vectors
are consecutive columns of the array VECTORS
(:,:). The eigenvalue-eigenvector relationship is
. Since
is real,
is also an eigenvalue; thus the conjugate relationship
will hold. For purposes of checking results the complex
residual
should be small in norm
relative to the norm of
. 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 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
for
:
.
is valid if the
denominator is non-zero. If
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.
We solve the generalized eigenvalue problem using the shift-invert category. The matrix
is tri-diagonal with the values 2 on the diagonal, -2 on
the sub-diagonal, and 3 on the super-diagonal. The matrix
is tri-diagonal with the values 4 on the diagonal and 1 on
the off-diagonals. We use the complex shift
and increase the factor for the number of Ritz vectors from 2.5
to 5. Two strategies of shift-invert are illustrated,
and
. In each case NEV=6
eigenvalues are obtained, each with 3 pairs of complex conjugate values.
Link to example source (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'
! Define matrix, A-SHIFT*B.
EX % B(1:N) = -2._DKIND-SHIFT
EX % A(1:N) = 2._DKIND-4._DKIND*SHIFT
EX % C(1:N) = 3._DKIND-SHIFT
! Factor the matrix with LU and partial pivoting.
IJOB = 3
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'
! Give output some ZVALUES to satisfy compiler.
Y=0._DKIND
NV=0
CASE DEFAULT
WRITE(NOUT,*) TASK, ': INVALID OPERATION REQUESTED'
STOP 'IMSL_ERROR_WRONG_OPERATION'
END SELECT
END ASSOCIATE
END SELECT
END FUNCTION
END MODULE
! Suppose we want to solve A*x = lambda*B*x in -invert mode
! The matrix A is the tridiagonal matrix with 2 on the diagonal,
! -2 on the subdiagonal and 3 on the superdiagonal. The matrix
! is the tridiagonal matrix with 4 on the diagonal and 1 on the
! off-diagonals.
! The sigma is a complex number (sigmar, sigmai).
! OP = Real_Part{invA-(SIGMAR,SIGMAI)*B*B.
USE ARPACK_NONSYMMETRIC_EX1_INT
USE UMACH_INT
USE WRCRN_INT
INTEGER, PARAMETER :: NEV=6, N=100
COMPLEX(DKIND) :: ZVALUES(NEV), RES(N),U(N),V(N),W(N)
REAL(DKIND), ALLOCATABLE :: VECTORS(:,:)
REAL(DKIND) NORM
LOGICAL SKIP, SMALL, SOLVED
INTEGER J, STRATEGY, NOUT
CHARACTER(LEN=12) TAG
CHARACTER(LEN=60) TITLE
TYPE(ARPACKBASE_EXT) EX
ASSOCIATE(NX => EX % NX, &
NV => EX % NV, &
SHIFT => EX % SHIFT,&
FACTOR => EX % FACTOR_MAXNCV,&
NACC => EX % NACC)
! Note that VECTORS(:,:) does not need to be allocated
! in the calling program. That happens within the
! routine ARPACK_NONSYMMETRIC(). It is OK to do this but
! the sizes (N,NCV) are determined in ARPACK_NONSYMMETRIC.
CALL UMACH(2, NOUT)
SOLVED=.TRUE.
DO STRATEGY=1,2
SHIFT=CMPLX(0.4_DKIND,0.6_DKIND,DKIND)
FACTOR=5._DKIND
EX % SHIFT_STRATEGY=STRATEGY
CALL ARPACK_NONSYMMETRIC(N, FCN, ZVALUES, &
TYPE=ARPACK_Generalized, &
CATEGORY=ARPACK_Complex_Part_Shift_Invert, &
EXTYPE=EX, VECTORS=VECTORS)
WRITE(NOUT, *) &
'Number of Matrix-Vector Products Required, NS Ex-1'
WRITE(NOUT, *) &
'--------------------------------------------------'
WRITE(NOUT, '(5X,I4)') NV
WRITE(NOUT, *) 'Number of accurate values determined'
WRITE(NOUT, *) '------------------------------------'
WRITE(NOUT, '(5X, I4)') NACC
! Check residuals, A*vectors = ZVALUES*M*vectors:
SKIP=.FALSE.
DO J=1,NACC
IF(SKIP) THEN
SKIP=.FALSE.
CYCLE
END IF
! The eigenvalue is complex and the pair of vectors
! for the complex eigenvector is returned.
IF(AIMAG(ZVALUES(J)) /= 0._DKIND)THEN
! Make calls for real and imaginary parts of eigenvectors
! applied to the operators A, M
U=CMPLX(FCN(VECTORS(:,J),ARPACK_A_x,EX),&
FCN(VECTORS(:,J+1),ARPACK_A_x,EX),DKIND)
V=CMPLX(FCN(VECTORS(:,J),ARPACK_B_x,EX),&
FCN(VECTORS(:,J+1),ARPACK_B_x,EX),DKIND)
! Since the matrix is real, there is an additional conjugate:
RES=U-ZVALUES(J)*V
SKIP=.TRUE.
ELSE
! The eigenvalue is real and the real eigenvector is returned.
RES=FCN(VECTORS(:,J),ARPACK_A_x,EX)-ZVALUES(J)*&
FCN(VECTORS(:,J),ARPACK_B_x,EX)
END IF
NORM=maxval(abs(RES))
SMALL=(NORM <= ABS(ZVALUES(J))*SQRT(EPSILON(NORM)))
SOLVED=SOLVED .and. SMALL
END DO
IF(STRATEGY==1) TAG='REAL SHIFT'
IF(STRATEGY==2) TAG='IMAG SHIFT'
TITLE = 'Largest Raleigh Quotient Eigenvalues,'//TAG
CALL WRCRN(TITLE, ZVALUES)
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 DO ! Shift strategy
END ASSOCIATE
END
Number of Matrix-Vector Products Required, NS Ex-1
--------------------------------------------------
280
Number of accurate values determined
------------------------------------
6
Largest Raleigh Quotient Eigenvalues,REAL SHIFT
1 ( 0.5000,-0.5958)
2 ( 0.5000, 0.5958)
3 ( 0.5000,-0.6331)
4 ( 0.5000, 0.6331)
5 ( 0.5000, 0.5583)
6 ( 0.5000,-0.5583)
All Ritz Values and Vectors have small residuals.
Number of Matrix-Vector Products Required, NS Ex-1
--------------------------------------------------
248
Number of accurate values determined
------------------------------------
6
Largest Raleigh Quotient Eigenvalues,IMAG SHIFT
1 ( 0.5000,-0.5958)
2 ( 0.5000, 0.5958)
3 ( 0.5000,-0.5583)
4 ( 0.5000, 0.5583)
5 ( 0.5000,-0.6331)
6 ( 0.5000, 0.6331)
All Ritz Values and Vectors have small residuals.
PHONE: 713.784.3131 FAX:713.781.9260 |