LAPACK: TESTING/EIG/cget52.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA,
00002      $                   WORK, RWORK, RESULT )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       LOGICAL            LEFT
00010       INTEGER            LDA, LDB, LDE, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       REAL               RESULT( 2 ), RWORK( * )
00014       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
00015      $                   BETA( * ), E( LDE, * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CGET52  does an eigenvector check for the generalized eigenvalue
00022 *  problem.
00023 *
00024 *  The basic test for right eigenvectors is:
00025 *
00026 *                            | b(i) A E(i) -  a(i) B E(i) |
00027 *          RESULT(1) = max   -------------------------------
00028 *                       i    n ulp max( |b(i) A|, |a(i) B| )
00029 *
00030 *  using the 1-norm.  Here, a(i)/b(i) = w is the i-th generalized
00031 *  eigenvalue of A - w B, or, equivalently, b(i)/a(i) = m is the i-th
00032 *  generalized eigenvalue of m A - B.
00033 *
00034 *                          H   H  _      _
00035 *  For left eigenvectors, A , B , a, and b  are used.
00036 *
00037 *  CGET52 also tests the normalization of E.  Each eigenvector is
00038 *  supposed to be normalized so that the maximum "absolute value"
00039 *  of its elements is 1, where in this case, "absolute value"
00040 *  of a complex value x is  |Re(x)| + |Im(x)| ; let us call this
00041 *  maximum "absolute value" norm of a vector v  M(v).  
00042 *  if a(i)=b(i)=0, then the eigenvector is set to be the jth coordinate
00043 *  vector. The normalization test is:
00044 *
00045 *          RESULT(2) =      max       | M(v(i)) - 1 | / ( n ulp )
00046 *                     eigenvectors v(i)
00047 *
00048 *  Arguments
00049 *  =========
00050 *
00051 *  LEFT    (input) LOGICAL
00052 *          =.TRUE.:  The eigenvectors in the columns of E are assumed
00053 *                    to be *left* eigenvectors.
00054 *          =.FALSE.: The eigenvectors in the columns of E are assumed
00055 *                    to be *right* eigenvectors.
00056 *
00057 *  N       (input) INTEGER
00058 *          The size of the matrices.  If it is zero, CGET52 does
00059 *          nothing.  It must be at least zero.
00060 *
00061 *  A       (input) COMPLEX array, dimension (LDA, N)
00062 *          The matrix A.
00063 *
00064 *  LDA     (input) INTEGER
00065 *          The leading dimension of A.  It must be at least 1
00066 *          and at least N.
00067 *
00068 *  B       (input) COMPLEX array, dimension (LDB, N)
00069 *          The matrix B.
00070 *
00071 *  LDB     (input) INTEGER
00072 *          The leading dimension of B.  It must be at least 1
00073 *          and at least N.
00074 *
00075 *  E       (input) COMPLEX array, dimension (LDE, N)
00076 *          The matrix of eigenvectors.  It must be O( 1 ).
00077 *
00078 *  LDE     (input) INTEGER
00079 *          The leading dimension of E.  It must be at least 1 and at
00080 *          least N.
00081 *
00082 *  ALPHA   (input) COMPLEX array, dimension (N)
00083 *          The values a(i) as described above, which, along with b(i),
00084 *          define the generalized eigenvalues.
00085 *
00086 *  BETA    (input) COMPLEX array, dimension (N)
00087 *          The values b(i) as described above, which, along with a(i),
00088 *          define the generalized eigenvalues.
00089 *
00090 *  WORK    (workspace) COMPLEX array, dimension (N**2)
00091 *
00092 *  RWORK   (workspace) REAL array, dimension (N)
00093 *
00094 *  RESULT  (output) REAL array, dimension (2)
00095 *          The values computed by the test described above.  If A E or
00096 *          B E is likely to overflow, then RESULT(1:2) is set to
00097 *          10 / ulp.
00098 *
00099 *  =====================================================================
00100 *
00101 *     .. Parameters ..
00102       REAL               ZERO, ONE
00103       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00104       COMPLEX            CZERO, CONE
00105       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00106      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00107 *     ..
00108 *     .. Local Scalars ..
00109       CHARACTER          NORMAB, TRANS
00110       INTEGER            J, JVEC
00111       REAL               ABMAX, ALFMAX, ANORM, BETMAX, BNORM, ENORM,
00112      $                   ENRMER, ERRNRM, SAFMAX, SAFMIN, SCALE, TEMP1,
00113      $                   ULP
00114       COMPLEX            ACOEFF, ALPHAI, BCOEFF, BETAI, X
00115 *     ..
00116 *     .. External Functions ..
00117       REAL               CLANGE, SLAMCH
00118       EXTERNAL           CLANGE, SLAMCH
00119 *     ..
00120 *     .. External Subroutines ..
00121       EXTERNAL           CGEMV
00122 *     ..
00123 *     .. Intrinsic Functions ..
00124       INTRINSIC          ABS, AIMAG, CONJG, MAX, REAL
00125 *     ..
00126 *     .. Statement Functions ..
00127       REAL               ABS1
00128 *     ..
00129 *     .. Statement Function definitions ..
00130       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
00131 *     ..
00132 *     .. Executable Statements ..
00133 *
00134       RESULT( 1 ) = ZERO
00135       RESULT( 2 ) = ZERO
00136       IF( N.LE.0 )
00137      $   RETURN
00138 *
00139       SAFMIN = SLAMCH( 'Safe minimum' )
00140       SAFMAX = ONE / SAFMIN
00141       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00142 *
00143       IF( LEFT ) THEN
00144          TRANS = 'C'
00145          NORMAB = 'I'
00146       ELSE
00147          TRANS = 'N'
00148          NORMAB = 'O'
00149       END IF
00150 *
00151 *     Norm of A, B, and E:
00152 *
00153       ANORM = MAX( CLANGE( NORMAB, N, N, A, LDA, RWORK ), SAFMIN )
00154       BNORM = MAX( CLANGE( NORMAB, N, N, B, LDB, RWORK ), SAFMIN )
00155       ENORM = MAX( CLANGE( 'O', N, N, E, LDE, RWORK ), ULP )
00156       ALFMAX = SAFMAX / MAX( ONE, BNORM )
00157       BETMAX = SAFMAX / MAX( ONE, ANORM )
00158 *
00159 *     Compute error matrix.
00160 *     Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
00161 *
00162       DO 10 JVEC = 1, N
00163          ALPHAI = ALPHA( JVEC )
00164          BETAI = BETA( JVEC )
00165          ABMAX = MAX( ABS1( ALPHAI ), ABS1( BETAI ) )
00166          IF( ABS1( ALPHAI ).GT.ALFMAX .OR. ABS1( BETAI ).GT.BETMAX .OR.
00167      $       ABMAX.LT.ONE ) THEN
00168             SCALE = ONE / MAX( ABMAX, SAFMIN )
00169             ALPHAI = SCALE*ALPHAI
00170             BETAI = SCALE*BETAI
00171          END IF
00172          SCALE = ONE / MAX( ABS1( ALPHAI )*BNORM, ABS1( BETAI )*ANORM,
00173      $           SAFMIN )
00174          ACOEFF = SCALE*BETAI
00175          BCOEFF = SCALE*ALPHAI
00176          IF( LEFT ) THEN
00177             ACOEFF = CONJG( ACOEFF )
00178             BCOEFF = CONJG( BCOEFF )
00179          END IF
00180          CALL CGEMV( TRANS, N, N, ACOEFF, A, LDA, E( 1, JVEC ), 1,
00181      $               CZERO, WORK( N*( JVEC-1 )+1 ), 1 )
00182          CALL CGEMV( TRANS, N, N, -BCOEFF, B, LDA, E( 1, JVEC ), 1,
00183      $               CONE, WORK( N*( JVEC-1 )+1 ), 1 )
00184    10 CONTINUE
00185 *
00186       ERRNRM = CLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM
00187 *
00188 *     Compute RESULT(1)
00189 *
00190       RESULT( 1 ) = ERRNRM / ULP
00191 *
00192 *     Normalization of E:
00193 *
00194       ENRMER = ZERO
00195       DO 30 JVEC = 1, N
00196          TEMP1 = ZERO
00197          DO 20 J = 1, N
00198             TEMP1 = MAX( TEMP1, ABS1( E( J, JVEC ) ) )
00199    20    CONTINUE
00200          ENRMER = MAX( ENRMER, TEMP1-ONE )
00201    30 CONTINUE
00202 *
00203 *     Compute RESULT(2) : the normalization error in E.
00204 *
00205       RESULT( 2 ) = ENRMER / ( REAL( N )*ULP )
00206 *
00207       RETURN
00208 *
00209 *     End of CGET52
00210 *
00211       END