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