Chapter 2: Eigensystem Analysis

ARPACK_NONSYMMETRIC

Compute some eigenvalues and eigenvectors of the generalized eigenvalue problem.  This can be used for the case.   The values forare real, but eigenvalues may be complex and occur in 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 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.

Optional Arguments

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.

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

FORTRAN 2003 Interface

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

Specific:         The specific interface name is D_ARPACK_NONSYMMETRIC.

FORTRAN 90 Interface

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

Description

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.

Comments

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 vectorsare consecutive columns of the array VECTORS (:,:).   The eigenvalue-eigenvector relationship is.  Sinceis real, is also an eigenvalue; thus the conjugate relationshipwill hold.  For purposes of checking results the complex residualshould 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 solvingfor: .   is valid if the denominator is non-zero.  Ifhas a non-zero imaginary part, then the complex conjugateis 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.

Example 1

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

Output

 

 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.



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