1
       2
       3
       4
       5
       6
       7
       8
       9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      47
      48
      49
      50
      51
      52
      53
      54
      55
      56
      57
      58
      59
      60
      61
      62
      63
      64
      65
      66
      67
      68
      69
      70
      71
      72
      73
      74
      75
      76
      77
      78
      79
      80
      81
      82
      83
      84
      85
      86
      87
      88
      89
      90
      91
      92
      93
      94
      95
      96
      97
      98
      99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112
     113
     114
     115
     116
     117
     118
     119
     120
     121
     122
     123
     124
     125
     126
     127
     128
     129
     130
     131
     132
     133
     134
     135
     136
     137
     138
     139
     140
     141
     142
     143
     144
     145
     146
     147
     148
     149
     150
     151
     152
     153
     154
     155
     156
     157
     158
     159
     160
     161
     162
     163
     164
     165
     166
     167
     168
     169
     170
     171
     172
     173
     174
     175
     176
     177
     178
     179
     180
     181
     182
     183
     184
     185
     186
     187
     188
     189
     190
     191
     192
     193
     194
     195
     196
     197
     198
     199
     200
     201
     202
     203
     204
     205
     206
     207
     208
     209
     210
     211
      SUBROUTINE CGET52LEFTNALDABLDBELDEALPHABETA,
     $                   WORKRWORKRESULT )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      LOGICAL            LEFT
      INTEGER            LDALDBLDEN
*     ..
*     .. Array Arguments ..
      REAL               RESULT2 ), RWORK* )
      COMPLEX            ALDA* ), ALPHA* ), BLDB* ),
     $                   BETA* ), ELDE* ), WORK* )
*     ..
*
*  Purpose
*  =======
*
*  CGET52  does an eigenvector check for the generalized eigenvalue
*  problem.
*
*  The basic test for right eigenvectors is:
*
*                            | b(i) A E(i) -  a(i) B E(i) |
*          RESULT(1) = max   -------------------------------
*                       i    n ulp max( |b(i) A|, |a(i) B| )
*
*  using the 1-norm.  Here, a(i)/b(i) = w is the i-th generalized
*  eigenvalue of A - w B, or, equivalently, b(i)/a(i) = m is the i-th
*  generalized eigenvalue of m A - B.
*
*                          H   H  _      _
*  For left eigenvectors, A , B , a, and b  are used.
*
*  CGET52 also tests the normalization of E.  Each eigenvector is
*  supposed to be normalized so that the maximum "absolute value"
*  of its elements is 1, where in this case, "absolute value"
*  of a complex value x is  |Re(x)| + |Im(x)| ; let us call this
*  maximum "absolute value" norm of a vector v  M(v).  
*  if a(i)=b(i)=0, then the eigenvector is set to be the jth coordinate
*  vector. The normalization test is:
*
*          RESULT(2) =      max       | M(v(i)) - 1 | / ( n ulp )
*                     eigenvectors v(i)
*
*  Arguments
*  =========
*
*  LEFT    (input) LOGICAL
*          =.TRUE.:  The eigenvectors in the columns of E are assumed
*                    to be *left* eigenvectors.
*          =.FALSE.: The eigenvectors in the columns of E are assumed
*                    to be *right* eigenvectors.
*
*  N       (input) INTEGER
*          The size of the matrices.  If it is zero, CGET52 does
*          nothing.  It must be at least zero.
*
*  A       (input) COMPLEX array, dimension (LDA, N)
*          The matrix A.
*
*  LDA     (input) INTEGER
*          The leading dimension of A.  It must be at least 1
*          and at least N.
*
*  B       (input) COMPLEX array, dimension (LDB, N)
*          The matrix B.
*
*  LDB     (input) INTEGER
*          The leading dimension of B.  It must be at least 1
*          and at least N.
*
*  E       (input) COMPLEX array, dimension (LDE, N)
*          The matrix of eigenvectors.  It must be O( 1 ).
*
*  LDE     (input) INTEGER
*          The leading dimension of E.  It must be at least 1 and at
*          least N.
*
*  ALPHA   (input) COMPLEX array, dimension (N)
*          The values a(i) as described above, which, along with b(i),
*          define the generalized eigenvalues.
*
*  BETA    (input) COMPLEX array, dimension (N)
*          The values b(i) as described above, which, along with a(i),
*          define the generalized eigenvalues.
*
*  WORK    (workspace) COMPLEX array, dimension (N**2)
*
*  RWORK   (workspace) REAL array, dimension (N)
*
*  RESULT  (output) REAL array, dimension (2)
*          The values computed by the test described above.  If A E or
*          B E is likely to overflow, then RESULT(1:2) is set to
*          10 / ulp.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZEROONE
      PARAMETER          ( ZERO = 0.0E+0ONE = 1.0E+0 )
      COMPLEX            CZEROCONE
      PARAMETER          ( CZERO = ( 0.0E+00.0E+0 ),
     $                   CONE = ( 1.0E+00.0E+0 ) )
