Computes some eigenvalues and eigenvectors of the generalized real symmetric eigenvalue problem Ax = λBx. This can be used for the case B = I.
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_SYMMETRIC 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.
VALUES — An array of eigenvalues. (Output) The value NEV=size(VALUES) defines the number of eigenvalues to be computed. The calling program declares or allocates the array VALUES(1:NEV). The number of eigenvalues computed accurately is optionally available as the component EXTYPE%NACC of the base class EXTYPE.
Optional Arguments
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_Algebraic
ARPACK_Smallest_Algebraic
ARPACK_Largest_Magnitude
ARPACK_inv_A_minus_Shift_x
ARPACK_Smallest_Magnitude
ARPACK_Both_Ends
Default: PLACE = ARPACK_Largest_Algebraic.
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
y=Ax
ARPACK_Regular_Inverse
y=Ax, y=Bx, y=B-1x
ARPACK_Shift_Invert
y=Ax, y= (A -σ I)-1x
ARPACK_Buckling
y=Ax, y=Bx, y= (A-σ B)-1x
ARPACK_Cayley
y=Ax, y=Bx, y= (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_SYMMETRIC 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_SYMMETRIC. 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(:,:) are the eigenvectors.
FORTRAN 2003 Interface
Generic: ARPACK_SYMMETRIC (N, F, VALUES [,…])
Specific: The specific interface name is D_ARPACK_SYMMETRIC.
Fortran 90 Interface
A Fortran 90 compiler version is not available for this routine.
Description
Routine ARPACK_SYMMETRIC calls ARPACK subroutines to compute partial eigenvalue-eigenvector decompositions for symmetric real matrices. The ARPACK routines are dsaupd and dseupd (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 user function F is written to supply requests for the matrix operations. The following psuedo-code outlines the required responses of F depending on the circumstances. Only those cases that follow from the settings of PLACE, TYPE and CATEGORY need to be provided in the user code. The enumerated constants, ARPACK_A_x, etc., are available by use-association from the module ARPACK_INT.
FUNCTION F (X, TASK, EXTYPE) RESULT(Y)
USE ARPACK_INT
IMPLICIT NONE
CLASS(ARPACKBASE), INTENT(INOUT) :: EXTYPE
REAL(DKIND), INTENT(INOUT) :: X(:)
INTEGER, INTENT(IN) :: TASK
REAL(DKIND) Y(SIZE(X))
SELECT CASE (TASK)
CASE (ARPACK_Prepare)
…{Take initial steps to prepare for the operations that follow.}
CASE (ARPACK_A_x)
CASE (ARPACK_B_x)
CASE (ARPACK_inv_A_minus_Shift_x)
CASE (ARPACK_inv_B_x)
CASE (ARPACK_inv_A_minus_Shift_B_x)
CASE DEFAULT
…{This is an error condition. }
END SELECT
END FUNCTION
Examples
Example 1
We approximate eigenvalues and eigenfunctions of the Laplacian operator
on the unit square, [0,1] × [0,1], with zero Dirichlet boundary values. The full set of eigenvalues and their eigenfunctions are given by the sequence λm,n= (m2+n2)π2, um,n(x,y) = sin(mπ) sin(nπ), where m,n are positive integers.
This provides a check on the accuracy of the numerical results. Equally spaced divided differences on the unit square are used to approximate -Δu. A few eigenvalues of smallest magnitude, and their eigenvectors, are requested. This application requires the optional argument PLACE=ARPACK_Smallest_Magnitude. The user function code provides the second order divided difference operator applied to an input vector. The problem is a symmetric matrix eigenvalue computation. It involves only the single TASK, ARPACK_A_x, in the user functions.
The function FCN defines a grid of values and provides the operation that derives from this eigenvalue problem. The class argument EXTYPE must be declared but need not be used. Within the main program, the function interface for the external function FCN is specified with the declaration PROCEDURE (DMV) FCN.
(Example arpack_symmetric_ex1.f90)
PROGRAM ARPACK_SYMMETRIC_EX1
USE ARPACK_INT
USE UMACH_INT
USE WRRRN_INT
IMPLICIT NONE
! Compute the smallest eigenvalues of a discrete Laplacian,
! based on second order divided differences.
! The matrix used is the 2 dimensional discrete Laplacian on
! the unit square with zero Dirichlet boundary condition.
INTEGER :: J, NOUT
INTEGER, parameter :: NEV=5 !number of Eigenvalues required
INTEGER, parameter :: NV=0, NX=10
INTEGER, parameter :: N=NX**2 !size of matrix problem
REAL(DKIND) :: VALUES(NEV), RES(N), EF(NX, NX)
REAL(DKIND), ALLOCATABLE :: VECTORS(:,:)
REAL(DKIND) :: NORM
LOGICAL :: SMALL, SOLVED
TYPE(ARPACKBASE) :: Q
PROCEDURE(DMV) :: FCN
CALL UMACH(2, NOUT)
! Note that VECTORS(:,:) does not need to be allocated
! in the calling program. That happens within the
! routine ARPACK_SYMMETRIC(). It is OK to do this but
! the sizes (N,NCV) are determined in ARPACK_SYMMETRIC.
CALL ARPACK_SYMMETRIC(N, FCN, VALUES, &
PLACE=ARPACK_Smallest_Magnitude, VECTORS=VECTORS)
WRITE(NOUT, *) 'Number of eigenvalues requested, and declared accurate'
We approximate eigenvalues and eigenfunctions of the 1D Laplacian operator on the unit interval, [0,1]. Equally spaced divided differences are used for the operator, which yields a tri-diagonal matrix. The eigenvalues and eigenfunctions are λn=n2π2, un(x) = sin(nπ), n = 1,2,.... This example shows that using inverse iteration for approximating the largest reciprocals of eigenvalues is more efficient than directly computing the smallest magnitude eigenvalues by products of the operator. This requires the optional argument CATEGORY=ARPACK_Shift_Invert. The user function, FCN, requires the solution of a tri-diagonal system of linear equations applied to an input vector. The base class ARPACKBASE is extended to the user’s type, TYPE(ARPACKBASE_EXT), defined in the user module ARPACK_SYMMETRIC_EX2_INT. This extension includes the number of intervals, a total kept in FCN for noting the number of operations, and allocatable arrays used to store the LU factorization of the tri-diagonal matrix. When FCN is entered with TASK=ARPACK_Prepare, these arrays are allocated, defined, and the LU factorization of the shifted matrix (A-σ I) is computed, here with σ= 0. When FCN is later entered with TASK=ARPACK_inv_A_minus_Shift_x, the LU factorization is available to efficiently compute y= (A‑σ I)-1x=A-1x. The function FCN is also entered with TASK=ARPACK_A_x, to compute Ax.
(Example arpack_symmetric_ex2.f90)
MODULE ARPACK_SYMMETRIC_EX2_INT
USE ARPACK_INT
USE LSLCR_INT
USE N1RTY_INT
IMPLICIT NONE
TYPE, EXTENDS(ARPACKBASE) :: ARPACKBASE_EXT
REAL(DKIND) :: HSQ=0._DKIND
INTEGER :: NX=0, NV=0
! This example extends the base type to
! information for solving tridiagonal systems.
REAL(DKIND), ALLOCATABLE :: A(:), B(:), C(:)
REAL(DKIND), ALLOCATABLE :: Y1(:), U(:)
INTEGER, ALLOCATABLE :: IR(:), IS(:)
END TYPE ARPACKBASE_EXT
CONTAINS
FUNCTION FCN(X, TASK, EX) RESULT(Y)
CLASS (ARPACKBASE), INTENT(INOUT) :: EX
REAL (DKIND), INTENT(INOUT) :: X(:)
INTEGER, INTENT(IN) :: TASK
REAL (DKIND) Y(size(X))
INTEGER J, IERR, IJOB, NSIZE
SELECT TYPE(EX)
TYPE IS(ARPACKBASE_EXT)
ASSOCIATE(N => EX%NCOLS, &
NV => EX%NV, &
HSQ => EX%HSQ, &
SHIFT => EX%SHIFT)
SELECT CASE(TASK)
CASE(ARPACK_A_x)
Y(1) = 2._DKIND*X(1) - X(2)
DO J = 2,N-1
Y(J) = - X(J-1) + 2._DKIND*X(J) - X(J+1)
END DO
Y(N) = - X(N-1) + 2._DKIND*X(N)
Y=Y/HSQ
CASE(ARPACK_inv_A_minus_Shift_x)
! Compute Y=inv(A-*I)*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.
EX%Y1(1:N) = X
IJOB = 2
CALL LSLCR (EX%C, EX%A, EX%B, EX%Y1, EX%U, &
EX%IR, EX%IS, N=N, IJOB=IJOB)
Y(1:N) = EX%Y1(1:N)
IERR= N1RTY(1)
IF (IERR==4 .OR. IERR==5) STOP 'IMSl_FATAL_ERROR_SOLVING'
! Total number of solve steps.
NV=NV+1
CASE(ARPACK_Prepare)
! Set up storage areas for factored tridiagonal matrix.
'All Ritz Values and Vectors have small residuals.'
ELSE
WRITE(nout,'(A///)') &
'Some Ritz Values and Vectors have large residuals.'
ENDIF
END ASSOCIATE
END
Output
Number of Matrix-Vector Products Required, EX-2
-----------------------------------------------
24
Largest Laplacian Eigenvalues Near Zero Shift
1 9.9
2 39.5
3 88.8
4 157.7
All Ritz Values and Vectors have small residuals.
Example 3
We compute the solution of a generalized problem. This comes from using equally spaced linear finite element test functions to solve eigenvalues and eigenfunctions of the 1D Laplacian operator on the unit interval, [0,1]. This is Example 2 but solved using finite elements. With matrix notation, we have the matrix problem Ax = λBx. Both A and B are tri-diagonal and symmetric. The matrix B is non-singular. We compute the smallest magnitude eigenvalues. This requires the optional arguments TYPE = ARPACK_Generalized, CATEGORY = ARPACK_Regular_Inverse, and PLACE = ARPACK_Smallest_Magnitude. The user function, FCN, requires the solution of a tri-diagonal system of linear equations applied to an input vector, y = B-1x. The base class ARPACKBASE is extended to the user’s type, TYPE(ARPACKBASE_EXT), defined in the user module ARPACK_SYMMETRIC_EX3_INT. This extension includes the number of intervals, a total kept in FCN for noting the number of operations, and allocatable arrays used to store the LU factorization of B. When FCN is entered with TASK=ARPACK_Prepare, these arrays are allocated, defined, and the LU factorization of the matrix B is computed. The function FCN is entered with the three values TASK=ARPACK_A_x, for y = Ax; TASK=ARPACK_B_x, for y = Bx; and TASK=ARPACK_inv_B_x, for y = B-1x.
Within the main program, the function interface for the external function FCN is specified with the declaration PROCEDURE (DMV) FCN.
(Example arpack_symmetric_ex3.f90)
MODULE ARPACK_SYMMETRIC_EX3_INT
USE ARPACK_INT
USE LSLCR_INT
USE N1RTY_INT
IMPLICIT NONE
TYPE, EXTENDS(ARPACKBASE) :: ARPACKBASE_EXT
REAL(DKIND) :: H=0._DKIND
INTEGER :: NX=0, NV=0
! This example extends the base type to
! information for solving tridiagonal systems.
REAL(DKIND), ALLOCATABLE :: A(:), B(:), C(:)
REAL(DKIND), ALLOCATABLE :: Y1(:), U(:)
INTEGER, ALLOCATABLE :: IR(:), IS(:)
END TYPE ARPACKBASE_EXT
END MODULE
PROGRAM ARPACK_SYMMETRIC_EX3
USE ARPACK_SYMMETRIC_EX3_INT
USE UMACH_INT
USE WRRRN_INT
IMPLICIT NONE
! We want to solve A*x = lambda*M*x in inverse mode,
! where A and M are obtained by the finite element method
! of the 1-dimensional discrete Laplacian
! d^2u / dx^2
! on the interval 0,1, with zero Dirichlet boundary conditions,
! using piecewise linear elements.
INTEGER,PARAMETER :: NEV=4, N=100
REAL(DKIND) :: VALUES(NEV), RES(N)
REAL(DKIND), ALLOCATABLE :: VECTORS(:,:)
REAL(DKIND) NORM
LOGICAL :: PRINTRESULTS = .FALSE.
LOGICAL SMALL, SOLVED
INTEGER J, NOUT
PROCEDURE(DMV) FCN
TYPE(ARPACKBASE_EXT) EX
ASSOCIATE(NX => EX%NX, &
NV => EX%NV)
EX%FACTOR_MAXNCV=5._DKIND
CALL UMACH(2, NOUT)
! Note that VECTORS(:,:) does not need to be allocated
! in the calling program. That happens within the
! routine ARPACK_SYMMETRIC(). It is OK to do this but
! the sizes (N,NCV) are determined in ARPACK_SYMMETRIC.
