1*> \brief \b SGET22
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 SGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR,
12*                          WI, WORK, RESULT )
13*
14*       .. Scalar Arguments ..
15*       CHARACTER          TRANSA, TRANSE, TRANSW
16*       INTEGER            LDA, LDE, N
17*       ..
18*       .. Array Arguments ..
19*       REAL               A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
20*      $                   WORK( * ), WR( * )
21*       ..
22*
23*
24*> \par Purpose:
25*  =============
26*>
27*> \verbatim
28*>
29*> SGET22 does an eigenvector check.
30*>
31*> The basic test is:
32*>
33*>    RESULT(1) = | A E  -  E W | / ( |A| |E| ulp )
34*>
35*> using the 1-norm.  It also tests the normalization of E:
36*>
37*>    RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
38*>                 j
39*>
40*> where E(j) is the j-th eigenvector, and m-norm is the max-norm of a
41*> vector.  If an eigenvector is complex, as determined from WI(j)
42*> nonzero, then the max-norm of the vector ( er + i*ei ) is the maximum
43*> of
44*>    |er(1)| + |ei(1)|, ... , |er(n)| + |ei(n)|
45*>
46*> W is a block diagonal matrix, with a 1 by 1 block for each real
47*> eigenvalue and a 2 by 2 block for each complex conjugate pair.
48*> If eigenvalues j and j+1 are a complex conjugate pair, so that
49*> WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 2 by 2
50*> block corresponding to the pair will be:
51*>
52*>    (  wr  wi  )
53*>    ( -wi  wr  )
54*>
55*> Such a block multiplying an n by 2 matrix ( ur ui ) on the right
56*> will be the same as multiplying  ur + i*ui  by  wr + i*wi.
57*>
58*> To handle various schemes for storage of left eigenvectors, there are
59*> options to use A-transpose instead of A, E-transpose instead of E,
60*> and/or W-transpose instead of W.
61*> \endverbatim
62*
63*  Arguments:
64*  ==========
65*
66*> \param[in] TRANSA
67*> \verbatim
68*>          TRANSA is CHARACTER*1
69*>          Specifies whether or not A is transposed.
70*>          = 'N':  No transpose
71*>          = 'T':  Transpose
72*>          = 'C':  Conjugate transpose (= Transpose)
73*> \endverbatim
74*>
75*> \param[in] TRANSE
76*> \verbatim
77*>          TRANSE is CHARACTER*1
78*>          Specifies whether or not E is transposed.
79*>          = 'N':  No transpose, eigenvectors are in columns of E
80*>          = 'T':  Transpose, eigenvectors are in rows of E
81*>          = 'C':  Conjugate transpose (= Transpose)
82*> \endverbatim
83*>
84*> \param[in] TRANSW
85*> \verbatim
86*>          TRANSW is CHARACTER*1
87*>          Specifies whether or not W is transposed.
88*>          = 'N':  No transpose
89*>          = 'T':  Transpose, use -WI(j) instead of WI(j)
90*>          = 'C':  Conjugate transpose, use -WI(j) instead of WI(j)
91*> \endverbatim
92*>
93*> \param[in] N
94*> \verbatim
95*>          N is INTEGER
96*>          The order of the matrix A.  N >= 0.
97*> \endverbatim
98*>
99*> \param[in] A
100*> \verbatim
101*>          A is REAL array, dimension (LDA,N)
102*>          The matrix whose eigenvectors are in E.
103*> \endverbatim
104*>
105*> \param[in] LDA
106*> \verbatim
107*>          LDA is INTEGER
108*>          The leading dimension of the array A.  LDA >= max(1,N).
109*> \endverbatim
110*>
111*> \param[in] E
112*> \verbatim
113*>          E is REAL array, dimension (LDE,N)
114*>          The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors
115*>          are stored in the columns of E, if TRANSE = 'T' or 'C', the
116*>          eigenvectors are stored in the rows of E.
117*> \endverbatim
118*>
119*> \param[in] LDE
120*> \verbatim
121*>          LDE is INTEGER
122*>          The leading dimension of the array E.  LDE >= max(1,N).
123*> \endverbatim
124*>
125*> \param[in] WR
126*> \verbatim
127*>          WR is REAL array, dimension (N)
128*> \endverbatim
129*>
130*> \param[in] WI
131*> \verbatim
132*>          WI is REAL array, dimension (N)
133*>
134*>          The real and imaginary parts of the eigenvalues of A.
135*>          Purely real eigenvalues are indicated by WI(j) = 0.
136*>          Complex conjugate pairs are indicated by WR(j)=WR(j+1) and
137*>          WI(j) = - WI(j+1) non-zero; the real part is assumed to be
138*>          stored in the j-th row/column and the imaginary part in
139*>          the (j+1)-th row/column.
140*> \endverbatim
141*>
142*> \param[out] WORK
143*> \verbatim
144*>          WORK is REAL array, dimension (N*(N+1))
145*> \endverbatim
146*>
147*> \param[out] RESULT
148*> \verbatim
149*>          RESULT is REAL array, dimension (2)
150*>          RESULT(1) = | A E  -  E W | / ( |A| |E| ulp )
151*>          RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
152*> \endverbatim
153*
154*  Authors:
155*  ========
156*
157*> \author Univ. of Tennessee
158*> \author Univ. of California Berkeley
159*> \author Univ. of Colorado Denver
160*> \author NAG Ltd.
161*
162*> \date November 2011
163*
164*> \ingroup single_eig
165*
166*  =====================================================================
167      SUBROUTINE SGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR,
168     $                   WI, WORK, RESULT )
169*
170*  -- LAPACK test routine (version 3.4.0) --
171*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
172*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
173*     November 2011
174*
175*     .. Scalar Arguments ..
176      CHARACTER          TRANSA, TRANSE, TRANSW
177      INTEGER            LDA, LDE, N
178*     ..
179*     .. Array Arguments ..
180      REAL               A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
181     $                   WORK( * ), WR( * )
182*     ..
183*
184*  =====================================================================
185*
186*     .. Parameters ..
187      REAL               ZERO, ONE
188      PARAMETER          ( ZERO = 0.0, ONE = 1.0 )
189*     ..
190*     .. Local Scalars ..
191      CHARACTER          NORMA, NORME
192      INTEGER            IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL,
193     $                   JVEC
194      REAL               ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
195     $                   ULP, UNFL
196*     ..
197*     .. Local Arrays ..
198      REAL               WMAT( 2, 2 )
199*     ..
200*     .. External Functions ..
201      LOGICAL            LSAME
202      REAL               SLAMCH, SLANGE
203      EXTERNAL           LSAME, SLAMCH, SLANGE
204*     ..
205*     .. External Subroutines ..
206      EXTERNAL           SAXPY, SGEMM, SLASET
207*     ..
208*     .. Intrinsic Functions ..
209      INTRINSIC          ABS, MAX, MIN, REAL
210*     ..
211*     .. Executable Statements ..
212*
213*     Initialize RESULT (in case N=0)
214*
215      RESULT( 1 ) = ZERO
216      RESULT( 2 ) = ZERO
217      IF( N.LE.0 )
218     $   RETURN
219*
220      UNFL = SLAMCH( 'Safe minimum' )
221      ULP = SLAMCH( 'Precision' )
222*
223      ITRNSE = 0
224      INCE = 1
225      NORMA = 'O'
226      NORME = 'O'
227*
228      IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN
229         NORMA = 'I'
230      END IF
231      IF( LSAME( TRANSE, 'T' ) .OR. LSAME( TRANSE, 'C' ) ) THEN
232         NORME = 'I'
233         ITRNSE = 1
234         INCE = LDE
235      END IF
236*
237*     Check normalization of E
238*
239      ENRMIN = ONE / ULP
240      ENRMAX = ZERO
241      IF( ITRNSE.EQ.0 ) THEN
242*
243*        Eigenvectors are column vectors.
244*
245         IPAIR = 0
246         DO 30 JVEC = 1, N
247            TEMP1 = ZERO
248            IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO )
249     $         IPAIR = 1
250            IF( IPAIR.EQ.1 ) THEN
251*
252*              Complex eigenvector
253*
254               DO 10 J = 1, N
255                  TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+
256     $                    ABS( E( J, JVEC+1 ) ) )
257   10          CONTINUE
258               ENRMIN = MIN( ENRMIN, TEMP1 )
259               ENRMAX = MAX( ENRMAX, TEMP1 )
260               IPAIR = 2
261            ELSE IF( IPAIR.EQ.2 ) THEN
262               IPAIR = 0
263            ELSE
264*
265*              Real eigenvector
266*
267               DO 20 J = 1, N
268                  TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) )
269   20          CONTINUE
270               ENRMIN = MIN( ENRMIN, TEMP1 )
271               ENRMAX = MAX( ENRMAX, TEMP1 )
272               IPAIR = 0
273            END IF
274   30    CONTINUE
275*
276      ELSE
277*
278*        Eigenvectors are row vectors.
279*
280         DO 40 JVEC = 1, N
281            WORK( JVEC ) = ZERO
282   40    CONTINUE
283*
284         DO 60 J = 1, N
285            IPAIR = 0
286            DO 50 JVEC = 1, N
287               IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO )
288     $            IPAIR = 1
289               IF( IPAIR.EQ.1 ) THEN
290                  WORK( JVEC ) = MAX( WORK( JVEC ),
291     $                           ABS( E( J, JVEC ) )+ABS( E( J,
292     $                           JVEC+1 ) ) )
293                  WORK( JVEC+1 ) = WORK( JVEC )
294               ELSE IF( IPAIR.EQ.2 ) THEN
295                  IPAIR = 0
296               ELSE
297                  WORK( JVEC ) = MAX( WORK( JVEC ),
298     $                           ABS( E( J, JVEC ) ) )
299                  IPAIR = 0
300               END IF
301   50       CONTINUE
302   60    CONTINUE
303*
304         DO 70 JVEC = 1, N
305            ENRMIN = MIN( ENRMIN, WORK( JVEC ) )
306            ENRMAX = MAX( ENRMAX, WORK( JVEC ) )
307   70    CONTINUE
308      END IF
309*
310*     Norm of A:
311*
312      ANORM = MAX( SLANGE( NORMA, N, N, A, LDA, WORK ), UNFL )
313*
314*     Norm of E:
315*
316      ENORM = MAX( SLANGE( NORME, N, N, E, LDE, WORK ), ULP )
317*
318*     Norm of error:
319*
320*     Error =  AE - EW
321*
322      CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
323*
324      IPAIR = 0
325      IEROW = 1
326      IECOL = 1
327*
328      DO 80 JCOL = 1, N
329         IF( ITRNSE.EQ.1 ) THEN
330            IEROW = JCOL
331         ELSE
332            IECOL = JCOL
333         END IF
334*
335         IF( IPAIR.EQ.0 .AND. WI( JCOL ).NE.ZERO )
336     $      IPAIR = 1
337*
338         IF( IPAIR.EQ.1 ) THEN
339            WMAT( 1, 1 ) = WR( JCOL )
340            WMAT( 2, 1 ) = -WI( JCOL )
341            WMAT( 1, 2 ) = WI( JCOL )
342            WMAT( 2, 2 ) = WR( JCOL )
343            CALL SGEMM( TRANSE, TRANSW, N, 2, 2, ONE, E( IEROW, IECOL ),
344     $                  LDE, WMAT, 2, ZERO, WORK( N*( JCOL-1 )+1 ), N )
345            IPAIR = 2
346         ELSE IF( IPAIR.EQ.2 ) THEN
347            IPAIR = 0
348*
349         ELSE
350*
351            CALL SAXPY( N, WR( JCOL ), E( IEROW, IECOL ), INCE,
352     $                  WORK( N*( JCOL-1 )+1 ), 1 )
353            IPAIR = 0
354         END IF
355*
356   80 CONTINUE
357*
358      CALL SGEMM( TRANSA, TRANSE, N, N, N, ONE, A, LDA, E, LDE, -ONE,
359     $            WORK, N )
360*
361      ERRNRM = SLANGE( 'One', N, N, WORK, N, WORK( N*N+1 ) ) / ENORM
362*
363*     Compute RESULT(1) (avoiding under/overflow)
364*
365      IF( ANORM.GT.ERRNRM ) THEN
366         RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP
367      ELSE
368         IF( ANORM.LT.ONE ) THEN
369            RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP
370         ELSE
371            RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP
372         END IF
373      END IF
374*
375*     Compute RESULT(2) : the normalization error in E.
376*
377      RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) /
378     $              ( REAL( N )*ULP )
379*
380      RETURN
381*
382*     End of SGET22
383*
384      END
385