*     ..
*     .. Local Scalars ..
      CHARACTER          NORMABTRANS
      INTEGER            JJVEC
      REAL               ABMAXALFMAXANORMBETMAXBNORMENORM,
     $                   ENRMERERRNRMSAFMAXSAFMINSCALETEMP1,
     $                   ULP
      COMPLEX            ACOEFFALPHAIBCOEFFBETAIX
*     ..
*     .. External Functions ..
      REAL               CLANGESLAMCH
      EXTERNAL           CLANGESLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           CGEMV
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABSAIMAGCONJGMAXREAL
*     ..
*     .. Statement Functions ..
      REAL               ABS1
*     ..
*     .. Statement Function definitions ..
      ABS1X ) = ABSREALX ) ) + ABSAIMAGX ) )
*     ..
*     .. Executable Statements ..
*
      RESULT1 ) = ZERO
      RESULT2 ) = ZERO
      IFN.LE.0 )
     $   RETURN
*
      SAFMIN = SLAMCH'Safe minimum' )
      SAFMAX = ONE / SAFMIN
      ULP = SLAMCH'Epsilon' )*SLAMCH'Base' )
*
      IFLEFT ) THEN
         TRANS = 'C'
         NORMAB = 'I'
      ELSE
         TRANS = 'N'
         NORMAB = 'O'
      END IF
*
*     Norm of A, B, and E:
*
      ANORM = MAXCLANGENORMABNNALDARWORK ), SAFMIN )
      BNORM = MAXCLANGENORMABNNBLDBRWORK ), SAFMIN )
      ENORM = MAXCLANGE'O'NNELDERWORK ), ULP )
      ALFMAX = SAFMAX / MAXONEBNORM )
      BETMAX = SAFMAX / MAXONEANORM )
*
*     Compute error matrix.
*     Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
*
      DO 10 JVEC = 1N
         ALPHAI = ALPHAJVEC )
         BETAI = BETAJVEC )
         ABMAX = MAXABS1ALPHAI ), ABS1BETAI ) )
         IFABS1ALPHAI ).GT.ALFMAX .OR. ABS1BETAI ).GT.BETMAX .OR.
     $       ABMAX.LT.ONE ) THEN
            SCALE = ONE / MAXABMAXSAFMIN )
            ALPHAI = SCALE*ALPHAI
            BETAI = SCALE*BETAI
         END IF
         SCALE = ONE / MAXABS1ALPHAI )*BNORMABS1BETAI )*ANORM,
     $           SAFMIN )
         ACOEFF = SCALE*BETAI
         BCOEFF = SCALE*ALPHAI
         IFLEFT ) THEN
            ACOEFF = CONJGACOEFF )
            BCOEFF = CONJGBCOEFF )
         END IF
         CALL CGEMVTRANSNNACOEFFALDAE1JVEC ), 1,
     $               CZEROWORKN*JVEC-1 )+1 ), 1 )
         CALL CGEMVTRANSNN-BCOEFFBLDAE1JVEC ), 1,
     $               CONEWORKN*JVEC-1 )+1 ), 1 )
   10 CONTINUE
*
      ERRNRM = CLANGE'One'NNWORKNRWORK ) / ENORM
*
*     Compute RESULT(1)
*
      RESULT1 ) = ERRNRM / ULP
*
*     Normalization of E:
*
      ENRMER = ZERO
      DO 30 JVEC = 1N
         TEMP1 = ZERO
         DO 20 J = 1N
            TEMP1 = MAXTEMP1ABS1EJJVEC ) ) )
   20    CONTINUE
         ENRMER = MAXENRMERTEMP1-ONE )
   30 CONTINUE
*
*     Compute RESULT(2) : the normalization error in E.
*
      RESULT2 ) = ENRMER / ( REALN )*ULP )
*
      RETURN
*
*     End of CGET52
*
      END