ARPACK_SYMMETRIC

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 USE ARPACK_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.

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'

WRITE(NOUT, *) '------------------------------------------------------'

WRITE(NOUT, '(5X, I4, 5X, I4)') NEV, Q%NACC

WRITE(NOUT, *) 'Number of Matrix-Vector Products Recorded, EX-11'

WRITE(NOUT, *) '------------------------------------------------'

WRITE(NOUT, '(5X, I4)') NV

CALL WRRRN('Smallest Laplacian Eigenvalues', VALUES)

! Check residuals, A*vectors = values*vectors:

DO J=1,NEV

! Necessary to have an unused TYPE(ARPACKBASE) :: Q as an argument:

RES=FCN(VECTORS(:,J), ARPACK_A_x, Q)-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.'

ENDIF

! The first eigenvector is scaled to be positive.

! It defines the eigenfunction for the smallest

! eigenvalue at the grid defined by the differencing.

EF=reshape(VECTORS(:,1),(/NX,NX/))

CALL WRRRN('First 2D Laplacian Eigenfunction at Grid Points', EF)

END

FUNCTION FCN(X, TASK, EX)RESULT(Y)

USE ARPACK_INT

CLASS(ARPACKBASE),INTENT(INOUT) :: EX

REAL(DKIND), INTENT(INOUT) :: X(:)

INTEGER, INTENT(IN) :: TASK

REAL(DKIND) Y(SIZE(X))

! Local variables:

INTEGER J

INTEGER, SAVE :: NX

REAL(DKIND), SAVE :: HSQ

SELECT CASE(TASK)

CASE(ARPACK_A_x)

! Computes y <-- A*x, where A is the N**2 by N**2 block

! tridiagonal matrix

! | T -I |

! |-I T -I |

! A = | -I T |

! | ... -I|

! | -I T|

Y(1:NX)=T(NX,X(1:NX)) - X(NX+1:2*NX)

DO J=NX+1,NX**2-NX,NX

Y(J:J+NX-1)=T(NX,X(J:J+NX-1))&

- X(J-NX:J-1)-X(J+NX:J+2*NX-1)

END DO

Y((NX-1)*NX+1:NX**2)= - X((NX-1)*NX-NX+1:(NX-1)*NX)&

+ T(NX,X((NX-1)*NX+1:NX**2))

! Note that HSQ is passed as a component of the extended type.

Y=(1._DKIND/HSQ)*Y

CASE(ARPACK_Prepare)

! Define NX, 1/H**2 so they are later available in the evaluator.

NX=10 ! This value is fixed in the evaluator.

HSQ = 1._DKIND/REAL(NX+1,DKIND)**2

Y=0._DKIND

CASE DEFAULT

WRITE(NOUT,*) TASK, ': INVALID TASK REQUESTED'

STOP 'IMSL_ERROR_WRONG_OPERATION'

END SELECT

CONTAINS

FUNCTION T(NX, X)RESULT(V)

INTEGER, INTENT(IN) :: NX

REAL(DKIND), INTENT(IN) :: X(:)

REAL(DKIND) :: V(NX)

REAL(DKIND) :: MONE=-1._DKIND, FOUR=4._DKIND

INTEGER J

V(1)=FOUR*X(1)+MONE*X(2)

DO J=2,NX-1

V(J)=MONE*X(J-1)+FOUR*X(J)+MONE*X(J+1)

END DO

V(NX)=MONE*X(NX-1)+FOUR*X(NX)

END FUNCTION

Output

Number of eigenvalues requested, and declared accurate

------------------------------------------------------

5 0

Number of Matrix-Vector Products Recorded, EX-11

------------------------------------------------

Smallest Laplacian Eigenvalues

1 19.61

2 48.22

3 48.22

4 76.83

5 93.33

All Ritz Values and Vectors have small residuals.

First 2D Laplacian Eigenfunction at Grid Points

1 2 3 4 5 6 7 8

1 0.0144 0.0277 0.0387 0.0466 0.0507 0.0507 0.0466 0.0387

