1 SUBROUTINE ZGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
2 $ H, HT, W, WT, WTMP, VS, LDVS, VS1, RCDEIN,
3 $ RCDVIN, NSLCT, ISLCT, ISRT, RESULT, WORK,
4 $ LWORK, RWORK, BWORK, INFO )
5 *
6 * -- LAPACK test routine (version 3.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 LOGICAL COMP
12 INTEGER INFO, ISRT, JTYPE, LDA, LDVS, LWORK, N, NOUNIT,
13 $ NSLCT
14 DOUBLE PRECISION RCDEIN, RCDVIN, THRESH
15 * ..
16 * .. Array Arguments ..
17 LOGICAL BWORK( * )
18 INTEGER ISEED( 4 ), ISLCT( * )
19 DOUBLE PRECISION RESULT( 17 ), RWORK( * )
20 COMPLEX*16 A( LDA, * ), H( LDA, * ), HT( LDA, * ),
21 $ VS( LDVS, * ), VS1( LDVS, * ), W( * ),
22 $ WORK( * ), WT( * ), WTMP( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * ZGET24 checks the nonsymmetric eigenvalue (Schur form) problem
29 * expert driver ZGEESX.
30 *
31 * If COMP = .FALSE., the first 13 of the following tests will be
32 * be performed on the input matrix A, and also tests 14 and 15
33 * if LWORK is sufficiently large.
34 * If COMP = .TRUE., all 17 test will be performed.
35 *
36 * (1) 0 if T is in Schur form, 1/ulp otherwise
37 * (no sorting of eigenvalues)
38 *
39 * (2) | A - VS T VS' | / ( n |A| ulp )
40 *
41 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
42 * form (no sorting of eigenvalues).
43 *
44 * (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
45 *
46 * (4) 0 if W are eigenvalues of T
47 * 1/ulp otherwise
48 * (no sorting of eigenvalues)
49 *
50 * (5) 0 if T(with VS) = T(without VS),
51 * 1/ulp otherwise
52 * (no sorting of eigenvalues)
53 *
54 * (6) 0 if eigenvalues(with VS) = eigenvalues(without VS),
55 * 1/ulp otherwise
56 * (no sorting of eigenvalues)
57 *
58 * (7) 0 if T is in Schur form, 1/ulp otherwise
59 * (with sorting of eigenvalues)
60 *
61 * (8) | A - VS T VS' | / ( n |A| ulp )
62 *
63 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
64 * form (with sorting of eigenvalues).
65 *
66 * (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
67 *
68 * (10) 0 if W are eigenvalues of T
69 * 1/ulp otherwise
70 * If workspace sufficient, also compare W with and
71 * without reciprocal condition numbers
72 * (with sorting of eigenvalues)
73 *
74 * (11) 0 if T(with VS) = T(without VS),
75 * 1/ulp otherwise
76 * If workspace sufficient, also compare T with and without
77 * reciprocal condition numbers
78 * (with sorting of eigenvalues)
79 *
80 * (12) 0 if eigenvalues(with VS) = eigenvalues(without VS),
81 * 1/ulp otherwise
82 * If workspace sufficient, also compare VS with and without
83 * reciprocal condition numbers
84 * (with sorting of eigenvalues)
85 *
86 * (13) if sorting worked and SDIM is the number of
87 * eigenvalues which were SELECTed
88 * If workspace sufficient, also compare SDIM with and
89 * without reciprocal condition numbers
90 *
91 * (14) if RCONDE the same no matter if VS and/or RCONDV computed
92 *
93 * (15) if RCONDV the same no matter if VS and/or RCONDE computed
94 *
95 * (16) |RCONDE - RCDEIN| / cond(RCONDE)
96 *
97 * RCONDE is the reciprocal average eigenvalue condition number
98 * computed by ZGEESX and RCDEIN (the precomputed true value)
99 * is supplied as input. cond(RCONDE) is the condition number
100 * of RCONDE, and takes errors in computing RCONDE into account,
101 * so that the resulting quantity should be O(ULP). cond(RCONDE)
102 * is essentially given by norm(A)/RCONDV.
103 *
104 * (17) |RCONDV - RCDVIN| / cond(RCONDV)
105 *
106 * RCONDV is the reciprocal right invariant subspace condition
107 * number computed by ZGEESX and RCDVIN (the precomputed true
108 * value) is supplied as input. cond(RCONDV) is the condition
109 * number of RCONDV, and takes errors in computing RCONDV into
110 * account, so that the resulting quantity should be O(ULP).
111 * cond(RCONDV) is essentially given by norm(A)/RCONDE.
112 *
113 * Arguments
114 * =========
115 *
116 * COMP (input) LOGICAL
117 * COMP describes which input tests to perform:
118 * = .FALSE. if the computed condition numbers are not to
119 * be tested against RCDVIN and RCDEIN
120 * = .TRUE. if they are to be compared
121 *
122 * JTYPE (input) INTEGER
123 * Type of input matrix. Used to label output if error occurs.
124 *
125 * ISEED (input) INTEGER array, dimension (4)
126 * If COMP = .FALSE., the random number generator seed
127 * used to produce matrix.
128 * If COMP = .TRUE., ISEED(1) = the number of the example.
129 * Used to label output if error occurs.
130 *
131 * THRESH (input) DOUBLE PRECISION
132 * A test will count as "failed" if the "error", computed as
133 * described above, exceeds THRESH. Note that the error
134 * is scaled to be O(1), so THRESH should be a reasonably
135 * small multiple of 1, e.g., 10 or 100. In particular,
136 * it should not depend on the precision (single vs. double)
137 * or the size of the matrix. It must be at least zero.
138 *
139 * NOUNIT (input) INTEGER
140 * The FORTRAN unit number for printing out error messages
141 * (e.g., if a routine returns INFO not equal to 0.)
142 *
143 * N (input) INTEGER
144 * The dimension of A. N must be at least 0.
145 *
146 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
147 * Used to hold the matrix whose eigenvalues are to be
148 * computed.
149 *
150 * LDA (input) INTEGER
151 * The leading dimension of A, and H. LDA must be at
152 * least 1 and at least N.
153 *
154 * H (workspace) COMPLEX*16 array, dimension (LDA, N)
155 * Another copy of the test matrix A, modified by ZGEESX.
156 *
157 * HT (workspace) COMPLEX*16 array, dimension (LDA, N)
158 * Yet another copy of the test matrix A, modified by ZGEESX.
159 *
160 * W (workspace) COMPLEX*16 array, dimension (N)
161 * The computed eigenvalues of A.
162 *
163 * WT (workspace) COMPLEX*16 array, dimension (N)
164 * Like W, this array contains the eigenvalues of A,
165 * but those computed when ZGEESX only computes a partial
166 * eigendecomposition, i.e. not Schur vectors
167 *
168 * WTMP (workspace) COMPLEX*16 array, dimension (N)
169 * Like W, this array contains the eigenvalues of A,
170 * but sorted by increasing real or imaginary part.
171 *
172 * VS (workspace) COMPLEX*16 array, dimension (LDVS, N)
173 * VS holds the computed Schur vectors.
174 *
175 * LDVS (input) INTEGER
176 * Leading dimension of VS. Must be at least max(1, N).
177 *
178 * VS1 (workspace) COMPLEX*16 array, dimension (LDVS, N)
179 * VS1 holds another copy of the computed Schur vectors.
180 *
181 * RCDEIN (input) DOUBLE PRECISION
182 * When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
183 * condition number for the average of selected eigenvalues.
184 *
185 * RCDVIN (input) DOUBLE PRECISION
186 * When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
187 * condition number for the selected right invariant subspace.
188 *
189 * NSLCT (input) INTEGER
190 * When COMP = .TRUE. the number of selected eigenvalues
191 * corresponding to the precomputed values RCDEIN and RCDVIN.
192 *
193 * ISLCT (input) INTEGER array, dimension (NSLCT)
194 * When COMP = .TRUE. ISLCT selects the eigenvalues of the
195 * input matrix corresponding to the precomputed values RCDEIN
196 * and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
197 * eigenvalue with the J-th largest real or imaginary part is
198 * selected. The real part is used if ISRT = 0, and the
199 * imaginary part if ISRT = 1.
200 * Not referenced if COMP = .FALSE.
201 *
202 * ISRT (input) INTEGER
203 * When COMP = .TRUE., ISRT describes how ISLCT is used to
204 * choose a subset of the spectrum.
205 * Not referenced if COMP = .FALSE.
206 *
207 * RESULT (output) DOUBLE PRECISION array, dimension (17)
208 * The values computed by the 17 tests described above.
209 * The values are currently limited to 1/ulp, to avoid
210 * overflow.
211 *
212 * WORK (workspace) COMPLEX*16 array, dimension (2*N*N)
213 *
214 * LWORK (input) INTEGER
215 * The number of entries in WORK to be passed to ZGEESX. This
216 * must be at least 2*N, and N*(N+1)/2 if tests 14--16 are to
217 * be performed.
218 *
219 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
220 *
221 * BWORK (workspace) LOGICAL array, dimension (N)
222 *
223 * INFO (output) INTEGER
224 * If 0, successful exit.
225 * If <0, input parameter -INFO had an incorrect value.
226 * If >0, ZGEESX returned an error code, the absolute
227 * value of which is returned.
228 *
229 * =====================================================================
230 *
231 * .. Parameters ..
232 COMPLEX*16 CZERO, CONE
233 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
234 $ CONE = ( 1.0D+0, 0.0D+0 ) )
235 DOUBLE PRECISION ZERO, ONE
236 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
237 DOUBLE PRECISION EPSIN
238 PARAMETER ( EPSIN = 5.9605D-8 )
239 * ..
240 * .. Local Scalars ..
241 CHARACTER SORT
242 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, RSUB,
243 $ SDIM, SDIM1
244 DOUBLE PRECISION ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
245 $ SMLNUM, TOL, TOLIN, ULP, ULPINV, V, VRICMP,
246 $ VRIMIN, WNORM
247 COMPLEX*16 CTMP
248 * ..
249 * .. Local Arrays ..
250 INTEGER IPNT( 20 )
251 * ..
252 * .. External Functions ..
253 LOGICAL ZSLECT
254 DOUBLE PRECISION DLAMCH, ZLANGE
255 EXTERNAL ZSLECT, DLAMCH, ZLANGE
256 * ..
257 * .. External Subroutines ..
258 EXTERNAL XERBLA, ZCOPY, ZGEESX, ZGEMM, ZLACPY, ZUNT01
259 * ..
260 * .. Intrinsic Functions ..
261 INTRINSIC ABS, DBLE, DIMAG, MAX, MIN
262 * ..
263 * .. Arrays in Common ..
264 LOGICAL SELVAL( 20 )
265 DOUBLE PRECISION SELWI( 20 ), SELWR( 20 )
266 * ..
267 * .. Scalars in Common ..
268 INTEGER SELDIM, SELOPT
269 * ..
270 * .. Common blocks ..
271 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
272 * ..
273 * .. Executable Statements ..
274 *
275 * Check for errors
276 *
277 INFO = 0
278 IF( THRESH.LT.ZERO ) THEN
279 INFO = -3
280 ELSE IF( NOUNIT.LE.0 ) THEN
281 INFO = -5
282 ELSE IF( N.LT.0 ) THEN
283 INFO = -6
284 ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN
285 INFO = -8
286 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.N ) THEN
287 INFO = -15
288 ELSE IF( LWORK.LT.2*N ) THEN
289 INFO = -24
290 END IF
291 *
292 IF( INFO.NE.0 ) THEN
293 CALL XERBLA( 'ZGET24', -INFO )
294 RETURN
295 END IF
296 *
297 * Quick return if nothing to do
298 *
299 DO 10 I = 1, 17
300 RESULT( I ) = -ONE
301 10 CONTINUE
302 *
303 IF( N.EQ.0 )
304 $ RETURN
305 *
306 * Important constants
307 *
308 SMLNUM = DLAMCH( 'Safe minimum' )
309 ULP = DLAMCH( 'Precision' )
310 ULPINV = ONE / ULP
311 *
312 * Perform tests (1)-(13)
313 *
314 SELOPT = 0
315 DO 90 ISORT = 0, 1
316 IF( ISORT.EQ.0 ) THEN
317 SORT = 'N'
318 RSUB = 0
319 ELSE
320 SORT = 'S'
321 RSUB = 6
322 END IF
323 *
324 * Compute Schur form and Schur vectors, and test them
325 *
326 CALL ZLACPY( 'F', N, N, A, LDA, H, LDA )
327 CALL ZGEESX( 'V', SORT, ZSLECT, 'N', N, H, LDA, SDIM, W, VS,
328 $ LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK,
329 $ IINFO )
330 IF( IINFO.NE.0 ) THEN
331 RESULT( 1+RSUB ) = ULPINV
332 IF( JTYPE.NE.22 ) THEN
333 WRITE( NOUNIT, FMT = 9998 )'ZGEESX1', IINFO, N, JTYPE,
334 $ ISEED
335 ELSE
336 WRITE( NOUNIT, FMT = 9999 )'ZGEESX1', IINFO, N,
337 $ ISEED( 1 )
338 END IF
339 INFO = ABS( IINFO )
340 RETURN
341 END IF
342 IF( ISORT.EQ.0 ) THEN
343 CALL ZCOPY( N, W, 1, WTMP, 1 )
344 END IF
345 *
346 * Do Test (1) or Test (7)
347 *
348 RESULT( 1+RSUB ) = ZERO
349 DO 30 J = 1, N - 1
350 DO 20 I = J + 1, N
351 IF( H( I, J ).NE.CZERO )
352 $ RESULT( 1+RSUB ) = ULPINV
353 20 CONTINUE
354 30 CONTINUE
355 *
356 * Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
357 *
358 * Copy A to VS1, used as workspace
359 *
360 CALL ZLACPY( ' ', N, N, A, LDA, VS1, LDVS )
361 *
362 * Compute Q*H and store in HT.
363 *
364 CALL ZGEMM( 'No transpose', 'No transpose', N, N, N, CONE, VS,
365 $ LDVS, H, LDA, CZERO, HT, LDA )
366 *
367 * Compute A - Q*H*Q'
368 *
369 CALL ZGEMM( 'No transpose', 'Conjugate transpose', N, N, N,
370 $ -CONE, HT, LDA, VS, LDVS, CONE, VS1, LDVS )
371 *
372 ANORM = MAX( ZLANGE( '1', N, N, A, LDA, RWORK ), SMLNUM )
373 WNORM = ZLANGE( '1', N, N, VS1, LDVS, RWORK )
374 *
375 IF( ANORM.GT.WNORM ) THEN
376 RESULT( 2+RSUB ) = ( WNORM / ANORM ) / ( N*ULP )
377 ELSE
378 IF( ANORM.LT.ONE ) THEN
379 RESULT( 2+RSUB ) = ( MIN( WNORM, N*ANORM ) / ANORM ) /
380 $ ( N*ULP )
381 ELSE
382 RESULT( 2+RSUB ) = MIN( WNORM / ANORM, DBLE( N ) ) /
383 $ ( N*ULP )
384 END IF
385 END IF
386 *
387 * Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP )
388 *
389 CALL ZUNT01( 'Columns', N, N, VS, LDVS, WORK, LWORK, RWORK,
390 $ RESULT( 3+RSUB ) )
391 *
392 * Do Test (4) or Test (10)
393 *
394 RESULT( 4+RSUB ) = ZERO
395 DO 40 I = 1, N
396 IF( H( I, I ).NE.W( I ) )
397 $ RESULT( 4+RSUB ) = ULPINV
398 40 CONTINUE
399 *
400 * Do Test (5) or Test (11)
401 *
402 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
403 CALL ZGEESX( 'N', SORT, ZSLECT, 'N', N, HT, LDA, SDIM, WT, VS,
404 $ LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK,
405 $ IINFO )
406 IF( IINFO.NE.0 ) THEN
407 RESULT( 5+RSUB ) = ULPINV
408 IF( JTYPE.NE.22 ) THEN
409 WRITE( NOUNIT, FMT = 9998 )'ZGEESX2', IINFO, N, JTYPE,
410 $ ISEED
411 ELSE
412 WRITE( NOUNIT, FMT = 9999 )'ZGEESX2', IINFO, N,
413 $ ISEED( 1 )
414 END IF
415 INFO = ABS( IINFO )
416 GO TO 220
417 END IF
418 *
419 RESULT( 5+RSUB ) = ZERO
420 DO 60 J = 1, N
421 DO 50 I = 1, N
422 IF( H( I, J ).NE.HT( I, J ) )
423 $ RESULT( 5+RSUB ) = ULPINV
424 50 CONTINUE
425 60 CONTINUE
426 *
427 * Do Test (6) or Test (12)
428 *
429 RESULT( 6+RSUB ) = ZERO
430 DO 70 I = 1, N
431 IF( W( I ).NE.WT( I ) )
432 $ RESULT( 6+RSUB ) = ULPINV
433 70 CONTINUE
434 *
435 * Do Test (13)
436 *
437 IF( ISORT.EQ.1 ) THEN
438 RESULT( 13 ) = ZERO
439 KNTEIG = 0
440 DO 80 I = 1, N
441 IF( ZSLECT( W( I ) ) )
442 $ KNTEIG = KNTEIG + 1
443 IF( I.LT.N ) THEN
444 IF( ZSLECT( W( I+1 ) ) .AND.
445 $ ( .NOT.ZSLECT( W( I ) ) ) )RESULT( 13 ) = ULPINV
446 END IF
447 80 CONTINUE
448 IF( SDIM.NE.KNTEIG )
449 $ RESULT( 13 ) = ULPINV
450 END IF
451 *
452 90 CONTINUE
453 *
454 * If there is enough workspace, perform tests (14) and (15)
455 * as well as (10) through (13)
456 *
457 IF( LWORK.GE.( N*( N+1 ) ) / 2 ) THEN
458 *
459 * Compute both RCONDE and RCONDV with VS
460 *
461 SORT = 'S'
462 RESULT( 14 ) = ZERO
463 RESULT( 15 ) = ZERO
464 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
465 CALL ZGEESX( 'V', SORT, ZSLECT, 'B', N, HT, LDA, SDIM1, WT,
466 $ VS1, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
467 $ BWORK, IINFO )
468 IF( IINFO.NE.0 ) THEN
469 RESULT( 14 ) = ULPINV
470 RESULT( 15 ) = ULPINV
471 IF( JTYPE.NE.22 ) THEN
472 WRITE( NOUNIT, FMT = 9998 )'ZGEESX3', IINFO, N, JTYPE,
473 $ ISEED
474 ELSE
475 WRITE( NOUNIT, FMT = 9999 )'ZGEESX3', IINFO, N,
476 $ ISEED( 1 )
477 END IF
478 INFO = ABS( IINFO )
479 GO TO 220
480 END IF
481 *
482 * Perform tests (10), (11), (12), and (13)
483 *
484 DO 110 I = 1, N
485 IF( W( I ).NE.WT( I ) )
486 $ RESULT( 10 ) = ULPINV
487 DO 100 J = 1, N
488 IF( H( I, J ).NE.HT( I, J ) )
489 $ RESULT( 11 ) = ULPINV
490 IF( VS( I, J ).NE.VS1( I, J ) )
491 $ RESULT( 12 ) = ULPINV
492 100 CONTINUE
493 110 CONTINUE
494 IF( SDIM.NE.SDIM1 )
495 $ RESULT( 13 ) = ULPINV
496 *
497 * Compute both RCONDE and RCONDV without VS, and compare
498 *
499 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
500 CALL ZGEESX( 'N', SORT, ZSLECT, 'B', N, HT, LDA, SDIM1, WT,
501 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
502 $ BWORK, IINFO )
503 IF( IINFO.NE.0 ) THEN
504 RESULT( 14 ) = ULPINV
505 RESULT( 15 ) = ULPINV
506 IF( JTYPE.NE.22 ) THEN
507 WRITE( NOUNIT, FMT = 9998 )'ZGEESX4', IINFO, N, JTYPE,
508 $ ISEED
509 ELSE
510 WRITE( NOUNIT, FMT = 9999 )'ZGEESX4', IINFO, N,
511 $ ISEED( 1 )
512 END IF
513 INFO = ABS( IINFO )
514 GO TO 220
515 END IF
516 *
517 * Perform tests (14) and (15)
518 *
519 IF( RCNDE1.NE.RCONDE )
520 $ RESULT( 14 ) = ULPINV
521 IF( RCNDV1.NE.RCONDV )
522 $ RESULT( 15 ) = ULPINV
523 *
524 * Perform tests (10), (11), (12), and (13)
525 *
526 DO 130 I = 1, N
527 IF( W( I ).NE.WT( I ) )
528 $ RESULT( 10 ) = ULPINV
529 DO 120 J = 1, N
530 IF( H( I, J ).NE.HT( I, J ) )
531 $ RESULT( 11 ) = ULPINV
532 IF( VS( I, J ).NE.VS1( I, J ) )
533 $ RESULT( 12 ) = ULPINV
534 120 CONTINUE
535 130 CONTINUE
536 IF( SDIM.NE.SDIM1 )
537 $ RESULT( 13 ) = ULPINV
538 *
539 * Compute RCONDE with VS, and compare
540 *
541 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
542 CALL ZGEESX( 'V', SORT, ZSLECT, 'E', N, HT, LDA, SDIM1, WT,
543 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
544 $ BWORK, IINFO )
545 IF( IINFO.NE.0 ) THEN
546 RESULT( 14 ) = ULPINV
547 IF( JTYPE.NE.22 ) THEN
548 WRITE( NOUNIT, FMT = 9998 )'ZGEESX5', IINFO, N, JTYPE,
549 $ ISEED
550 ELSE
551 WRITE( NOUNIT, FMT = 9999 )'ZGEESX5', IINFO, N,
552 $ ISEED( 1 )
553 END IF
554 INFO = ABS( IINFO )
555 GO TO 220
556 END IF
557 *
558 * Perform test (14)
559 *
560 IF( RCNDE1.NE.RCONDE )
561 $ RESULT( 14 ) = ULPINV
562 *
563 * Perform tests (10), (11), (12), and (13)
564 *
565 DO 150 I = 1, N
566 IF( W( I ).NE.WT( I ) )
567 $ RESULT( 10 ) = ULPINV
568 DO 140 J = 1, N
569 IF( H( I, J ).NE.HT( I, J ) )
570 $ RESULT( 11 ) = ULPINV
571 IF( VS( I, J ).NE.VS1( I, J ) )
572 $ RESULT( 12 ) = ULPINV
573 140 CONTINUE
574 150 CONTINUE
575 IF( SDIM.NE.SDIM1 )
576 $ RESULT( 13 ) = ULPINV
577 *
578 * Compute RCONDE without VS, and compare
579 *
580 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
581 CALL ZGEESX( 'N', SORT, ZSLECT, 'E', N, HT, LDA, SDIM1, WT,
582 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
583 $ BWORK, IINFO )
584 IF( IINFO.NE.0 ) THEN
585 RESULT( 14 ) = ULPINV
586 IF( JTYPE.NE.22 ) THEN
587 WRITE( NOUNIT, FMT = 9998 )'ZGEESX6', IINFO, N, JTYPE,
588 $ ISEED
589 ELSE
590 WRITE( NOUNIT, FMT = 9999 )'ZGEESX6', IINFO, N,
591 $ ISEED( 1 )
592 END IF
593 INFO = ABS( IINFO )
594 GO TO 220
595 END IF
596 *
597 * Perform test (14)
598 *
599 IF( RCNDE1.NE.RCONDE )
600 $ RESULT( 14 ) = ULPINV
601 *
602 * Perform tests (10), (11), (12), and (13)
603 *
604 DO 170 I = 1, N
605 IF( W( I ).NE.WT( I ) )
606 $ RESULT( 10 ) = ULPINV
607 DO 160 J = 1, N
608 IF( H( I, J ).NE.HT( I, J ) )
609 $ RESULT( 11 ) = ULPINV
610 IF( VS( I, J ).NE.VS1( I, J ) )
611 $ RESULT( 12 ) = ULPINV
612 160 CONTINUE
613 170 CONTINUE
614 IF( SDIM.NE.SDIM1 )
615 $ RESULT( 13 ) = ULPINV
616 *
617 * Compute RCONDV with VS, and compare
618 *
619 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
620 CALL ZGEESX( 'V', SORT, ZSLECT, 'V', N, HT, LDA, SDIM1, WT,
621 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
622 $ BWORK, IINFO )
623 IF( IINFO.NE.0 ) THEN
624 RESULT( 15 ) = ULPINV
625 IF( JTYPE.NE.22 ) THEN
626 WRITE( NOUNIT, FMT = 9998 )'ZGEESX7', IINFO, N, JTYPE,
627 $ ISEED
628 ELSE
629 WRITE( NOUNIT, FMT = 9999 )'ZGEESX7', IINFO, N,
630 $ ISEED( 1 )
631 END IF
632 INFO = ABS( IINFO )
633 GO TO 220
634 END IF
635 *
636 * Perform test (15)
637 *
638 IF( RCNDV1.NE.RCONDV )
639 $ RESULT( 15 ) = ULPINV
640 *
641 * Perform tests (10), (11), (12), and (13)
642 *
643 DO 190 I = 1, N
644 IF( W( I ).NE.WT( I ) )
645 $ RESULT( 10 ) = ULPINV
646 DO 180 J = 1, N
647 IF( H( I, J ).NE.HT( I, J ) )
648 $ RESULT( 11 ) = ULPINV
649 IF( VS( I, J ).NE.VS1( I, J ) )
650 $ RESULT( 12 ) = ULPINV
651 180 CONTINUE
652 190 CONTINUE
653 IF( SDIM.NE.SDIM1 )
654 $ RESULT( 13 ) = ULPINV
655 *
656 * Compute RCONDV without VS, and compare
657 *
658 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
659 CALL ZGEESX( 'N', SORT, ZSLECT, 'V', N, HT, LDA, SDIM1, WT,
660 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
661 $ BWORK, IINFO )
662 IF( IINFO.NE.0 ) THEN
663 RESULT( 15 ) = ULPINV
664 IF( JTYPE.NE.22 ) THEN
665 WRITE( NOUNIT, FMT = 9998 )'ZGEESX8', IINFO, N, JTYPE,
666 $ ISEED
667 ELSE
668 WRITE( NOUNIT, FMT = 9999 )'ZGEESX8', IINFO, N,
669 $ ISEED( 1 )
670 END IF
671 INFO = ABS( IINFO )
672 GO TO 220
673 END IF
674 *
675 * Perform test (15)
676 *
677 IF( RCNDV1.NE.RCONDV )
678 $ RESULT( 15 ) = ULPINV
679 *
680 * Perform tests (10), (11), (12), and (13)
681 *
682 DO 210 I = 1, N
683 IF( W( I ).NE.WT( I ) )
684 $ RESULT( 10 ) = ULPINV
685 DO 200 J = 1, N
686 IF( H( I, J ).NE.HT( I, J ) )
687 $ RESULT( 11 ) = ULPINV
688 IF( VS( I, J ).NE.VS1( I, J ) )
689 $ RESULT( 12 ) = ULPINV
690 200 CONTINUE
691 210 CONTINUE
692 IF( SDIM.NE.SDIM1 )
693 $ RESULT( 13 ) = ULPINV
694 *
695 END IF
696 *
697 220 CONTINUE
698 *
699 * If there are precomputed reciprocal condition numbers, compare
700 * computed values with them.
701 *
702 IF( COMP ) THEN
703 *
704 * First set up SELOPT, SELDIM, SELVAL, SELWR and SELWI so that
705 * the logical function ZSLECT selects the eigenvalues specified
706 * by NSLCT, ISLCT and ISRT.
707 *
708 SELDIM = N
709 SELOPT = 1
710 EPS = MAX( ULP, EPSIN )
711 DO 230 I = 1, N
712 IPNT( I ) = I
713 SELVAL( I ) = .FALSE.
714 SELWR( I ) = DBLE( WTMP( I ) )
715 SELWI( I ) = DIMAG( WTMP( I ) )
716 230 CONTINUE
717 DO 250 I = 1, N - 1
718 KMIN = I
719 IF( ISRT.EQ.0 ) THEN
720 VRIMIN = DBLE( WTMP( I ) )
721 ELSE
722 VRIMIN = DIMAG( WTMP( I ) )
723 END IF
724 DO 240 J = I + 1, N
725 IF( ISRT.EQ.0 ) THEN
726 VRICMP = DBLE( WTMP( J ) )
727 ELSE
728 VRICMP = DIMAG( WTMP( J ) )
729 END IF
730 IF( VRICMP.LT.VRIMIN ) THEN
731 KMIN = J
732 VRIMIN = VRICMP
733 END IF
734 240 CONTINUE
735 CTMP = WTMP( KMIN )
736 WTMP( KMIN ) = WTMP( I )
737 WTMP( I ) = CTMP
738 ITMP = IPNT( I )
739 IPNT( I ) = IPNT( KMIN )
740 IPNT( KMIN ) = ITMP
741 250 CONTINUE
742 DO 260 I = 1, NSLCT
743 SELVAL( IPNT( ISLCT( I ) ) ) = .TRUE.
744 260 CONTINUE
745 *
746 * Compute condition numbers
747 *
748 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
749 CALL ZGEESX( 'N', 'S', ZSLECT, 'B', N, HT, LDA, SDIM1, WT, VS1,
750 $ LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK,
751 $ IINFO )
752 IF( IINFO.NE.0 ) THEN
753 RESULT( 16 ) = ULPINV
754 RESULT( 17 ) = ULPINV
755 WRITE( NOUNIT, FMT = 9999 )'ZGEESX9', IINFO, N, ISEED( 1 )
756 INFO = ABS( IINFO )
757 GO TO 270
758 END IF
759 *
760 * Compare condition number for average of selected eigenvalues
761 * taking its condition number into account
762 *
763 ANORM = ZLANGE( '1', N, N, A, LDA, RWORK )
764 V = MAX( DBLE( N )*EPS*ANORM, SMLNUM )
765 IF( ANORM.EQ.ZERO )
766 $ V = ONE
767 IF( V.GT.RCONDV ) THEN
768 TOL = ONE
769 ELSE
770 TOL = V / RCONDV
771 END IF
772 IF( V.GT.RCDVIN ) THEN
773 TOLIN = ONE
774 ELSE
775 TOLIN = V / RCDVIN
776 END IF
777 TOL = MAX( TOL, SMLNUM / EPS )
778 TOLIN = MAX( TOLIN, SMLNUM / EPS )
779 IF( EPS*( RCDEIN-TOLIN ).GT.RCONDE+TOL ) THEN
780 RESULT( 16 ) = ULPINV
781 ELSE IF( RCDEIN-TOLIN.GT.RCONDE+TOL ) THEN
782 RESULT( 16 ) = ( RCDEIN-TOLIN ) / ( RCONDE+TOL )
783 ELSE IF( RCDEIN+TOLIN.LT.EPS*( RCONDE-TOL ) ) THEN
784 RESULT( 16 ) = ULPINV
785 ELSE IF( RCDEIN+TOLIN.LT.RCONDE-TOL ) THEN
786 RESULT( 16 ) = ( RCONDE-TOL ) / ( RCDEIN+TOLIN )
787 ELSE
788 RESULT( 16 ) = ONE
789 END IF
790 *
791 * Compare condition numbers for right invariant subspace
792 * taking its condition number into account
793 *
794 IF( V.GT.RCONDV*RCONDE ) THEN
795 TOL = RCONDV
796 ELSE
797 TOL = V / RCONDE
798 END IF
799 IF( V.GT.RCDVIN*RCDEIN ) THEN
800 TOLIN = RCDVIN
801 ELSE
802 TOLIN = V / RCDEIN
803 END IF
804 TOL = MAX( TOL, SMLNUM / EPS )
805 TOLIN = MAX( TOLIN, SMLNUM / EPS )
806 IF( EPS*( RCDVIN-TOLIN ).GT.RCONDV+TOL ) THEN
807 RESULT( 17 ) = ULPINV
808 ELSE IF( RCDVIN-TOLIN.GT.RCONDV+TOL ) THEN
809 RESULT( 17 ) = ( RCDVIN-TOLIN ) / ( RCONDV+TOL )
810 ELSE IF( RCDVIN+TOLIN.LT.EPS*( RCONDV-TOL ) ) THEN
811 RESULT( 17 ) = ULPINV
812 ELSE IF( RCDVIN+TOLIN.LT.RCONDV-TOL ) THEN
813 RESULT( 17 ) = ( RCONDV-TOL ) / ( RCDVIN+TOLIN )
814 ELSE
815 RESULT( 17 ) = ONE
816 END IF
817 *
818 270 CONTINUE
819 *
820 END IF
821 *
822 9999 FORMAT( ' ZGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
823 $ I6, ', INPUT EXAMPLE NUMBER = ', I4 )
824 9998 FORMAT( ' ZGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
825 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
826 *
827 RETURN
828 *
829 * End of ZGET24
830 *
831 END
2 $ H, HT, W, WT, WTMP, VS, LDVS, VS1, RCDEIN,
3 $ RCDVIN, NSLCT, ISLCT, ISRT, RESULT, WORK,
4 $ LWORK, RWORK, BWORK, INFO )
5 *
6 * -- LAPACK test routine (version 3.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 LOGICAL COMP
12 INTEGER INFO, ISRT, JTYPE, LDA, LDVS, LWORK, N, NOUNIT,
13 $ NSLCT
14 DOUBLE PRECISION RCDEIN, RCDVIN, THRESH
15 * ..
16 * .. Array Arguments ..
17 LOGICAL BWORK( * )
18 INTEGER ISEED( 4 ), ISLCT( * )
19 DOUBLE PRECISION RESULT( 17 ), RWORK( * )
20 COMPLEX*16 A( LDA, * ), H( LDA, * ), HT( LDA, * ),
21 $ VS( LDVS, * ), VS1( LDVS, * ), W( * ),
22 $ WORK( * ), WT( * ), WTMP( * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * ZGET24 checks the nonsymmetric eigenvalue (Schur form) problem
29 * expert driver ZGEESX.
30 *
31 * If COMP = .FALSE., the first 13 of the following tests will be
32 * be performed on the input matrix A, and also tests 14 and 15
33 * if LWORK is sufficiently large.
34 * If COMP = .TRUE., all 17 test will be performed.
35 *
36 * (1) 0 if T is in Schur form, 1/ulp otherwise
37 * (no sorting of eigenvalues)
38 *
39 * (2) | A - VS T VS' | / ( n |A| ulp )
40 *
41 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
42 * form (no sorting of eigenvalues).
43 *
44 * (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
45 *
46 * (4) 0 if W are eigenvalues of T
47 * 1/ulp otherwise
48 * (no sorting of eigenvalues)
49 *
50 * (5) 0 if T(with VS) = T(without VS),
51 * 1/ulp otherwise
52 * (no sorting of eigenvalues)
53 *
54 * (6) 0 if eigenvalues(with VS) = eigenvalues(without VS),
55 * 1/ulp otherwise
56 * (no sorting of eigenvalues)
57 *
58 * (7) 0 if T is in Schur form, 1/ulp otherwise
59 * (with sorting of eigenvalues)
60 *
61 * (8) | A - VS T VS' | / ( n |A| ulp )
62 *
63 * Here VS is the matrix of Schur eigenvectors, and T is in Schur
64 * form (with sorting of eigenvalues).
65 *
66 * (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
67 *
68 * (10) 0 if W are eigenvalues of T
69 * 1/ulp otherwise
70 * If workspace sufficient, also compare W with and
71 * without reciprocal condition numbers
72 * (with sorting of eigenvalues)
73 *
74 * (11) 0 if T(with VS) = T(without VS),
75 * 1/ulp otherwise
76 * If workspace sufficient, also compare T with and without
77 * reciprocal condition numbers
78 * (with sorting of eigenvalues)
79 *
80 * (12) 0 if eigenvalues(with VS) = eigenvalues(without VS),
81 * 1/ulp otherwise
82 * If workspace sufficient, also compare VS with and without
83 * reciprocal condition numbers
84 * (with sorting of eigenvalues)
85 *
86 * (13) if sorting worked and SDIM is the number of
87 * eigenvalues which were SELECTed
88 * If workspace sufficient, also compare SDIM with and
89 * without reciprocal condition numbers
90 *
91 * (14) if RCONDE the same no matter if VS and/or RCONDV computed
92 *
93 * (15) if RCONDV the same no matter if VS and/or RCONDE computed
94 *
95 * (16) |RCONDE - RCDEIN| / cond(RCONDE)
96 *
97 * RCONDE is the reciprocal average eigenvalue condition number
98 * computed by ZGEESX and RCDEIN (the precomputed true value)
99 * is supplied as input. cond(RCONDE) is the condition number
100 * of RCONDE, and takes errors in computing RCONDE into account,
101 * so that the resulting quantity should be O(ULP). cond(RCONDE)
102 * is essentially given by norm(A)/RCONDV.
103 *
104 * (17) |RCONDV - RCDVIN| / cond(RCONDV)
105 *
106 * RCONDV is the reciprocal right invariant subspace condition
107 * number computed by ZGEESX and RCDVIN (the precomputed true
108 * value) is supplied as input. cond(RCONDV) is the condition
109 * number of RCONDV, and takes errors in computing RCONDV into
110 * account, so that the resulting quantity should be O(ULP).
111 * cond(RCONDV) is essentially given by norm(A)/RCONDE.
112 *
113 * Arguments
114 * =========
115 *
116 * COMP (input) LOGICAL
117 * COMP describes which input tests to perform:
118 * = .FALSE. if the computed condition numbers are not to
119 * be tested against RCDVIN and RCDEIN
120 * = .TRUE. if they are to be compared
121 *
122 * JTYPE (input) INTEGER
123 * Type of input matrix. Used to label output if error occurs.
124 *
125 * ISEED (input) INTEGER array, dimension (4)
126 * If COMP = .FALSE., the random number generator seed
127 * used to produce matrix.
128 * If COMP = .TRUE., ISEED(1) = the number of the example.
129 * Used to label output if error occurs.
130 *
131 * THRESH (input) DOUBLE PRECISION
132 * A test will count as "failed" if the "error", computed as
133 * described above, exceeds THRESH. Note that the error
134 * is scaled to be O(1), so THRESH should be a reasonably
135 * small multiple of 1, e.g., 10 or 100. In particular,
136 * it should not depend on the precision (single vs. double)
137 * or the size of the matrix. It must be at least zero.
138 *
139 * NOUNIT (input) INTEGER
140 * The FORTRAN unit number for printing out error messages
141 * (e.g., if a routine returns INFO not equal to 0.)
142 *
143 * N (input) INTEGER
144 * The dimension of A. N must be at least 0.
145 *
146 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
147 * Used to hold the matrix whose eigenvalues are to be
148 * computed.
149 *
150 * LDA (input) INTEGER
151 * The leading dimension of A, and H. LDA must be at
152 * least 1 and at least N.
153 *
154 * H (workspace) COMPLEX*16 array, dimension (LDA, N)
155 * Another copy of the test matrix A, modified by ZGEESX.
156 *
157 * HT (workspace) COMPLEX*16 array, dimension (LDA, N)
158 * Yet another copy of the test matrix A, modified by ZGEESX.
159 *
160 * W (workspace) COMPLEX*16 array, dimension (N)
161 * The computed eigenvalues of A.
162 *
163 * WT (workspace) COMPLEX*16 array, dimension (N)
164 * Like W, this array contains the eigenvalues of A,
165 * but those computed when ZGEESX only computes a partial
166 * eigendecomposition, i.e. not Schur vectors
167 *
168 * WTMP (workspace) COMPLEX*16 array, dimension (N)
169 * Like W, this array contains the eigenvalues of A,
170 * but sorted by increasing real or imaginary part.
171 *
172 * VS (workspace) COMPLEX*16 array, dimension (LDVS, N)
173 * VS holds the computed Schur vectors.
174 *
175 * LDVS (input) INTEGER
176 * Leading dimension of VS. Must be at least max(1, N).
177 *
178 * VS1 (workspace) COMPLEX*16 array, dimension (LDVS, N)
179 * VS1 holds another copy of the computed Schur vectors.
180 *
181 * RCDEIN (input) DOUBLE PRECISION
182 * When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
183 * condition number for the average of selected eigenvalues.
184 *
185 * RCDVIN (input) DOUBLE PRECISION
186 * When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
187 * condition number for the selected right invariant subspace.
188 *
189 * NSLCT (input) INTEGER
190 * When COMP = .TRUE. the number of selected eigenvalues
191 * corresponding to the precomputed values RCDEIN and RCDVIN.
192 *
193 * ISLCT (input) INTEGER array, dimension (NSLCT)
194 * When COMP = .TRUE. ISLCT selects the eigenvalues of the
195 * input matrix corresponding to the precomputed values RCDEIN
196 * and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
197 * eigenvalue with the J-th largest real or imaginary part is
198 * selected. The real part is used if ISRT = 0, and the
199 * imaginary part if ISRT = 1.
200 * Not referenced if COMP = .FALSE.
201 *
202 * ISRT (input) INTEGER
203 * When COMP = .TRUE., ISRT describes how ISLCT is used to
204 * choose a subset of the spectrum.
205 * Not referenced if COMP = .FALSE.
206 *
207 * RESULT (output) DOUBLE PRECISION array, dimension (17)
208 * The values computed by the 17 tests described above.
209 * The values are currently limited to 1/ulp, to avoid
210 * overflow.
211 *
212 * WORK (workspace) COMPLEX*16 array, dimension (2*N*N)
213 *
214 * LWORK (input) INTEGER
215 * The number of entries in WORK to be passed to ZGEESX. This
216 * must be at least 2*N, and N*(N+1)/2 if tests 14--16 are to
217 * be performed.
218 *
219 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
220 *
221 * BWORK (workspace) LOGICAL array, dimension (N)
222 *
223 * INFO (output) INTEGER
224 * If 0, successful exit.
225 * If <0, input parameter -INFO had an incorrect value.
226 * If >0, ZGEESX returned an error code, the absolute
227 * value of which is returned.
228 *
229 * =====================================================================
230 *
231 * .. Parameters ..
232 COMPLEX*16 CZERO, CONE
233 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
234 $ CONE = ( 1.0D+0, 0.0D+0 ) )
235 DOUBLE PRECISION ZERO, ONE
236 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
237 DOUBLE PRECISION EPSIN
238 PARAMETER ( EPSIN = 5.9605D-8 )
239 * ..
240 * .. Local Scalars ..
241 CHARACTER SORT
242 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, RSUB,
243 $ SDIM, SDIM1
244 DOUBLE PRECISION ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
245 $ SMLNUM, TOL, TOLIN, ULP, ULPINV, V, VRICMP,
246 $ VRIMIN, WNORM
247 COMPLEX*16 CTMP
248 * ..
249 * .. Local Arrays ..
250 INTEGER IPNT( 20 )
251 * ..
252 * .. External Functions ..
253 LOGICAL ZSLECT
254 DOUBLE PRECISION DLAMCH, ZLANGE
255 EXTERNAL ZSLECT, DLAMCH, ZLANGE
256 * ..
257 * .. External Subroutines ..
258 EXTERNAL XERBLA, ZCOPY, ZGEESX, ZGEMM, ZLACPY, ZUNT01
259 * ..
260 * .. Intrinsic Functions ..
261 INTRINSIC ABS, DBLE, DIMAG, MAX, MIN
262 * ..
263 * .. Arrays in Common ..
264 LOGICAL SELVAL( 20 )
265 DOUBLE PRECISION SELWI( 20 ), SELWR( 20 )
266 * ..
267 * .. Scalars in Common ..
268 INTEGER SELDIM, SELOPT
269 * ..
270 * .. Common blocks ..
271 COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
272 * ..
273 * .. Executable Statements ..
274 *
275 * Check for errors
276 *
277 INFO = 0
278 IF( THRESH.LT.ZERO ) THEN
279 INFO = -3
280 ELSE IF( NOUNIT.LE.0 ) THEN
281 INFO = -5
282 ELSE IF( N.LT.0 ) THEN
283 INFO = -6
284 ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN
285 INFO = -8
286 ELSE IF( LDVS.LT.1 .OR. LDVS.LT.N ) THEN
287 INFO = -15
288 ELSE IF( LWORK.LT.2*N ) THEN
289 INFO = -24
290 END IF
291 *
292 IF( INFO.NE.0 ) THEN
293 CALL XERBLA( 'ZGET24', -INFO )
294 RETURN
295 END IF
296 *
297 * Quick return if nothing to do
298 *
299 DO 10 I = 1, 17
300 RESULT( I ) = -ONE
301 10 CONTINUE
302 *
303 IF( N.EQ.0 )
304 $ RETURN
305 *
306 * Important constants
307 *
308 SMLNUM = DLAMCH( 'Safe minimum' )
309 ULP = DLAMCH( 'Precision' )
310 ULPINV = ONE / ULP
311 *
312 * Perform tests (1)-(13)
313 *
314 SELOPT = 0
315 DO 90 ISORT = 0, 1
316 IF( ISORT.EQ.0 ) THEN
317 SORT = 'N'
318 RSUB = 0
319 ELSE
320 SORT = 'S'
321 RSUB = 6
322 END IF
323 *
324 * Compute Schur form and Schur vectors, and test them
325 *
326 CALL ZLACPY( 'F', N, N, A, LDA, H, LDA )
327 CALL ZGEESX( 'V', SORT, ZSLECT, 'N', N, H, LDA, SDIM, W, VS,
328 $ LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK,
329 $ IINFO )
330 IF( IINFO.NE.0 ) THEN
331 RESULT( 1+RSUB ) = ULPINV
332 IF( JTYPE.NE.22 ) THEN
333 WRITE( NOUNIT, FMT = 9998 )'ZGEESX1', IINFO, N, JTYPE,
334 $ ISEED
335 ELSE
336 WRITE( NOUNIT, FMT = 9999 )'ZGEESX1', IINFO, N,
337 $ ISEED( 1 )
338 END IF
339 INFO = ABS( IINFO )
340 RETURN
341 END IF
342 IF( ISORT.EQ.0 ) THEN
343 CALL ZCOPY( N, W, 1, WTMP, 1 )
344 END IF
345 *
346 * Do Test (1) or Test (7)
347 *
348 RESULT( 1+RSUB ) = ZERO
349 DO 30 J = 1, N - 1
350 DO 20 I = J + 1, N
351 IF( H( I, J ).NE.CZERO )
352 $ RESULT( 1+RSUB ) = ULPINV
353 20 CONTINUE
354 30 CONTINUE
355 *
356 * Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
357 *
358 * Copy A to VS1, used as workspace
359 *
360 CALL ZLACPY( ' ', N, N, A, LDA, VS1, LDVS )
361 *
362 * Compute Q*H and store in HT.
363 *
364 CALL ZGEMM( 'No transpose', 'No transpose', N, N, N, CONE, VS,
365 $ LDVS, H, LDA, CZERO, HT, LDA )
366 *
367 * Compute A - Q*H*Q'
368 *
369 CALL ZGEMM( 'No transpose', 'Conjugate transpose', N, N, N,
370 $ -CONE, HT, LDA, VS, LDVS, CONE, VS1, LDVS )
371 *
372 ANORM = MAX( ZLANGE( '1', N, N, A, LDA, RWORK ), SMLNUM )
373 WNORM = ZLANGE( '1', N, N, VS1, LDVS, RWORK )
374 *
375 IF( ANORM.GT.WNORM ) THEN
376 RESULT( 2+RSUB ) = ( WNORM / ANORM ) / ( N*ULP )
377 ELSE
378 IF( ANORM.LT.ONE ) THEN
379 RESULT( 2+RSUB ) = ( MIN( WNORM, N*ANORM ) / ANORM ) /
380 $ ( N*ULP )
381 ELSE
382 RESULT( 2+RSUB ) = MIN( WNORM / ANORM, DBLE( N ) ) /
383 $ ( N*ULP )
384 END IF
385 END IF
386 *
387 * Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP )
388 *
389 CALL ZUNT01( 'Columns', N, N, VS, LDVS, WORK, LWORK, RWORK,
390 $ RESULT( 3+RSUB ) )
391 *
392 * Do Test (4) or Test (10)
393 *
394 RESULT( 4+RSUB ) = ZERO
395 DO 40 I = 1, N
396 IF( H( I, I ).NE.W( I ) )
397 $ RESULT( 4+RSUB ) = ULPINV
398 40 CONTINUE
399 *
400 * Do Test (5) or Test (11)
401 *
402 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
403 CALL ZGEESX( 'N', SORT, ZSLECT, 'N', N, HT, LDA, SDIM, WT, VS,
404 $ LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK,
405 $ IINFO )
406 IF( IINFO.NE.0 ) THEN
407 RESULT( 5+RSUB ) = ULPINV
408 IF( JTYPE.NE.22 ) THEN
409 WRITE( NOUNIT, FMT = 9998 )'ZGEESX2', IINFO, N, JTYPE,
410 $ ISEED
411 ELSE
412 WRITE( NOUNIT, FMT = 9999 )'ZGEESX2', IINFO, N,
413 $ ISEED( 1 )
414 END IF
415 INFO = ABS( IINFO )
416 GO TO 220
417 END IF
418 *
419 RESULT( 5+RSUB ) = ZERO
420 DO 60 J = 1, N
421 DO 50 I = 1, N
422 IF( H( I, J ).NE.HT( I, J ) )
423 $ RESULT( 5+RSUB ) = ULPINV
424 50 CONTINUE
425 60 CONTINUE
426 *
427 * Do Test (6) or Test (12)
428 *
429 RESULT( 6+RSUB ) = ZERO
430 DO 70 I = 1, N
431 IF( W( I ).NE.WT( I ) )
432 $ RESULT( 6+RSUB ) = ULPINV
433 70 CONTINUE
434 *
435 * Do Test (13)
436 *
437 IF( ISORT.EQ.1 ) THEN
438 RESULT( 13 ) = ZERO
439 KNTEIG = 0
440 DO 80 I = 1, N
441 IF( ZSLECT( W( I ) ) )
442 $ KNTEIG = KNTEIG + 1
443 IF( I.LT.N ) THEN
444 IF( ZSLECT( W( I+1 ) ) .AND.
445 $ ( .NOT.ZSLECT( W( I ) ) ) )RESULT( 13 ) = ULPINV
446 END IF
447 80 CONTINUE
448 IF( SDIM.NE.KNTEIG )
449 $ RESULT( 13 ) = ULPINV
450 END IF
451 *
452 90 CONTINUE
453 *
454 * If there is enough workspace, perform tests (14) and (15)
455 * as well as (10) through (13)
456 *
457 IF( LWORK.GE.( N*( N+1 ) ) / 2 ) THEN
458 *
459 * Compute both RCONDE and RCONDV with VS
460 *
461 SORT = 'S'
462 RESULT( 14 ) = ZERO
463 RESULT( 15 ) = ZERO
464 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
465 CALL ZGEESX( 'V', SORT, ZSLECT, 'B', N, HT, LDA, SDIM1, WT,
466 $ VS1, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
467 $ BWORK, IINFO )
468 IF( IINFO.NE.0 ) THEN
469 RESULT( 14 ) = ULPINV
470 RESULT( 15 ) = ULPINV
471 IF( JTYPE.NE.22 ) THEN
472 WRITE( NOUNIT, FMT = 9998 )'ZGEESX3', IINFO, N, JTYPE,
473 $ ISEED
474 ELSE
475 WRITE( NOUNIT, FMT = 9999 )'ZGEESX3', IINFO, N,
476 $ ISEED( 1 )
477 END IF
478 INFO = ABS( IINFO )
479 GO TO 220
480 END IF
481 *
482 * Perform tests (10), (11), (12), and (13)
483 *
484 DO 110 I = 1, N
485 IF( W( I ).NE.WT( I ) )
486 $ RESULT( 10 ) = ULPINV
487 DO 100 J = 1, N
488 IF( H( I, J ).NE.HT( I, J ) )
489 $ RESULT( 11 ) = ULPINV
490 IF( VS( I, J ).NE.VS1( I, J ) )
491 $ RESULT( 12 ) = ULPINV
492 100 CONTINUE
493 110 CONTINUE
494 IF( SDIM.NE.SDIM1 )
495 $ RESULT( 13 ) = ULPINV
496 *
497 * Compute both RCONDE and RCONDV without VS, and compare
498 *
499 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
500 CALL ZGEESX( 'N', SORT, ZSLECT, 'B', N, HT, LDA, SDIM1, WT,
501 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
502 $ BWORK, IINFO )
503 IF( IINFO.NE.0 ) THEN
504 RESULT( 14 ) = ULPINV
505 RESULT( 15 ) = ULPINV
506 IF( JTYPE.NE.22 ) THEN
507 WRITE( NOUNIT, FMT = 9998 )'ZGEESX4', IINFO, N, JTYPE,
508 $ ISEED
509 ELSE
510 WRITE( NOUNIT, FMT = 9999 )'ZGEESX4', IINFO, N,
511 $ ISEED( 1 )
512 END IF
513 INFO = ABS( IINFO )
514 GO TO 220
515 END IF
516 *
517 * Perform tests (14) and (15)
518 *
519 IF( RCNDE1.NE.RCONDE )
520 $ RESULT( 14 ) = ULPINV
521 IF( RCNDV1.NE.RCONDV )
522 $ RESULT( 15 ) = ULPINV
523 *
524 * Perform tests (10), (11), (12), and (13)
525 *
526 DO 130 I = 1, N
527 IF( W( I ).NE.WT( I ) )
528 $ RESULT( 10 ) = ULPINV
529 DO 120 J = 1, N
530 IF( H( I, J ).NE.HT( I, J ) )
531 $ RESULT( 11 ) = ULPINV
532 IF( VS( I, J ).NE.VS1( I, J ) )
533 $ RESULT( 12 ) = ULPINV
534 120 CONTINUE
535 130 CONTINUE
536 IF( SDIM.NE.SDIM1 )
537 $ RESULT( 13 ) = ULPINV
538 *
539 * Compute RCONDE with VS, and compare
540 *
541 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
542 CALL ZGEESX( 'V', SORT, ZSLECT, 'E', N, HT, LDA, SDIM1, WT,
543 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
544 $ BWORK, IINFO )
545 IF( IINFO.NE.0 ) THEN
546 RESULT( 14 ) = ULPINV
547 IF( JTYPE.NE.22 ) THEN
548 WRITE( NOUNIT, FMT = 9998 )'ZGEESX5', IINFO, N, JTYPE,
549 $ ISEED
550 ELSE
551 WRITE( NOUNIT, FMT = 9999 )'ZGEESX5', IINFO, N,
552 $ ISEED( 1 )
553 END IF
554 INFO = ABS( IINFO )
555 GO TO 220
556 END IF
557 *
558 * Perform test (14)
559 *
560 IF( RCNDE1.NE.RCONDE )
561 $ RESULT( 14 ) = ULPINV
562 *
563 * Perform tests (10), (11), (12), and (13)
564 *
565 DO 150 I = 1, N
566 IF( W( I ).NE.WT( I ) )
567 $ RESULT( 10 ) = ULPINV
568 DO 140 J = 1, N
569 IF( H( I, J ).NE.HT( I, J ) )
570 $ RESULT( 11 ) = ULPINV
571 IF( VS( I, J ).NE.VS1( I, J ) )
572 $ RESULT( 12 ) = ULPINV
573 140 CONTINUE
574 150 CONTINUE
575 IF( SDIM.NE.SDIM1 )
576 $ RESULT( 13 ) = ULPINV
577 *
578 * Compute RCONDE without VS, and compare
579 *
580 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
581 CALL ZGEESX( 'N', SORT, ZSLECT, 'E', N, HT, LDA, SDIM1, WT,
582 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
583 $ BWORK, IINFO )
584 IF( IINFO.NE.0 ) THEN
585 RESULT( 14 ) = ULPINV
586 IF( JTYPE.NE.22 ) THEN
587 WRITE( NOUNIT, FMT = 9998 )'ZGEESX6', IINFO, N, JTYPE,
588 $ ISEED
589 ELSE
590 WRITE( NOUNIT, FMT = 9999 )'ZGEESX6', IINFO, N,
591 $ ISEED( 1 )
592 END IF
593 INFO = ABS( IINFO )
594 GO TO 220
595 END IF
596 *
597 * Perform test (14)
598 *
599 IF( RCNDE1.NE.RCONDE )
600 $ RESULT( 14 ) = ULPINV
601 *
602 * Perform tests (10), (11), (12), and (13)
603 *
604 DO 170 I = 1, N
605 IF( W( I ).NE.WT( I ) )
606 $ RESULT( 10 ) = ULPINV
607 DO 160 J = 1, N
608 IF( H( I, J ).NE.HT( I, J ) )
609 $ RESULT( 11 ) = ULPINV
610 IF( VS( I, J ).NE.VS1( I, J ) )
611 $ RESULT( 12 ) = ULPINV
612 160 CONTINUE
613 170 CONTINUE
614 IF( SDIM.NE.SDIM1 )
615 $ RESULT( 13 ) = ULPINV
616 *
617 * Compute RCONDV with VS, and compare
618 *
619 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
620 CALL ZGEESX( 'V', SORT, ZSLECT, 'V', N, HT, LDA, SDIM1, WT,
621 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
622 $ BWORK, IINFO )
623 IF( IINFO.NE.0 ) THEN
624 RESULT( 15 ) = ULPINV
625 IF( JTYPE.NE.22 ) THEN
626 WRITE( NOUNIT, FMT = 9998 )'ZGEESX7', IINFO, N, JTYPE,
627 $ ISEED
628 ELSE
629 WRITE( NOUNIT, FMT = 9999 )'ZGEESX7', IINFO, N,
630 $ ISEED( 1 )
631 END IF
632 INFO = ABS( IINFO )
633 GO TO 220
634 END IF
635 *
636 * Perform test (15)
637 *
638 IF( RCNDV1.NE.RCONDV )
639 $ RESULT( 15 ) = ULPINV
640 *
641 * Perform tests (10), (11), (12), and (13)
642 *
643 DO 190 I = 1, N
644 IF( W( I ).NE.WT( I ) )
645 $ RESULT( 10 ) = ULPINV
646 DO 180 J = 1, N
647 IF( H( I, J ).NE.HT( I, J ) )
648 $ RESULT( 11 ) = ULPINV
649 IF( VS( I, J ).NE.VS1( I, J ) )
650 $ RESULT( 12 ) = ULPINV
651 180 CONTINUE
652 190 CONTINUE
653 IF( SDIM.NE.SDIM1 )
654 $ RESULT( 13 ) = ULPINV
655 *
656 * Compute RCONDV without VS, and compare
657 *
658 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
659 CALL ZGEESX( 'N', SORT, ZSLECT, 'V', N, HT, LDA, SDIM1, WT,
660 $ VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK, RWORK,
661 $ BWORK, IINFO )
662 IF( IINFO.NE.0 ) THEN
663 RESULT( 15 ) = ULPINV
664 IF( JTYPE.NE.22 ) THEN
665 WRITE( NOUNIT, FMT = 9998 )'ZGEESX8', IINFO, N, JTYPE,
666 $ ISEED
667 ELSE
668 WRITE( NOUNIT, FMT = 9999 )'ZGEESX8', IINFO, N,
669 $ ISEED( 1 )
670 END IF
671 INFO = ABS( IINFO )
672 GO TO 220
673 END IF
674 *
675 * Perform test (15)
676 *
677 IF( RCNDV1.NE.RCONDV )
678 $ RESULT( 15 ) = ULPINV
679 *
680 * Perform tests (10), (11), (12), and (13)
681 *
682 DO 210 I = 1, N
683 IF( W( I ).NE.WT( I ) )
684 $ RESULT( 10 ) = ULPINV
685 DO 200 J = 1, N
686 IF( H( I, J ).NE.HT( I, J ) )
687 $ RESULT( 11 ) = ULPINV
688 IF( VS( I, J ).NE.VS1( I, J ) )
689 $ RESULT( 12 ) = ULPINV
690 200 CONTINUE
691 210 CONTINUE
692 IF( SDIM.NE.SDIM1 )
693 $ RESULT( 13 ) = ULPINV
694 *
695 END IF
696 *
697 220 CONTINUE
698 *
699 * If there are precomputed reciprocal condition numbers, compare
700 * computed values with them.
701 *
702 IF( COMP ) THEN
703 *
704 * First set up SELOPT, SELDIM, SELVAL, SELWR and SELWI so that
705 * the logical function ZSLECT selects the eigenvalues specified
706 * by NSLCT, ISLCT and ISRT.
707 *
708 SELDIM = N
709 SELOPT = 1
710 EPS = MAX( ULP, EPSIN )
711 DO 230 I = 1, N
712 IPNT( I ) = I
713 SELVAL( I ) = .FALSE.
714 SELWR( I ) = DBLE( WTMP( I ) )
715 SELWI( I ) = DIMAG( WTMP( I ) )
716 230 CONTINUE
717 DO 250 I = 1, N - 1
718 KMIN = I
719 IF( ISRT.EQ.0 ) THEN
720 VRIMIN = DBLE( WTMP( I ) )
721 ELSE
722 VRIMIN = DIMAG( WTMP( I ) )
723 END IF
724 DO 240 J = I + 1, N
725 IF( ISRT.EQ.0 ) THEN
726 VRICMP = DBLE( WTMP( J ) )
727 ELSE
728 VRICMP = DIMAG( WTMP( J ) )
729 END IF
730 IF( VRICMP.LT.VRIMIN ) THEN
731 KMIN = J
732 VRIMIN = VRICMP
733 END IF
734 240 CONTINUE
735 CTMP = WTMP( KMIN )
736 WTMP( KMIN ) = WTMP( I )
737 WTMP( I ) = CTMP
738 ITMP = IPNT( I )
739 IPNT( I ) = IPNT( KMIN )
740 IPNT( KMIN ) = ITMP
741 250 CONTINUE
742 DO 260 I = 1, NSLCT
743 SELVAL( IPNT( ISLCT( I ) ) ) = .TRUE.
744 260 CONTINUE
745 *
746 * Compute condition numbers
747 *
748 CALL ZLACPY( 'F', N, N, A, LDA, HT, LDA )
749 CALL ZGEESX( 'N', 'S', ZSLECT, 'B', N, HT, LDA, SDIM1, WT, VS1,
750 $ LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK,
751 $ IINFO )
752 IF( IINFO.NE.0 ) THEN
753 RESULT( 16 ) = ULPINV
754 RESULT( 17 ) = ULPINV
755 WRITE( NOUNIT, FMT = 9999 )'ZGEESX9', IINFO, N, ISEED( 1 )
756 INFO = ABS( IINFO )
757 GO TO 270
758 END IF
759 *
760 * Compare condition number for average of selected eigenvalues
761 * taking its condition number into account
762 *
763 ANORM = ZLANGE( '1', N, N, A, LDA, RWORK )
764 V = MAX( DBLE( N )*EPS*ANORM, SMLNUM )
765 IF( ANORM.EQ.ZERO )
766 $ V = ONE
767 IF( V.GT.RCONDV ) THEN
768 TOL = ONE
769 ELSE
770 TOL = V / RCONDV
771 END IF
772 IF( V.GT.RCDVIN ) THEN
773 TOLIN = ONE
774 ELSE
775 TOLIN = V / RCDVIN
776 END IF
777 TOL = MAX( TOL, SMLNUM / EPS )
778 TOLIN = MAX( TOLIN, SMLNUM / EPS )
779 IF( EPS*( RCDEIN-TOLIN ).GT.RCONDE+TOL ) THEN
780 RESULT( 16 ) = ULPINV
781 ELSE IF( RCDEIN-TOLIN.GT.RCONDE+TOL ) THEN
782 RESULT( 16 ) = ( RCDEIN-TOLIN ) / ( RCONDE+TOL )
783 ELSE IF( RCDEIN+TOLIN.LT.EPS*( RCONDE-TOL ) ) THEN
784 RESULT( 16 ) = ULPINV
785 ELSE IF( RCDEIN+TOLIN.LT.RCONDE-TOL ) THEN
786 RESULT( 16 ) = ( RCONDE-TOL ) / ( RCDEIN+TOLIN )
787 ELSE
788 RESULT( 16 ) = ONE
789 END IF
790 *
791 * Compare condition numbers for right invariant subspace
792 * taking its condition number into account
793 *
794 IF( V.GT.RCONDV*RCONDE ) THEN
795 TOL = RCONDV
796 ELSE
797 TOL = V / RCONDE
798 END IF
799 IF( V.GT.RCDVIN*RCDEIN ) THEN
800 TOLIN = RCDVIN
801 ELSE
802 TOLIN = V / RCDEIN
803 END IF
804 TOL = MAX( TOL, SMLNUM / EPS )
805 TOLIN = MAX( TOLIN, SMLNUM / EPS )
806 IF( EPS*( RCDVIN-TOLIN ).GT.RCONDV+TOL ) THEN
807 RESULT( 17 ) = ULPINV
808 ELSE IF( RCDVIN-TOLIN.GT.RCONDV+TOL ) THEN
809 RESULT( 17 ) = ( RCDVIN-TOLIN ) / ( RCONDV+TOL )
810 ELSE IF( RCDVIN+TOLIN.LT.EPS*( RCONDV-TOL ) ) THEN
811 RESULT( 17 ) = ULPINV
812 ELSE IF( RCDVIN+TOLIN.LT.RCONDV-TOL ) THEN
813 RESULT( 17 ) = ( RCONDV-TOL ) / ( RCDVIN+TOLIN )
814 ELSE
815 RESULT( 17 ) = ONE
816 END IF
817 *
818 270 CONTINUE
819 *
820 END IF
821 *
822 9999 FORMAT( ' ZGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
823 $ I6, ', INPUT EXAMPLE NUMBER = ', I4 )
824 9998 FORMAT( ' ZGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
825 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
826 *
827 RETURN
828 *
829 * End of ZGET24
830 *
831 END