1       SUBROUTINE SGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR,
  2      $                   ALPHAI, BETA, WORK, RESULT )
  3 *
  4 *  -- LAPACK test routine (version 3.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     November 2006
  7 *
  8 *     .. Scalar Arguments ..
  9       LOGICAL            LEFT
 10       INTEGER            LDA, LDB, LDE, N
 11 *     ..
 12 *     .. Array Arguments ..
 13       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
 14      $                   B( LDB, * ), BETA( * ), E( LDE, * ),
 15      $                   RESULT2 ), WORK( * )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  SGET52  does an eigenvector check for the generalized eigenvalue
 22 *  problem.
 23 *
 24 *  The basic test for right eigenvectors is:
 25 *
 26 *                            | b(j) A E(j) -  a(j) B E(j) |
 27 *          RESULT(1) = max   -------------------------------
 28 *                       j    n ulp max( |b(j) A|, |a(j) B| )
 29 *
 30 *  using the 1-norm.  Here, a(j)/b(j) = w is the j-th generalized
 31 *  eigenvalue of A - w B, or, equivalently, b(j)/a(j) = m is the j-th
 32 *  generalized eigenvalue of m A - B.
 33 *
 34 *  For real eigenvalues, the test is straightforward.  For complex
 35 *  eigenvalues, E(j) and a(j) are complex, represented by
 36 *  Er(j) + i*Ei(j) and ar(j) + i*ai(j), resp., so the test for that
 37 *  eigenvector becomes
 38 *
 39 *                  max( |Wr|, |Wi| )
 40 *      --------------------------------------------
 41 *      n ulp max( |b(j) A|, (|ar(j)|+|ai(j)|) |B| )
 42 *
 43 *  where
 44 *
 45 *      Wr = b(j) A Er(j) - ar(j) B Er(j) + ai(j) B Ei(j)
 46 *
 47 *      Wi = b(j) A Ei(j) - ai(j) B Er(j) - ar(j) B Ei(j)
 48 *
 49 *                          T   T  _
 50 *  For left eigenvectors, A , B , a, and b  are used.
 51 *
 52 *  SGET52 also tests the normalization of E.  Each eigenvector is
 53 *  supposed to be normalized so that the maximum "absolute value"
 54 *  of its elements is 1, where in this case, "absolute value"
 55 *  of a complex value x is  |Re(x)| + |Im(x)| ; let us call this
 56 *  maximum "absolute value" norm of a vector v  M(v). 
 57 *  if a(j)=b(j)=0, then the eigenvector is set to be the jth coordinate
 58 *  vector.  The normalization test is:
 59 *
 60 *          RESULT(2) =      max       | M(v(j)) - 1 | / ( n ulp )
 61 *                     eigenvectors v(j)
 62 *
 63 *  Arguments
 64 *  =========
 65 *
 66 *  LEFT    (input) LOGICAL
 67 *          =.TRUE.:  The eigenvectors in the columns of E are assumed
 68 *                    to be *left* eigenvectors.
 69 *          =.FALSE.: The eigenvectors in the columns of E are assumed
 70 *                    to be *right* eigenvectors.
 71 *
 72 *  N       (input) INTEGER
 73 *          The size of the matrices.  If it is zero, SGET52 does
 74 *          nothing.  It must be at least zero.
 75 *
 76 *  A       (input) REAL array, dimension (LDA, N)
 77 *          The matrix A.
 78 *
 79 *  LDA     (input) INTEGER
 80 *          The leading dimension of A.  It must be at least 1
 81 *          and at least N.
 82 *
 83 *  B       (input) REAL array, dimension (LDB, N)
 84 *          The matrix B.
 85 *
 86 *  LDB     (input) INTEGER
 87 *          The leading dimension of B.  It must be at least 1
 88 *          and at least N.
 89 *
 90 *  E       (input) REAL array, dimension (LDE, N)
 91 *          The matrix of eigenvectors.  It must be O( 1 ).  Complex
 92 *          eigenvalues and eigenvectors always come in pairs, the
 93 *          eigenvalue and its conjugate being stored in adjacent
 94 *          elements of ALPHAR, ALPHAI, and BETA.  Thus, if a(j)/b(j)
 95 *          and a(j+1)/b(j+1) are a complex conjugate pair of
 96 *          generalized eigenvalues, then E(,j) contains the real part
 97 *          of the eigenvector and E(,j+1) contains the imaginary part.
 98 *          Note that whether E(,j) is a real eigenvector or part of a
 99 *          complex one is specified by whether ALPHAI(j) is zero or not.
