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