Chapter 2: Eigensystem Analysis

ARPACK_COMPLEX

Compute some eigenvalues and eigenvectors of the generalized eigenvalue problem.  This can be used for the case.  The values forare real or complex.  When the values are complex the eigenvalues may be complex and are not expected to occur in complex conjugate pairs.

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

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_Magnitude

ARPACK_Smallest_Magnitude

ARPACK_Largest_Real_Parts

ARPACK_Smallest_Real_Parts

ARPACK_Largest_Imag_Parts

ARPACK_Smallest_Imag_Parts

 

            Default: PLACEARPACK_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).

FORTRAN 2003 Interface

Generic:          ARPACK_COMPLEX (N, F,ZVALUES [,…])

Specific:         The specific interface name is Z_ARPACK_COMPLEX.

FORTRAN 90 Interface

A Fortran 90 compiler version is not available for this routine.

Description

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.

Comments

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 solvingfor,

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.

Example 1

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

Output

 

 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.



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