LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
sget23.f
Go to the documentation of this file.
1 *> \brief \b SGET23
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 SGET23( COMP, BALANC, JTYPE, THRESH, ISEED, NOUNIT, N,
12 * A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, VR,
13 * LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
14 * RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
15 * WORK, LWORK, IWORK, INFO )
16 *
17 * .. Scalar Arguments ..
18 * LOGICAL COMP
19 * CHARACTER BALANC
20 * INTEGER INFO, JTYPE, LDA, LDLRE, LDVL, LDVR, LWORK, N,
21 * $ NOUNIT
22 * REAL THRESH
23 * ..
24 * .. Array Arguments ..
25 * INTEGER ISEED( 4 ), IWORK( * )
26 * REAL A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
27 * $ RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
28 * $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
29 * $ RESULT( 11 ), SCALE( * ), SCALE1( * ),
30 * $ VL( LDVL, * ), VR( LDVR, * ), WI( * ),
31 * $ WI1( * ), WORK( * ), WR( * ), WR1( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> SGET23 checks the nonsymmetric eigenvalue problem driver SGEEVX.
41 *> If COMP = .FALSE., the first 8 of the following tests will be
42 *> performed on the input matrix A, and also test 9 if LWORK is
43 *> sufficiently large.
44 *> if COMP is .TRUE. all 11 tests will be performed.
45 *>
46 *> (1) | A * VR - VR * W | / ( n |A| ulp )
47 *>
48 *> Here VR is the matrix of unit right eigenvectors.
49 *> W is a block diagonal matrix, with a 1x1 block for each
50 *> real eigenvalue and a 2x2 block for each complex conjugate
51 *> pair. If eigenvalues j and j+1 are a complex conjugate pair,
52 *> so WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the
53 *> 2 x 2 block corresponding to the pair will be:
54 *>
55 *> ( wr wi )
56 *> ( -wi wr )
57 *>
58 *> Such a block multiplying an n x 2 matrix ( ur ui ) on the
59 *> right will be the same as multiplying ur + i*ui by wr + i*wi.
60 *>
61 *> (2) | A**H * VL - VL * W**H | / ( n |A| ulp )
62 *>
63 *> Here VL is the matrix of unit left eigenvectors, A**H is the
64 *> conjugate transpose of A, and W is as above.
65 *>
66 *> (3) | |VR(i)| - 1 | / ulp and largest component real
67 *>
68 *> VR(i) denotes the i-th column of VR.
69 *>
70 *> (4) | |VL(i)| - 1 | / ulp and largest component real
71 *>
72 *> VL(i) denotes the i-th column of VL.
73 *>
74 *> (5) 0 if W(full) = W(partial), 1/ulp otherwise
75 *>
76 *> W(full) denotes the eigenvalues computed when VR, VL, RCONDV
77 *> and RCONDE are also computed, and W(partial) denotes the
78 *> eigenvalues computed when only some of VR, VL, RCONDV, and
79 *> RCONDE are computed.
80 *>
81 *> (6) 0 if VR(full) = VR(partial), 1/ulp otherwise
82 *>
83 *> VR(full) denotes the right eigenvectors computed when VL, RCONDV
84 *> and RCONDE are computed, and VR(partial) denotes the result
85 *> when only some of VL and RCONDV are computed.
86 *>
87 *> (7) 0 if VL(full) = VL(partial), 1/ulp otherwise
88 *>
89 *> VL(full) denotes the left eigenvectors computed when VR, RCONDV
90 *> and RCONDE are computed, and VL(partial) denotes the result
91 *> when only some of VR and RCONDV are computed.
92 *>
93 *> (8) 0 if SCALE, ILO, IHI, ABNRM (full) =
94 *> SCALE, ILO, IHI, ABNRM (partial)
95 *> 1/ulp otherwise
96 *>
97 *> SCALE, ILO, IHI and ABNRM describe how the matrix is balanced.
98 *> (full) is when VR, VL, RCONDE and RCONDV are also computed, and
99 *> (partial) is when some are not computed.
100 *>
101 *> (9) 0 if RCONDV(full) = RCONDV(partial), 1/ulp otherwise
102 *>
103 *> RCONDV(full) denotes the reciprocal condition numbers of the
104 *> right eigenvectors computed when VR, VL and RCONDE are also
105 *> computed. RCONDV(partial) denotes the reciprocal condition
106 *> numbers when only some of VR, VL and RCONDE are computed.
107 *>
108 *> (10) |RCONDV - RCDVIN| / cond(RCONDV)
109 *>
110 *> RCONDV is the reciprocal right eigenvector condition number
111 *> computed by SGEEVX and RCDVIN (the precomputed true value)
112 *> is supplied as input. cond(RCONDV) is the condition number of
113 *> RCONDV, and takes errors in computing RCONDV into account, so
114 *> that the resulting quantity should be O(ULP). cond(RCONDV) is
115 *> essentially given by norm(A)/RCONDE.
116 *>
117 *> (11) |RCONDE - RCDEIN| / cond(RCONDE)
118 *>
119 *> RCONDE is the reciprocal eigenvalue condition number
120 *> computed by SGEEVX and RCDEIN (the precomputed true value)
121 *> is supplied as input. cond(RCONDE) is the condition number
122 *> of RCONDE, and takes errors in computing RCONDE into account,
123 *> so that the resulting quantity should be O(ULP). cond(RCONDE)
124 *> is essentially given by norm(A)/RCONDV.
125 *> \endverbatim
126 *
127 * Arguments:
128 * ==========
129 *
130 *> \param[in] COMP
131 *> \verbatim
132 *> COMP is LOGICAL
133 *> COMP describes which input tests to perform:
134 *> = .FALSE. if the computed condition numbers are not to
135 *> be tested against RCDVIN and RCDEIN
136 *> = .TRUE. if they are to be compared
137 *> \endverbatim
138 *>
139 *> \param[in] BALANC
140 *> \verbatim
141 *> BALANC is CHARACTER
142 *> Describes the balancing option to be tested.
143 *> = 'N' for no permuting or diagonal scaling
144 *> = 'P' for permuting but no diagonal scaling
145 *> = 'S' for no permuting but diagonal scaling
146 *> = 'B' for permuting and diagonal scaling
147 *> \endverbatim
148 *>
149 *> \param[in] JTYPE
150 *> \verbatim
151 *> JTYPE is INTEGER
152 *> Type of input matrix. Used to label output if error occurs.
153 *> \endverbatim
154 *>
155 *> \param[in] THRESH
156 *> \verbatim
157 *> THRESH is REAL
158 *> A test will count as "failed" if the "error", computed as
159 *> described above, exceeds THRESH. Note that the error
160 *> is scaled to be O(1), so THRESH should be a reasonably
161 *> small multiple of 1, e.g., 10 or 100. In particular,
162 *> it should not depend on the precision (single vs. double)
163 *> or the size of the matrix. It must be at least zero.
164 *> \endverbatim
165 *>
166 *> \param[in] ISEED
167 *> \verbatim
168 *> ISEED is INTEGER array, dimension (4)
169 *> If COMP = .FALSE., the random number generator seed
170 *> used to produce matrix.
171 *> If COMP = .TRUE., ISEED(1) = the number of the example.
172 *> Used to label output if error occurs.
173 *> \endverbatim
174 *>
175 *> \param[in] NOUNIT
176 *> \verbatim
177 *> NOUNIT is INTEGER
178 *> The FORTRAN unit number for printing out error messages
179 *> (e.g., if a routine returns INFO not equal to 0.)
180 *> \endverbatim
181 *>
182 *> \param[in] N
183 *> \verbatim
184 *> N is INTEGER
185 *> The dimension of A. N must be at least 0.
186 *> \endverbatim
187 *>
188 *> \param[in,out] A
189 *> \verbatim
190 *> A is REAL array, dimension (LDA,N)
191 *> Used to hold the matrix whose eigenvalues are to be
192 *> computed.
193 *> \endverbatim
194 *>
195 *> \param[in] LDA
196 *> \verbatim
197 *> LDA is INTEGER
198 *> The leading dimension of A, and H. LDA must be at
199 *> least 1 and at least N.
200 *> \endverbatim
201 *>
202 *> \param[out] H
203 *> \verbatim
204 *> H is REAL array, dimension (LDA,N)
205 *> Another copy of the test matrix A, modified by SGEEVX.
206 *> \endverbatim
207 *>
208 *> \param[out] WR
209 *> \verbatim
210 *> WR is REAL array, dimension (N)
211 *> \endverbatim
212 *>
213 *> \param[out] WI
214 *> \verbatim
215 *> WI is REAL array, dimension (N)
216 *>
217 *> The real and imaginary parts of the eigenvalues of A.
218 *> On exit, WR + WI*i are the eigenvalues of the matrix in A.
219 *> \endverbatim
220 *>
221 *> \param[out] WR1
222 *> \verbatim
223 *> WR1 is REAL array, dimension (N)
224 *> \endverbatim
225 *>
226 *> \param[out] WI1
227 *> \verbatim
228 *> WI1 is REAL array, dimension (N)
229 *>
230 *> Like WR, WI, these arrays contain the eigenvalues of A,
231 *> but those computed when SGEEVX only computes a partial
232 *> eigendecomposition, i.e. not the eigenvalues and left
233 *> and right eigenvectors.
234 *> \endverbatim
235 *>
236 *> \param[out] VL
237 *> \verbatim
238 *> VL is REAL array, dimension (LDVL,N)
239 *> VL holds the computed left eigenvectors.
240 *> \endverbatim
241 *>
242 *> \param[in] LDVL
243 *> \verbatim
244 *> LDVL is INTEGER
245 *> Leading dimension of VL. Must be at least max(1,N).
246 *> \endverbatim
247 *>
248 *> \param[out] VR
249 *> \verbatim
250 *> VR is REAL array, dimension (LDVR,N)
251 *> VR holds the computed right eigenvectors.
252 *> \endverbatim
253 *>
254 *> \param[in] LDVR
255 *> \verbatim
256 *> LDVR is INTEGER
257 *> Leading dimension of VR. Must be at least max(1,N).
258 *> \endverbatim
259 *>
260 *> \param[out] LRE
261 *> \verbatim
262 *> LRE is REAL array, dimension (LDLRE,N)
263 *> LRE holds the computed right or left eigenvectors.
264 *> \endverbatim
265 *>
266 *> \param[in] LDLRE
267 *> \verbatim
268 *> LDLRE is INTEGER
269 *> Leading dimension of LRE. Must be at least max(1,N).
270 *> \endverbatim
271 *>
272 *> \param[out] RCONDV
273 *> \verbatim
274 *> RCONDV is REAL array, dimension (N)
275 *> RCONDV holds the computed reciprocal condition numbers
276 *> for eigenvectors.
277 *> \endverbatim
278 *>
279 *> \param[out] RCNDV1
280 *> \verbatim
281 *> RCNDV1 is REAL array, dimension (N)
282 *> RCNDV1 holds more computed reciprocal condition numbers
283 *> for eigenvectors.
284 *> \endverbatim
285 *>
286 *> \param[in] RCDVIN
287 *> \verbatim
288 *> RCDVIN is REAL array, dimension (N)
289 *> When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
290 *> condition numbers for eigenvectors to be compared with
291 *> RCONDV.
292 *> \endverbatim
293 *>
294 *> \param[out] RCONDE
295 *> \verbatim
296 *> RCONDE is REAL array, dimension (N)
297 *> RCONDE holds the computed reciprocal condition numbers
298 *> for eigenvalues.
299 *> \endverbatim
300 *>
301 *> \param[out] RCNDE1
302 *> \verbatim
303 *> RCNDE1 is REAL array, dimension (N)
304 *> RCNDE1 holds more computed reciprocal condition numbers
305 *> for eigenvalues.
306 *> \endverbatim
307 *>
308 *> \param[in] RCDEIN
309 *> \verbatim
310 *> RCDEIN is REAL array, dimension (N)
311 *> When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
312 *> condition numbers for eigenvalues to be compared with
313 *> RCONDE.
314 *> \endverbatim
315 *>
316 *> \param[out] SCALE
317 *> \verbatim
318 *> SCALE is REAL array, dimension (N)
319 *> Holds information describing balancing of matrix.
320 *> \endverbatim
321 *>
322 *> \param[out] SCALE1
323 *> \verbatim
324 *> SCALE1 is REAL array, dimension (N)
325 *> Holds information describing balancing of matrix.
326 *> \endverbatim
327 *>
328 *> \param[out] RESULT
329 *> \verbatim
330 *> RESULT is REAL array, dimension (11)
331 *> The values computed by the 11 tests described above.
332 *> The values are currently limited to 1/ulp, to avoid
333 *> overflow.
334 *> \endverbatim
335 *>
336 *> \param[out] WORK
337 *> \verbatim
338 *> WORK is REAL array, dimension (LWORK)
339 *> \endverbatim
340 *>
341 *> \param[in] LWORK
342 *> \verbatim
343 *> LWORK is INTEGER
344 *> The number of entries in WORK. This must be at least
345 *> 3*N, and 6*N+N**2 if tests 9, 10 or 11 are to be performed.
346 *> \endverbatim
347 *>
348 *> \param[out] IWORK
349 *> \verbatim
350 *> IWORK is INTEGER array, dimension (2*N)
351 *> \endverbatim
352 *>
353 *> \param[out] INFO
354 *> \verbatim
355 *> INFO is INTEGER
356 *> If 0, successful exit.
357 *> If <0, input parameter -INFO had an incorrect value.
358 *> If >0, SGEEVX returned an error code, the absolute
359 *> value of which is returned.
360 *> \endverbatim
361 *
362 * Authors:
363 * ========
364 *
365 *> \author Univ. of Tennessee
366 *> \author Univ. of California Berkeley
367 *> \author Univ. of Colorado Denver
368 *> \author NAG Ltd.
369 *
370 *> \ingroup single_eig
371 *
372 * =====================================================================
373  SUBROUTINE sget23( COMP, BALANC, JTYPE, THRESH, ISEED, NOUNIT, N,
374  $ A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, VR,
375  $ LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
376  $ RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
377  $ WORK, LWORK, IWORK, INFO )
378 *
379 * -- LAPACK test routine --
380 * -- LAPACK is a software package provided by Univ. of Tennessee, --
381 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
382 *
383 * .. Scalar Arguments ..
384  LOGICAL COMP
385  CHARACTER BALANC
386  INTEGER INFO, JTYPE, LDA, LDLRE, LDVL, LDVR, LWORK, N,
387  $ NOUNIT
388  REAL THRESH
389 * ..
390 * .. Array Arguments ..
391  INTEGER ISEED( 4 ), IWORK( * )
392  REAL A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
393  $ RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
394  $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
395  $ result( 11 ), scale( * ), scale1( * ),
396  $ vl( ldvl, * ), vr( ldvr, * ), wi( * ),
397  $ wi1( * ), work( * ), wr( * ), wr1( * )
398 * ..
399 *
400 * =====================================================================
401 *
402 *
403 * .. Parameters ..
404  REAL ZERO, ONE, TWO
405  PARAMETER ( ZERO = 0.0e0, one = 1.0e0, two = 2.0e0 )
406  REAL EPSIN
407  PARAMETER ( EPSIN = 5.9605e-8 )
408 * ..
409 * .. Local Scalars ..
410  LOGICAL BALOK, NOBAL
411  CHARACTER SENSE
412  INTEGER I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM,
413  $ J, JJ, KMIN
414  REAL ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN,
415  $ ulp, ulpinv, v, vimin, vmax, vmx, vrmin, vrmx,
416  $ vtst
417 * ..
418 * .. Local Arrays ..
419  CHARACTER SENS( 2 )
420  REAL DUM( 1 ), RES( 2 )
421 * ..
422 * .. External Functions ..
423  LOGICAL LSAME
424  REAL SLAMCH, SLAPY2, SNRM2
425  EXTERNAL LSAME, SLAMCH, SLAPY2, SNRM2
426 * ..
427 * .. External Subroutines ..
428  EXTERNAL sgeevx, sget22, slacpy, xerbla
429 * ..
430 * .. Intrinsic Functions ..
431  INTRINSIC abs, max, min, real
432 * ..
433 * .. Data statements ..
434  DATA sens / 'N', 'V' /
435 * ..
436 * .. Executable Statements ..
437 *
438 * Check for errors
439 *
440  nobal = lsame( balanc, 'N' )
441  balok = nobal .OR. lsame( balanc, 'P' ) .OR.
442  $ lsame( balanc, 'S' ) .OR. lsame( balanc, 'B' )
443  info = 0
444  IF( .NOT.balok ) THEN
445  info = -2
446  ELSE IF( thresh.LT.zero ) THEN
447  info = -4
448  ELSE IF( nounit.LE.0 ) THEN
449  info = -6
450  ELSE IF( n.LT.0 ) THEN
451  info = -7
452  ELSE IF( lda.LT.1 .OR. lda.LT.n ) THEN
453  info = -9
454  ELSE IF( ldvl.LT.1 .OR. ldvl.LT.n ) THEN
455  info = -16
456  ELSE IF( ldvr.LT.1 .OR. ldvr.LT.n ) THEN
457  info = -18
458  ELSE IF( ldlre.LT.1 .OR. ldlre.LT.n ) THEN
459  info = -20
460  ELSE IF( lwork.LT.3*n .OR. ( comp .AND. lwork.LT.6*n+n*n ) ) THEN
461  info = -31
462  END IF
463 *
464  IF( info.NE.0 ) THEN
465  CALL xerbla( 'SGET23', -info )
466  RETURN
467  END IF
468 *
469 * Quick return if nothing to do
470 *
471  DO 10 i = 1, 11
472  result( i ) = -one
473  10 CONTINUE
474 *
475  IF( n.EQ.0 )
476  $ RETURN
477 *
478 * More Important constants
479 *
480  ulp = slamch( 'Precision' )
481  smlnum = slamch( 'S' )
482  ulpinv = one / ulp
483 *
484 * Compute eigenvalues and eigenvectors, and test them
485 *
486  IF( lwork.GE.6*n+n*n ) THEN
487  sense = 'B'
488  isensm = 2
489  ELSE
490  sense = 'E'
491  isensm = 1
492  END IF
493  CALL slacpy( 'F', n, n, a, lda, h, lda )
494  CALL sgeevx( balanc, 'V', 'V', sense, n, h, lda, wr, wi, vl, ldvl,
495  $ vr, ldvr, ilo, ihi, scale, abnrm, rconde, rcondv,
496  $ work, lwork, iwork, iinfo )
497  IF( iinfo.NE.0 ) THEN
498  result( 1 ) = ulpinv
499  IF( jtype.NE.22 ) THEN
500  WRITE( nounit, fmt = 9998 )'SGEEVX1', iinfo, n, jtype,
501  $ balanc, iseed
502  ELSE
503  WRITE( nounit, fmt = 9999 )'SGEEVX1', iinfo, n, iseed( 1 )
504  END IF
505  info = abs( iinfo )
506  RETURN
507  END IF
508 *
509 * Do Test (1)
510 *
511  CALL sget22( 'N', 'N', 'N', n, a, lda, vr, ldvr, wr, wi, work,
512  $ res )
513  result( 1 ) = res( 1 )
514 *
515 * Do Test (2)
516 *
517  CALL sget22( 'T', 'N', 'T', n, a, lda, vl, ldvl, wr, wi, work,
518  $ res )
519  result( 2 ) = res( 1 )
520 *
521 * Do Test (3)
522 *
523  DO 30 j = 1, n
524  tnrm = one
525  IF( wi( j ).EQ.zero ) THEN
526  tnrm = snrm2( n, vr( 1, j ), 1 )
527  ELSE IF( wi( j ).GT.zero ) THEN
528  tnrm = slapy2( snrm2( n, vr( 1, j ), 1 ),
529  $ snrm2( n, vr( 1, j+1 ), 1 ) )
530  END IF
531  result( 3 ) = max( result( 3 ),
532  $ min( ulpinv, abs( tnrm-one ) / ulp ) )
533  IF( wi( j ).GT.zero ) THEN
534  vmx = zero
535  vrmx = zero
536  DO 20 jj = 1, n
537  vtst = slapy2( vr( jj, j ), vr( jj, j+1 ) )
538  IF( vtst.GT.vmx )
539  $ vmx = vtst
540  IF( vr( jj, j+1 ).EQ.zero .AND. abs( vr( jj, j ) ).GT.
541  $ vrmx )vrmx = abs( vr( jj, j ) )
542  20 CONTINUE
543  IF( vrmx / vmx.LT.one-two*ulp )
544  $ result( 3 ) = ulpinv
545  END IF
546  30 CONTINUE
547 *
548 * Do Test (4)
549 *
550  DO 50 j = 1, n
551  tnrm = one
552  IF( wi( j ).EQ.zero ) THEN
553  tnrm = snrm2( n, vl( 1, j ), 1 )
554  ELSE IF( wi( j ).GT.zero ) THEN
555  tnrm = slapy2( snrm2( n, vl( 1, j ), 1 ),
556  $ snrm2( n, vl( 1, j+1 ), 1 ) )
557  END IF
558  result( 4 ) = max( result( 4 ),
559  $ min( ulpinv, abs( tnrm-one ) / ulp ) )
560  IF( wi( j ).GT.zero ) THEN
561  vmx = zero
562  vrmx = zero
563  DO 40 jj = 1, n
564  vtst = slapy2( vl( jj, j ), vl( jj, j+1 ) )
565  IF( vtst.GT.vmx )
566  $ vmx = vtst
567  IF( vl( jj, j+1 ).EQ.zero .AND. abs( vl( jj, j ) ).GT.
568  $ vrmx )vrmx = abs( vl( jj, j ) )
569  40 CONTINUE
570  IF( vrmx / vmx.LT.one-two*ulp )
571  $ result( 4 ) = ulpinv
572  END IF
573  50 CONTINUE
574 *
575 * Test for all options of computing condition numbers
576 *
577  DO 200 isens = 1, isensm
578 *
579  sense = sens( isens )
580 *
581 * Compute eigenvalues only, and test them
582 *
583  CALL slacpy( 'F', n, n, a, lda, h, lda )
584  CALL sgeevx( balanc, 'N', 'N', sense, n, h, lda, wr1, wi1, dum,
585  $ 1, dum, 1, ilo1, ihi1, scale1, abnrm1, rcnde1,
586  $ rcndv1, work, lwork, iwork, iinfo )
587  IF( iinfo.NE.0 ) THEN
588  result( 1 ) = ulpinv
589  IF( jtype.NE.22 ) THEN
590  WRITE( nounit, fmt = 9998 )'SGEEVX2', iinfo, n, jtype,
591  $ balanc, iseed
592  ELSE
593  WRITE( nounit, fmt = 9999 )'SGEEVX2', iinfo, n,
594  $ iseed( 1 )
595  END IF
596  info = abs( iinfo )
597  GO TO 190
598  END IF
599 *
600 * Do Test (5)
601 *
602  DO 60 j = 1, n
603  IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
604  $ result( 5 ) = ulpinv
605  60 CONTINUE
606 *
607 * Do Test (8)
608 *
609  IF( .NOT.nobal ) THEN
610  DO 70 j = 1, n
611  IF( scale( j ).NE.scale1( j ) )
612  $ result( 8 ) = ulpinv
613  70 CONTINUE
614  IF( ilo.NE.ilo1 )
615  $ result( 8 ) = ulpinv
616  IF( ihi.NE.ihi1 )
617  $ result( 8 ) = ulpinv
618  IF( abnrm.NE.abnrm1 )
619  $ result( 8 ) = ulpinv
620  END IF
621 *
622 * Do Test (9)
623 *
624  IF( isens.EQ.2 .AND. n.GT.1 ) THEN
625  DO 80 j = 1, n
626  IF( rcondv( j ).NE.rcndv1( j ) )
627  $ result( 9 ) = ulpinv
628  80 CONTINUE
629  END IF
630 *
631 * Compute eigenvalues and right eigenvectors, and test them
632 *
633  CALL slacpy( 'F', n, n, a, lda, h, lda )
634  CALL sgeevx( balanc, 'N', 'V', sense, n, h, lda, wr1, wi1, dum,
635  $ 1, lre, ldlre, ilo1, ihi1, scale1, abnrm1, rcnde1,
636  $ rcndv1, work, lwork, iwork, iinfo )
637  IF( iinfo.NE.0 ) THEN
638  result( 1 ) = ulpinv
639  IF( jtype.NE.22 ) THEN
640  WRITE( nounit, fmt = 9998 )'SGEEVX3', iinfo, n, jtype,
641  $ balanc, iseed
642  ELSE
643  WRITE( nounit, fmt = 9999 )'SGEEVX3', iinfo, n,
644  $ iseed( 1 )
645  END IF
646  info = abs( iinfo )
647  GO TO 190
648  END IF
649 *
650 * Do Test (5) again
651 *
652  DO 90 j = 1, n
653  IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
654  $ result( 5 ) = ulpinv
655  90 CONTINUE
656 *
657 * Do Test (6)
658 *
659  DO 110 j = 1, n
660  DO 100 jj = 1, n
661  IF( vr( j, jj ).NE.lre( j, jj ) )
662  $ result( 6 ) = ulpinv
663  100 CONTINUE
664  110 CONTINUE
665 *
666 * Do Test (8) again
667 *
668  IF( .NOT.nobal ) THEN
669  DO 120 j = 1, n
670  IF( scale( j ).NE.scale1( j ) )
671  $ result( 8 ) = ulpinv
672  120 CONTINUE
673  IF( ilo.NE.ilo1 )
674  $ result( 8 ) = ulpinv
675  IF( ihi.NE.ihi1 )
676  $ result( 8 ) = ulpinv
677  IF( abnrm.NE.abnrm1 )
678  $ result( 8 ) = ulpinv
679  END IF
680 *
681 * Do Test (9) again
682 *
683  IF( isens.EQ.2 .AND. n.GT.1 ) THEN
684  DO 130 j = 1, n
685  IF( rcondv( j ).NE.rcndv1( j ) )
686  $ result( 9 ) = ulpinv
687  130 CONTINUE
688  END IF
689 *
690 * Compute eigenvalues and left eigenvectors, and test them
691 *
692  CALL slacpy( 'F', n, n, a, lda, h, lda )
693  CALL sgeevx( balanc, 'V', 'N', sense, n, h, lda, wr1, wi1, lre,
694  $ ldlre, dum, 1, ilo1, ihi1, scale1, abnrm1, rcnde1,
695  $ rcndv1, work, lwork, iwork, iinfo )
696  IF( iinfo.NE.0 ) THEN
697  result( 1 ) = ulpinv
698  IF( jtype.NE.22 ) THEN
699  WRITE( nounit, fmt = 9998 )'SGEEVX4', iinfo, n, jtype,
700  $ balanc, iseed
701  ELSE
702  WRITE( nounit, fmt = 9999 )'SGEEVX4', iinfo, n,
703  $ iseed( 1 )
704  END IF
705  info = abs( iinfo )
706  GO TO 190
707  END IF
708 *
709 * Do Test (5) again
710 *
711  DO 140 j = 1, n
712  IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
713  $ result( 5 ) = ulpinv
714  140 CONTINUE
715 *
716 * Do Test (7)
717 *
718  DO 160 j = 1, n
719  DO 150 jj = 1, n
720  IF( vl( j, jj ).NE.lre( j, jj ) )
721  $ result( 7 ) = ulpinv
722  150 CONTINUE
723  160 CONTINUE
724 *
725 * Do Test (8) again
726 *
727  IF( .NOT.nobal ) THEN
728  DO 170 j = 1, n
729  IF( scale( j ).NE.scale1( j ) )
730  $ result( 8 ) = ulpinv
731  170 CONTINUE
732  IF( ilo.NE.ilo1 )
733  $ result( 8 ) = ulpinv
734  IF( ihi.NE.ihi1 )
735  $ result( 8 ) = ulpinv
736  IF( abnrm.NE.abnrm1 )
737  $ result( 8 ) = ulpinv
738  END IF
739 *
740 * Do Test (9) again
741 *
742  IF( isens.EQ.2 .AND. n.GT.1 ) THEN
743  DO 180 j = 1, n
744  IF( rcondv( j ).NE.rcndv1( j ) )
745  $ result( 9 ) = ulpinv
746  180 CONTINUE
747  END IF
748 *
749  190 CONTINUE
750 *
751  200 CONTINUE
752 *
753 * If COMP, compare condition numbers to precomputed ones
754 *
755  IF( comp ) THEN
756  CALL slacpy( 'F', n, n, a, lda, h, lda )
757  CALL sgeevx( 'N', 'V', 'V', 'B', n, h, lda, wr, wi, vl, ldvl,
758  $ vr, ldvr, ilo, ihi, scale, abnrm, rconde, rcondv,
759  $ work, lwork, iwork, iinfo )
760  IF( iinfo.NE.0 ) THEN
761  result( 1 ) = ulpinv
762  WRITE( nounit, fmt = 9999 )'SGEEVX5', iinfo, n, iseed( 1 )
763  info = abs( iinfo )
764  GO TO 250
765  END IF
766 *
767 * Sort eigenvalues and condition numbers lexicographically
768 * to compare with inputs
769 *
770  DO 220 i = 1, n - 1
771  kmin = i
772  vrmin = wr( i )
773  vimin = wi( i )
774  DO 210 j = i + 1, n
775  IF( wr( j ).LT.vrmin ) THEN
776  kmin = j
777  vrmin = wr( j )
778  vimin = wi( j )
779  END IF
780  210 CONTINUE
781  wr( kmin ) = wr( i )
782  wi( kmin ) = wi( i )
783  wr( i ) = vrmin
784  wi( i ) = vimin
785  vrmin = rconde( kmin )
786  rconde( kmin ) = rconde( i )
787  rconde( i ) = vrmin
788  vrmin = rcondv( kmin )
789  rcondv( kmin ) = rcondv( i )
790  rcondv( i ) = vrmin
791  220 CONTINUE
792 *
793 * Compare condition numbers for eigenvectors
794 * taking their condition numbers into account
795 *
796  result( 10 ) = zero
797  eps = max( epsin, ulp )
798  v = max( real( n )*eps*abnrm, smlnum )
799  IF( abnrm.EQ.zero )
800  $ v = one
801  DO 230 i = 1, n
802  IF( v.GT.rcondv( i )*rconde( i ) ) THEN
803  tol = rcondv( i )
804  ELSE
805  tol = v / rconde( i )
806  END IF
807  IF( v.GT.rcdvin( i )*rcdein( i ) ) THEN
808  tolin = rcdvin( i )
809  ELSE
810  tolin = v / rcdein( i )
811  END IF
812  tol = max( tol, smlnum / eps )
813  tolin = max( tolin, smlnum / eps )
814  IF( eps*( rcdvin( i )-tolin ).GT.rcondv( i )+tol ) THEN
815  vmax = one / eps
816  ELSE IF( rcdvin( i )-tolin.GT.rcondv( i )+tol ) THEN
817  vmax = ( rcdvin( i )-tolin ) / ( rcondv( i )+tol )
818  ELSE IF( rcdvin( i )+tolin.LT.eps*( rcondv( i )-tol ) ) THEN
819  vmax = one / eps
820  ELSE IF( rcdvin( i )+tolin.LT.rcondv( i )-tol ) THEN
821  vmax = ( rcondv( i )-tol ) / ( rcdvin( i )+tolin )
822  ELSE
823  vmax = one
824  END IF
825  result( 10 ) = max( result( 10 ), vmax )
826  230 CONTINUE
827 *
828 * Compare condition numbers for eigenvalues
829 * taking their condition numbers into account
830 *
831  result( 11 ) = zero
832  DO 240 i = 1, n
833  IF( v.GT.rcondv( i ) ) THEN
834  tol = one
835  ELSE
836  tol = v / rcondv( i )
837  END IF
838  IF( v.GT.rcdvin( i ) ) THEN
839  tolin = one
840  ELSE
841  tolin = v / rcdvin( i )
842  END IF
843  tol = max( tol, smlnum / eps )
844  tolin = max( tolin, smlnum / eps )
845  IF( eps*( rcdein( i )-tolin ).GT.rconde( i )+tol ) THEN
846  vmax = one / eps
847  ELSE IF( rcdein( i )-tolin.GT.rconde( i )+tol ) THEN
848  vmax = ( rcdein( i )-tolin ) / ( rconde( i )+tol )
849  ELSE IF( rcdein( i )+tolin.LT.eps*( rconde( i )-tol ) ) THEN
850  vmax = one / eps
851  ELSE IF( rcdein( i )+tolin.LT.rconde( i )-tol ) THEN
852  vmax = ( rconde( i )-tol ) / ( rcdein( i )+tolin )
853  ELSE
854  vmax = one
855  END IF
856  result( 11 ) = max( result( 11 ), vmax )
857  240 CONTINUE
858  250 CONTINUE
859 *
860  END IF
861 *
862  9999 FORMAT( ' SGET23: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
863  $ i6, ', INPUT EXAMPLE NUMBER = ', i4 )
864  9998 FORMAT( ' SGET23: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
865  $ i6, ', JTYPE=', i6, ', BALANC = ', a, ', ISEED=(',
866  $ 3( i5, ',' ), i5, ')' )
867 *
868  RETURN
869 *
870 * End of SGET23
871 *
872  END
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sgeevx(BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV, WORK, LWORK, IWORK, INFO)
SGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition: sgeevx.f:306
subroutine sget22(TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, WI, WORK, RESULT)
SGET22
Definition: sget22.f:168
subroutine sget23(COMP, BALANC, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, VR, LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT, WORK, LWORK, IWORK, INFO)
SGET23
Definition: sget23.f:378