2 0.0277 0.0531 0.0743 0.0894 0.0973 0.0973 0.0894 0.0743

3 0.0387 0.0743 0.1038 0.1250 0.1360 0.1360 0.1250 0.1038

4 0.0466 0.0894 0.1250 0.1504 0.1637 0.1637 0.1504 0.1250

5 0.0507 0.0973 0.1360 0.1637 0.1781 0.1781 0.1637 0.1360

6 0.0507 0.0973 0.1360 0.1637 0.1781 0.1781 0.1637 0.1360

7 0.0466 0.0894 0.1250 0.1504 0.1637 0.1637 0.1504 0.1250

8 0.0387 0.0743 0.1038 0.1250 0.1360 0.1360 0.1250 0.1038

9 0.0277 0.0531 0.0743 0.0894 0.0973 0.0973 0.0894 0.0743

10 0.0144 0.0277 0.0387 0.0466 0.0507 0.0507 0.0466 0.0387

9 10

1 0.0277 0.0144

2 0.0531 0.0277

3 0.0743 0.0387

4 0.0894 0.0466

5 0.0973 0.0507

6 0.0973 0.0507

7 0.0894 0.0466

8 0.0743 0.0387

9 0.0531 0.0277

10 0.0277 0.0144

Example 2

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)-1 x = 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.

IF (ALLOCATED(EX%A)) DEALLOCATE(EX%A)

IF (ALLOCATED(EX%B)) DEALLOCATE(EX%B)

IF (ALLOCATED(EX%C)) DEALLOCATE(EX%C)

IF (ALLOCATED(EX%Y1)) DEALLOCATE(EX%Y1)

IF (ALLOCATED(EX%U)) DEALLOCATE(EX%U)

IF (ALLOCATED(EX%IR)) DEALLOCATE(EX%IR)

IF (ALLOCATED(EX%IS)) DEALLOCATE(EX%IS)

NSIZE = (log(dble(N))/log(2.0)) + 5

ALLOCATE(EX%A(2*N), EX%B(2*N), EX%C(2*N), EX%Y1(2*N), &

EX%U(2*N), EX%IR(NSIZE), &

EX%IS(NSIZE), STAT=IERR)

IF (IERR /= 0) STOP 'IMSL_ERROR_ALLOCATING_WORKSPACE'

! Define matrix values.

HSQ=1._DKIND/REAL((N+1)**2,DKIND)

EX%B(1:N) = -1._DKIND/HSQ

EX%A(1:N) = 2._DKIND/HSQ - SHIFT

EX%C(1:N) = EX%B(1:N)

EX%Y1(:) = 0.0_DKIND

! Factor the matrix with LU and partial pivoting.

IJOB = 3

