Chapter 2: Eigensystem Analysis

ARPACK_SYMMETRIC

Computes some eigenvalues and eigenvectors of the generalized real symmetric eigenvalue problem.  This can be used for the case.

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

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_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

ARPACK_Generalized

 

            Default: TYPE = ARPACK_Standard.

CATEGORYCATEGORY 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_Buckling

ARPACK_Cayley

 

            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_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

Example 1

We approximate eigenvalues and eigenfunctions of the Laplacian operator

,

on the unit square, , with zero Dirichlet boundary values.  The full set of eigenvalues and their eigenfunctions are given by the sequence, whereare 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. 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.

Link to example source (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   

      END FUNCTION

Output

 

Number of eigenvalues requested, and declared accurate

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

        5        0

 Number of Matrix-Vector Products Recorded, EX-11

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

        0

 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 operatoron the unit interval, .  Equally spaced divided differences are used for the operator, which yields a tri-diagonal matrix.  The eigenvalues and eigenfunctions are . 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 is computed, here with.  When FCN is later entered with TASK=ARPACK_inv_A_minus_Shift_x, the LU factorization is available to efficiently compute.  The function FCN is also entered with TASK=ARPACK_A_x, to compute .

Link to example source (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

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

       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 operatoron the unit interval, .  This is Example 2 but solved using finite elements.  With matrix notation, we have the matrix problem.  Both andare tri-diagonal and symmetric.  The matrix 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,.  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.  When FCN is entered with TASK=ARPACK_Prepare, these arrays are allocated, defined, and the LU factorization of the matrix is computed. The function FCN is entered with the three values TASK=ARPACK_A_x, for; TASK=ARPACK_B_x, for; and TASK=ARPACK_inv_B_x, for .

Within the main program, the function interface for the external function FCN is specified with the declaration PROCEDURE (DMV) FCN.

Link to example source (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.



http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260