100 *
101 *  LDE     (input) INTEGER
102 *          The leading dimension of E.  It must be at least 1 and at
103 *          least N.
104 *
105 *  ALPHAR  (input) REAL array, dimension (N)
106 *          The real parts of the values a(j) as described above, which,
107 *          along with b(j), define the generalized eigenvalues.
108 *          Complex eigenvalues always come in complex conjugate pairs
109 *          a(j)/b(j) and a(j+1)/b(j+1), which are stored in adjacent
110 *          elements in ALPHAR, ALPHAI, and BETA.  Thus, if the j-th
111 *          and (j+1)-st eigenvalues form a pair, ALPHAR(j+1)/BETA(j+1)
112 *          is assumed to be equal to ALPHAR(j)/BETA(j).
113 *
114 *  ALPHAI  (input) REAL array, dimension (N)
115 *          The imaginary parts of the values a(j) as described above,
116 *          which, along with b(j), define the generalized eigenvalues.
117 *          If ALPHAI(j)=0, then the eigenvalue is real, otherwise it
118 *          is part of a complex conjugate pair.  Complex eigenvalues
119 *          always come in complex conjugate pairs a(j)/b(j) and
120 *          a(j+1)/b(j+1), which are stored in adjacent elements in
121 *          ALPHAR, ALPHAI, and BETA.  Thus, if the j-th and (j+1)-st
122 *          eigenvalues form a pair, ALPHAI(j+1)/BETA(j+1) is assumed to
123 *          be equal to  -ALPHAI(j)/BETA(j).  Also, nonzero values in
124 *          ALPHAI are assumed to always come in adjacent pairs.
125 *
126 *  BETA    (input) REAL array, dimension (N)
127 *          The values b(j) as described above, which, along with a(j),
128 *          define the generalized eigenvalues.
129 *
130 *  WORK    (workspace) REAL array, dimension (N**2+N)
131 *
132 *  RESULT  (output) REAL array, dimension (2)
133 *          The values computed by the test described above.  If A E or
134 *          B E is likely to overflow, then RESULT(1:2) is set to
135 *          10 / ulp.
136 *
137 *  =====================================================================
138 *
139 *     .. Parameters ..
140       REAL               ZERO, ONE, TEN
141       PARAMETER          ( ZERO = 0.0, ONE = 1.0, TEN = 10.0 )
142 *     ..
143 *     .. Local Scalars ..
144       LOGICAL            ILCPLX
145       CHARACTER          NORMAB, TRANS
146       INTEGER            J, JVEC
147       REAL               ABMAX, ACOEF, ALFMAX, ANORM, BCOEFI, BCOEFR,
148      $                   BETMAX, BNORM, ENORM, ENRMER, ERRNRM, SAFMAX,
149      $                   SAFMIN, SALFI, SALFR, SBETA, SCALE, TEMP1, ULP
150 *     ..
151 *     .. External Functions ..
152       REAL               SLAMCH, SLANGE
153       EXTERNAL           SLAMCH, SLANGE
154 *     ..
155 *     .. External Subroutines ..
156       EXTERNAL           SGEMV
157 *     ..
158 *     .. Intrinsic Functions ..
159       INTRINSIC          ABSMAX, REAL
160 *     ..
161 *     .. Executable Statements ..
162 *
163       RESULT1 ) = ZERO
164       RESULT2 ) = ZERO
165       IF( N.LE.0 )
166      $   RETURN
167 *
168       SAFMIN = SLAMCH( 'Safe minimum' )
169       SAFMAX = ONE / SAFMIN
170       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
171 *
172       IF( LEFT ) THEN
173          TRANS = 'T'
174          NORMAB = 'I'
175       ELSE
176          TRANS = 'N'
177          NORMAB = 'O'
178       END IF
179 *
180 *     Norm of A, B, and E:
181 *
182       ANORM = MAX( SLANGE( NORMAB, N, N, A, LDA, WORK ), SAFMIN )
183       BNORM = MAX( SLANGE( NORMAB, N, N, B, LDB, WORK ), SAFMIN )
184       ENORM = MAX( SLANGE( 'O', N, N, E, LDE, WORK ), ULP )
185       ALFMAX = SAFMAX / MAX( ONE, BNORM )
186       BETMAX = SAFMAX / MAX( ONE, ANORM )
187 *
188 *     Compute error matrix.
189 *     Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
190 *
191       ILCPLX = .FALSE.
192       DO 10 JVEC = 1, N
193          IF( ILCPLX ) THEN
194 *
195 *           2nd Eigenvalue/-vector of pair -- do nothing
196 *
197             ILCPLX = .FALSE.
198          ELSE
199             SALFR = ALPHAR( JVEC )
200             SALFI = ALPHAI( JVEC )
201             SBETA = BETA( JVEC )
202             IF( SALFI.EQ.ZERO ) THEN
203 *
204 *              Real eigenvalue and -vector
205 *
206                ABMAX = MAXABS( SALFR ), ABS( SBETA ) )
207                IFABS( SALFR ).GT.ALFMAX .OR. ABS( SBETA ).GT.
208      $             BETMAX .OR. ABMAX.LT.ONE ) THEN
209                   SCALE = ONE / MAX( ABMAX, SAFMIN )
210                   SALFR = SCALE*SALFR
211                   SBETA = SCALE*SBETA
212                END IF
213                SCALE = ONE / MAXABS( SALFR )*BNORM,
214      $                 ABS( SBETA )*ANORM, SAFMIN )
215                ACOEF = SCALE*SBETA
216                BCOEFR = SCALE*SALFR
217                CALL SGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC ), 1,
218      $                     ZERO, WORK( N*( JVEC-1 )+1 ), 1 )
219                CALL SGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC ),
220      $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
221             ELSE
222 *
223 *              Complex conjugate pair
224 *
225                ILCPLX = .TRUE.
226                IF( JVEC.EQ.N ) THEN
227                   RESULT1 ) = TEN / ULP
228                   RETURN
229                END IF
230                ABMAX = MAXABS( SALFR )+ABS( SALFI ), ABS( SBETA ) )
231                IFABS( SALFR )+ABS( SALFI ).GT.ALFMAX .OR.
232      $             ABS( SBETA ).GT.BETMAX .OR. ABMAX.LT.ONE ) THEN
233                   SCALE = ONE / MAX( ABMAX, SAFMIN )
234                   SALFR = SCALE*SALFR
235                   SALFI = SCALE*SALFI
236                   SBETA = SCALE*SBETA
237                END IF
238                SCALE = ONE / MAX( ( ABS( SALFR )+ABS( SALFI ) )*BNORM,
239      $                 ABS( SBETA )*ANORM, SAFMIN )
240                ACOEF = SCALE*SBETA
241                BCOEFR = SCALE*SALFR
242                BCOEFI = SCALE*SALFI
243                IF( LEFT ) THEN
244                   BCOEFI = -BCOEFI
245                END IF
246 *
247                CALL SGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC ), 1,
248      $                     ZERO, WORK( N*( JVEC-1 )+1 ), 1 )
249                CALL SGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC ),
250      $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
251                CALL SGEMV( TRANS, N, N, BCOEFI, B, LDA, E( 1, JVEC+1 ),
252      $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
253 *
254                CALL SGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC+1 ),
255      $                     1, ZERO, WORK( N*JVEC+1 ), 1 )
256                CALL SGEMV( TRANS, N, N, -BCOEFI, B, LDA, E( 1, JVEC ),
257      $                     1, ONE, WORK( N*JVEC+1 ), 1 )
258                CALL SGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC+1 ),
259      $                     1, ONE, WORK( N*JVEC+1 ), 1 )
260             END IF
261          END IF
262    10 CONTINUE
263 *
264       ERRNRM = SLANGE( 'One', N, N, WORK, N, WORK( N**2+1 ) ) / ENORM
265 *
266 *     Compute RESULT(1)
267 *
268       RESULT1 ) = ERRNRM / ULP
269 *
270 *     Normalization of E:
271 *
272       ENRMER = ZERO
273       ILCPLX = .FALSE.
274       DO 40 JVEC = 1, N
275          IF( ILCPLX ) THEN
276             ILCPLX = .FALSE.
277          ELSE
278             TEMP1 = ZERO
279             IF( ALPHAI( JVEC ).EQ.ZERO ) THEN
280                DO 20 J = 1, N
281                   TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) )
282    20          CONTINUE
283                ENRMER = MAX( ENRMER, TEMP1-ONE )
284             ELSE
285                ILCPLX = .TRUE.
286                DO 30 J = 1, N
287                   TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+
288      $                    ABS( E( J, JVEC+1 ) ) )
289    30          CONTINUE
290                ENRMER = MAX( ENRMER, TEMP1-ONE )
291             END IF
292          END IF
293    40 CONTINUE
294 *
295 *     Compute RESULT(2) : the normalization error in E.
296 *
297       RESULT2 ) = ENRMER / ( REAL( N )*ULP )
298 *
299       RETURN
300 *
301 *     End of SGET52
302 *
303       END