CALL LSLCR (EX%C, EX%A, EX%B, EX%Y1, EX%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 values to satisfy compiler.

Y=0._DKIND

NV=0

CASE DEFAULT

STOP 'IMSL_ERROR_WRONG_OPERATION'

END SELECT

END ASSOCIATE

END SELECT

END FUNCTION

END MODULE

USE ARPACK_SYMMETRIC_EX2_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 is the 1 dimensional discrete Laplacian on

! the interval 0,1 with zero Dirichlet boundary condition.

INTEGER, PARAMETER :: NEV=4, N=100

REAL(DKIND) :: VALUES(NEV), RES(N)

REAL(DKIND), ALLOCATABLE :: VECTORS(:,:)

REAL(DKIND) NORM

LOGICAL SMALL, SOLVED

INTEGER J, NOUT

TYPE(ARPACKBASE_EXT) EX

ASSOCIATE(NX => EX%NX, &

NV => EX%NV, &

SIGMA => EX%SHIFT)

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.

SIGMA=0._DKIND

CALL ARPACK_SYMMETRIC(N, FCN, VALUES,&

CATEGORY=ARPACK_Shift_Invert, EXTYPE=EX, VECTORS=VECTORS)

WRITE(NOUT,*) 'Number of Matrix-Vector Products Required, EX-2'

WRITE(NOUT,*) '-----------------------------------------------'

WRITE(NOUT, '(5X, I4)') NV

CALL WRRRN('Largest Laplacian Eigenvalues Near Zero Shift', &

VALUES)

! Check residuals, A*vectors = values*vectors:

DO J=1,NEV

RES=FCN(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.'

ENDIF

END ASSOCIATE

END

Output

Number of Matrix-Vector Products Required, EX-2

-----------------------------------------------

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.

CALL ARPACK_SYMMETRIC(N, FCN, VALUES, &

TYPE=ARPACK_Generalized, &

CATEGORY=ARPACK_Regular_Inverse, &

PLACE=ARPACK_Smallest_Magnitude, EXTYPE=EX, VECTORS=VECTORS)

WRITE(NOUT,*) 'Number of Matrix-Vector Products Required, EX-3'

WRITE(NOUT,*) '-----------------------------------------------'

WRITE(NOUT, '(5X, I4)') NV

CALL WRRRN('Largest Laplacian Eigenvalues', VALUES)

! Check residuals, A*vectors = values*B*vectors:

DO J=1,NEV

RES=FCN(VECTORS(:,J),ARPACK_A_x,EX)-&

VALUES(J)*FCN(VECTORS(:,J),ARPACK_B_x,EX)

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.'

ENDIF

END ASSOCIATE

END PROGRAM

FUNCTION FCN(X, TASK, EX) RESULT(Y)

USE ARPACK_SYMMETRIC_EX3_INT

CLASS (ARPACKBASE), INTENT(INOUT) :: EX

REAL (DKIND), INTENT(INOUT) :: X(:)

INTEGER, INTENT(IN) :: TASK

REAL (DKIND) Y(SIZE(X)), PI

INTEGER J, IERR, IJOB, NSIZE

SELECT TYPE(EX)

TYPE IS(ARPACKBASE_EXT)

ASSOCIATE(N => EX%NCOLS, &

NV => EX%NV, &

H => EX%H)

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/H

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)

Y=Y*(H/6._DKIND)

CASE(ARPACK_inv_B_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.

IF (ALLOCATED(EX%A)) DEALLOCATE(EX%A)

IF (ALLOCATED(EX%B)) DEALLOCATE(EX%B)

IF (ALLOCATED(EX%C)) DEALLOCATE(EX%C)

IF (ALLOCATED(EX%Y1)) DEALLOCATE(EX%Y1)

IF (ALLOCATED(EX%U)) DEALLOCATE(EX%U)

IF (ALLOCATED(EX%IR)) DEALLOCATE(EX%IR)

IF (ALLOCATED(EX%IS)) DEALLOCATE(EX%IS)

NSIZE = (log(dble(N))/log(2.0d0)) + 5

ALLOCATE(EX%A(2*N), EX%B(2*N), EX%C(2*N), EX%Y1(2*N), &

EX%U(2*N), EX%IR(NSIZE), &

EX%IS(NSIZE), STAT=IERR)

IF (IERR /= 0) STOP 'IMSL_ERROR_ALLOCATING_WORKSPACE'

! Define matrix values.

PI=ATAN(1._DKIND)*4._DKIND

H=PI/REAL(N+1,DKIND)

EX%B(1:N) = (1._DKIND/6._DKIND)*H

EX%A(1:N) = (2._DKIND/3._DKIND)*H

EX%C(1:N) = EX%B(1:N)

EX%Y1(:) = 0.0_DKIND

! Factor the matrix with LU and partial pivoting.

IJOB = 3

CALL LSLCR (EX%C, EX%A, EX%B, EX%Y1, EX%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 values to satisfy compiler.

Y=0._DKIND

NV=0

CASE DEFAULT

write(*,*)TASK

STOP 'IMSL_ERROR_WRONG_OPERATION'

END SELECT

END ASSOCIATE

END SELECT

END FUNCTION

Output

Number of Matrix-Vector Products Required, EX-3

-----------------------------------------------

1126

Largest Laplacian Eigenvalues

1 1.00

2 4.00

3 9.01

4 16.02

All Ritz Values and Vectors have small residuals.