3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
11 * SUBROUTINE ZCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
13 * S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1,
14 * ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK,
15 * RWORK, LLWORK, RESULT, INFO )
17 * .. Scalar Arguments ..
19 * INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
20 * DOUBLE PRECISION THRESH, THRSHN
22 * .. Array Arguments ..
23 * LOGICAL DOTYPE( * ), LLWORK( * )
24 * INTEGER ISEED( 4 ), NN( * )
25 * DOUBLE PRECISION RESULT( 15 ), RWORK( * )
26 * COMPLEX*16 A( LDA, * ), ALPHA1( * ), ALPHA3( * ),
27 * $ B( LDA, * ), BETA1( * ), BETA3( * ),
28 * $ EVECTL( LDU, * ), EVECTR( LDU, * ),
29 * $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ),
30 * $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ),
31 * $ T( LDA, * ), U( LDU, * ), V( LDU, * ),
32 * $ WORK( * ), Z( LDU, * )
41 *> ZCHKGG checks the nonsymmetric generalized eigenvalue problem
44 *> ZGGHRD factors A and B as U H V and U T V , where means conjugate
45 *> transpose, H is hessenberg, T is triangular and U and V are unitary.
48 *> ZHGEQZ factors H and T as Q S Z and Q P Z , where P and S are upper
49 *> triangular and Q and Z are unitary. It also computes the generalized
50 *> eigenvalues (alpha(1),beta(1)),...,(alpha(n),beta(n)), where
51 *> alpha(j)=S(j,j) and beta(j)=P(j,j) -- thus, w(j) = alpha(j)/beta(j)
52 *> is a root of the generalized eigenvalue problem
54 *> det( A - w(j) B ) = 0
56 *> and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent
59 *> det( m(j) A - B ) = 0
61 *> ZTGEVC computes the matrix L of left eigenvectors and the matrix R
62 *> of right eigenvectors for the matrix pair ( S, P ). In the
63 *> description below, l and r are left and right eigenvectors
64 *> corresponding to the generalized eigenvalues (alpha,beta).
66 *> When ZCHKGG is called, a number of matrix "sizes" ("n's") and a
67 *> number of matrix "types" are specified. For each size ("n")
68 *> and each type of matrix, one matrix will be generated and used
69 *> to test the nonsymmetric eigenroutines. For each matrix, 13
70 *> tests will be performed. The first twelve "test ratios" should be
71 *> small -- O(1). They will be compared with the threshold THRESH:
74 *> (1) | A - U H V | / ( |A| n ulp )
77 *> (2) | B - U T V | / ( |B| n ulp )
80 *> (3) | I - UU | / ( n ulp )
83 *> (4) | I - VV | / ( n ulp )
86 *> (5) | H - Q S Z | / ( |H| n ulp )
89 *> (6) | T - Q P Z | / ( |T| n ulp )
92 *> (7) | I - QQ | / ( n ulp )
95 *> (8) | I - ZZ | / ( n ulp )
97 *> (9) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
99 *> | (beta A - alpha B) l | / ( ulp max( |beta A|, |alpha B| ) )
101 *> (10) max over all left eigenvalue/-vector pairs (beta/alpha,l') of
103 *> | (beta H - alpha T) l' | / ( ulp max( |beta H|, |alpha T| ) )
105 *> where the eigenvectors l' are the result of passing Q to
106 *> DTGEVC and back transforming (JOB='B').
108 *> (11) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
110 *> | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
112 *> (12) max over all right eigenvalue/-vector pairs (beta/alpha,r') of
114 *> | (beta H - alpha T) r' | / ( ulp max( |beta H|, |alpha T| ) )
116 *> where the eigenvectors r' are the result of passing Z to
117 *> DTGEVC and back transforming (JOB='B').
119 *> The last three test ratios will usually be small, but there is no
120 *> mathematical requirement that they be so. They are therefore
121 *> compared with THRESH only if TSTDIF is .TRUE.
123 *> (13) | S(Q,Z computed) - S(Q,Z not computed) | / ( |S| ulp )
125 *> (14) | P(Q,Z computed) - P(Q,Z not computed) | / ( |P| ulp )
127 *> (15) max( |alpha(Q,Z computed) - alpha(Q,Z not computed)|/|S| ,
128 *> |beta(Q,Z computed) - beta(Q,Z not computed)|/|P| ) / ulp
130 *> In addition, the normalization of L and R are checked, and compared
131 *> with the threshold THRSHN.
136 *> The sizes of the test matrices are specified by an array
137 *> NN(1:NSIZES); the value of each element NN(j) specifies one size.
138 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
139 *> DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
140 *> Currently, the list of possible types is:
142 *> (1) ( 0, 0 ) (a pair of zero matrices)
144 *> (2) ( I, 0 ) (an identity and a zero matrix)
146 *> (3) ( 0, I ) (an identity and a zero matrix)
148 *> (4) ( I, I ) (a pair of identity matrices)
151 *> (5) ( J , J ) (a pair of transposed Jordan blocks)
154 *> (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
156 *> and I is a k x k identity and J a (k+1)x(k+1)
157 *> Jordan block; k=(N-1)/2
159 *> (7) ( D, I ) where D is P*D1, P is a random unitary diagonal
160 *> matrix (i.e., with random magnitude 1 entries
161 *> on the diagonal), and D1=diag( 0, 1,..., N-1 )
162 *> (i.e., a diagonal matrix with D1(1,1)=0,
163 *> D1(2,2)=1, ..., D1(N,N)=N-1.)
166 *> (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
168 *> (10) ( small*D, big*I )
170 *> (11) ( big*I, small*D )
172 *> (12) ( small*I, big*D )
174 *> (13) ( big*D, big*I )
176 *> (14) ( small*D, small*I )
178 *> (15) ( D1, D2 ) where D1=P*diag( 0, 0, 1, ..., N-3, 0 ) and
179 *> D2=Q*diag( 0, N-3, N-4,..., 1, 0, 0 ), and
180 *> P and Q are random unitary diagonal matrices.
182 *> (16) U ( J , J ) V where U and V are random unitary matrices.
184 *> (17) U ( T1, T2 ) V where T1 and T2 are upper triangular matrices
185 *> with random O(1) entries above the diagonal
186 *> and diagonal entries diag(T1) =
187 *> P*( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
188 *> Q*( 0, N-3, N-4,..., 1, 0, 0 )
190 *> (18) U ( T1, T2 ) V diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
191 *> diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
192 *> s = machine precision.
194 *> (19) U ( T1, T2 ) V diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
195 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
198 *> (20) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
199 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
201 *> (21) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
202 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
203 *> where r1,..., r(N-4) are random.
205 *> (22) U ( big*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
206 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
208 *> (23) U ( small*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
209 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
211 *> (24) U ( small*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
212 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
214 *> (25) U ( big*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
215 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
217 *> (26) U ( T1, T2 ) V where T1 and T2 are random upper-triangular
227 *> The number of sizes of matrices to use. If it is zero,
228 *> ZCHKGG does nothing. It must be at least zero.
233 *> NN is INTEGER array, dimension (NSIZES)
234 *> An array containing the sizes to be used for the matrices.
235 *> Zero values will be skipped. The values must be at least
242 *> The number of elements in DOTYPE. If it is zero, ZCHKGG
243 *> does nothing. It must be at least zero. If it is MAXTYP+1
244 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
245 *> defined, which is to use whatever matrix is in A. This
246 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
247 *> DOTYPE(MAXTYP+1) is .TRUE. .
252 *> DOTYPE is LOGICAL array, dimension (NTYPES)
253 *> If DOTYPE(j) is .TRUE., then for each size in NN a
254 *> matrix of that size and of type j will be generated.
255 *> If NTYPES is smaller than the maximum number of types
256 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
257 *> MAXTYP will not be generated. If NTYPES is larger
258 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
262 *> \param[in,out] ISEED
264 *> ISEED is INTEGER array, dimension (4)
265 *> On entry ISEED specifies the seed of the random number
266 *> generator. The array elements should be between 0 and 4095;
267 *> if not they will be reduced mod 4096. Also, ISEED(4) must
268 *> be odd. The random number generator uses a linear
269 *> congruential sequence limited to small integers, and so
270 *> should produce machine independent random numbers. The
271 *> values of ISEED are changed on exit, and can be used in the
272 *> next call to ZCHKGG to continue the same random number
278 *> THRESH is DOUBLE PRECISION
279 *> A test will count as "failed" if the "error", computed as
280 *> described above, exceeds THRESH. Note that the error
281 *> is scaled to be O(1), so THRESH should be a reasonably
282 *> small multiple of 1, e.g., 10 or 100. In particular,
283 *> it should not depend on the precision (single vs. double)
284 *> or the size of the matrix. It must be at least zero.
290 *> Specifies whether test ratios 13-15 will be computed and
291 *> compared with THRESH.
292 *> = .FALSE.: Only test ratios 1-12 will be computed and tested.
293 *> Ratios 13-15 will be set to zero.
294 *> = .TRUE.: All the test ratios 1-15 will be computed and
300 *> THRSHN is DOUBLE PRECISION
301 *> Threshold for reporting eigenvector normalization error.
302 *> If the normalization of any eigenvector differs from 1 by
303 *> more than THRSHN*ulp, then a special error message will be
304 *> printed. (This is handled separately from the other tests,
305 *> since only a compiler or programming error should cause an
306 *> error message, at least if THRSHN is at least 5--10.)
312 *> The FORTRAN unit number for printing out error messages
313 *> (e.g., if a routine returns IINFO not equal to 0.)
318 *> A is COMPLEX*16 array, dimension (LDA, max(NN))
319 *> Used to hold the original A matrix. Used as input only
320 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
321 *> DOTYPE(MAXTYP+1)=.TRUE.
327 *> The leading dimension of A, B, H, T, S1, P1, S2, and P2.
328 *> It must be at least 1 and at least max( NN ).
333 *> B is COMPLEX*16 array, dimension (LDA, max(NN))
334 *> Used to hold the original B matrix. Used as input only
335 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
336 *> DOTYPE(MAXTYP+1)=.TRUE.
341 *> H is COMPLEX*16 array, dimension (LDA, max(NN))
342 *> The upper Hessenberg matrix computed from A by ZGGHRD.
347 *> T is COMPLEX*16 array, dimension (LDA, max(NN))
348 *> The upper triangular matrix computed from B by ZGGHRD.
353 *> S1 is COMPLEX*16 array, dimension (LDA, max(NN))
354 *> The Schur (upper triangular) matrix computed from H by ZHGEQZ
355 *> when Q and Z are also computed.
360 *> S2 is COMPLEX*16 array, dimension (LDA, max(NN))
361 *> The Schur (upper triangular) matrix computed from H by ZHGEQZ
362 *> when Q and Z are not computed.
367 *> P1 is COMPLEX*16 array, dimension (LDA, max(NN))
368 *> The upper triangular matrix computed from T by ZHGEQZ
369 *> when Q and Z are also computed.
374 *> P2 is COMPLEX*16 array, dimension (LDA, max(NN))
375 *> The upper triangular matrix computed from T by ZHGEQZ
376 *> when Q and Z are not computed.
381 *> U is COMPLEX*16 array, dimension (LDU, max(NN))
382 *> The (left) unitary matrix computed by ZGGHRD.
388 *> The leading dimension of U, V, Q, Z, EVECTL, and EVEZTR. It
389 *> must be at least 1 and at least max( NN ).
394 *> V is COMPLEX*16 array, dimension (LDU, max(NN))
395 *> The (right) unitary matrix computed by ZGGHRD.
400 *> Q is COMPLEX*16 array, dimension (LDU, max(NN))
401 *> The (left) unitary matrix computed by ZHGEQZ.
406 *> Z is COMPLEX*16 array, dimension (LDU, max(NN))
407 *> The (left) unitary matrix computed by ZHGEQZ.
410 *> \param[out] ALPHA1
412 *> ALPHA1 is COMPLEX*16 array, dimension (max(NN))
417 *> BETA1 is COMPLEX*16 array, dimension (max(NN))
418 *> The generalized eigenvalues of (A,B) computed by ZHGEQZ
419 *> when Q, Z, and the full Schur matrices are computed.
422 *> \param[out] ALPHA3
424 *> ALPHA3 is COMPLEX*16 array, dimension (max(NN))
429 *> BETA3 is COMPLEX*16 array, dimension (max(NN))
430 *> The generalized eigenvalues of (A,B) computed by ZHGEQZ
431 *> when neither Q, Z, nor the Schur matrices are computed.
434 *> \param[out] EVECTL
436 *> EVECTL is COMPLEX*16 array, dimension (LDU, max(NN))
437 *> The (lower triangular) left eigenvector matrix for the
438 *> matrices in S1 and P1.
441 *> \param[out] EVECTR
443 *> EVECTR is COMPLEX*16 array, dimension (LDU, max(NN))
444 *> The (upper triangular) right eigenvector matrix for the
445 *> matrices in S1 and P1.
450 *> WORK is COMPLEX*16 array, dimension (LWORK)
456 *> The number of entries in WORK. This must be at least
457 *> max( 4*N, 2 * N**2, 1 ), for all N=NN(j).
462 *> RWORK is DOUBLE PRECISION array, dimension (2*max(NN))
465 *> \param[out] LLWORK
467 *> LLWORK is LOGICAL array, dimension (max(NN))
470 *> \param[out] RESULT
472 *> RESULT is DOUBLE PRECISION array, dimension (15)
473 *> The values computed by the tests described above.
474 *> The values are currently limited to 1/ulp, to avoid
481 *> = 0: successful exit.
482 *> < 0: if INFO = -i, the i-th argument had an illegal value.
483 *> > 0: A routine returned an error code. INFO is the
484 *> absolute value of the INFO value returned.
490 *> \author Univ. of Tennessee
491 *> \author Univ. of California Berkeley
492 *> \author Univ. of Colorado Denver
497 *> \ingroup complex16_eig
499 * =====================================================================
500 SUBROUTINE ZCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
501 $ TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
502 $ S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1,
503 $ ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK,
504 $ RWORK, LLWORK, RESULT, INFO )
506 * -- LAPACK test routine (version 3.6.1) --
507 * -- LAPACK is a software package provided by Univ. of Tennessee, --
508 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
511 * .. Scalar Arguments ..
513 INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
514 DOUBLE PRECISION THRESH, THRSHN
516 * .. Array Arguments ..
517 LOGICAL DOTYPE( * ), LLWORK( * )
518 INTEGER ISEED( 4 ), NN( * )
519 DOUBLE PRECISION RESULT( 15 ), RWORK( * )
520 COMPLEX*16 A( LDA, * ), ALPHA1( * ), ALPHA3( * ),
521 $ B( LDA, * ), BETA1( * ), BETA3( * ),
522 $ EVECTL( LDU, * ), EVECTR( LDU, * ),
523 $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ),
524 $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ),
525 $ T( LDA, * ), U( LDU, * ), V( LDU, * ),
526 $ WORK( * ), Z( LDU, * )
529 * =====================================================================
532 DOUBLE PRECISION ZERO, ONE
533 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
534 COMPLEX*16 CZERO, CONE
535 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
536 $ CONE = ( 1.0D+0, 0.0D+0 ) )
538 PARAMETER ( MAXTYP = 26 )
540 * .. Local Scalars ..
542 INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE,
543 $ LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX,
545 DOUBLE PRECISION ANORM, BNORM, SAFMAX, SAFMIN, TEMP1, TEMP2,
550 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
551 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
552 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
553 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
554 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ),
555 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
556 DOUBLE PRECISION DUMMA( 4 ), RMAGN( 0: 3 )
557 COMPLEX*16 CDUMMA( 4 )
559 * .. External Functions ..
560 DOUBLE PRECISION DLAMCH, ZLANGE
562 EXTERNAL DLAMCH, ZLANGE, ZLARND
564 * .. External Subroutines ..
565 EXTERNAL DLABAD, DLASUM, XERBLA, ZGEQR2, ZGET51, ZGET52,
566 $ ZGGHRD, ZHGEQZ, ZLACPY, ZLARFG, ZLASET, ZLATM4,
569 * .. Intrinsic Functions ..
570 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, SIGN
572 * .. Data statements ..
573 DATA KCLASS / 15*1, 10*2, 1*3 /
574 DATA KZ1 / 0, 1, 2, 1, 3, 3 /
575 DATA KZ2 / 0, 0, 1, 2, 1, 1 /
576 DATA KADD / 0, 0, 0, 0, 3, 2 /
577 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
578 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
579 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
580 $ 1, 1, -4, 2, -4, 8*8, 0 /
581 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
583 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
585 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
587 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
589 DATA KTRIAN / 16*0, 10*1 /
590 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE.,
591 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE.,
592 $ 3*.FALSE., 5*.TRUE., .FALSE. /
593 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE.,
594 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE.,
597 * .. Executable Statements ..
606 NMAX = MAX( NMAX, NN( J ) )
611 LWKOPT = MAX( 2*NMAX*NMAX, 4*NMAX, 1 )
615 IF( NSIZES.LT.0 ) THEN
617 ELSE IF( BADNN ) THEN
619 ELSE IF( NTYPES.LT.0 ) THEN
621 ELSE IF( THRESH.LT.ZERO ) THEN
623 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
625 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN
627 ELSE IF( LWKOPT.GT.LWORK ) THEN
632 CALL XERBLA( 'ZCHKGG', -INFO )
636 * Quick return if possible
638 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
641 SAFMIN = DLAMCH( 'Safe minimum' )
642 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
643 SAFMIN = SAFMIN / ULP
644 SAFMAX = ONE / SAFMIN
645 CALL DLABAD( SAFMIN, SAFMAX )
648 * The values RMAGN(2:3) depend on N, see below.
653 * Loop over sizes, types
659 DO 240 JSIZE = 1, NSIZES
662 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 )
663 RMAGN( 3 ) = SAFMIN*ULPINV*N1
665 IF( NSIZES.NE.1 ) THEN
666 MTYPES = MIN( MAXTYP, NTYPES )
668 MTYPES = MIN( MAXTYP+1, NTYPES )
671 DO 230 JTYPE = 1, MTYPES
672 IF( .NOT.DOTYPE( JTYPE ) )
677 * Save ISEED in case of an error.
680 IOLDSD( J ) = ISEED( J )
691 * Description of control parameters:
693 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
695 * KATYPE: the "type" to be passed to ZLATM4 for computing A.
696 * KAZERO: the pattern of zeros on the diagonal for A:
697 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
698 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
699 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
701 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
702 * =2: large, =3: small.
703 * LASIGN: .TRUE. if the diagonal elements of A are to be
704 * multiplied by a random magnitude 1 number.
705 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
706 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
707 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
708 * RMAGN: used to implement KAMAGN and KBMAGN.
710 IF( MTYPES.GT.MAXTYP )
713 IF( KCLASS( JTYPE ).LT.3 ) THEN
715 * Generate A (w/o rotation)
717 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
718 IN = 2*( ( N-1 ) / 2 ) + 1
720 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
724 CALL ZLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
725 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ),
726 $ RMAGN( KAMAGN( JTYPE ) ), ULP,
727 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 4,
729 IADD = KADD( KAZERO( JTYPE ) )
730 IF( IADD.GT.0 .AND. IADD.LE.N )
731 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) )
733 * Generate B (w/o rotation)
735 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
736 IN = 2*( ( N-1 ) / 2 ) + 1
738 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, B, LDA )
742 CALL ZLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
743 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ),
744 $ RMAGN( KBMAGN( JTYPE ) ), ONE,
745 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 4,
747 IADD = KADD( KBZERO( JTYPE ) )
749 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) )
751 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
755 * Generate U, V as Householder transformations times a
756 * diagonal matrix. (Note that ZLARFG makes U(j,j) and
761 U( JR, JC ) = ZLARND( 3, ISEED )
762 V( JR, JC ) = ZLARND( 3, ISEED )
764 CALL ZLARFG( N+1-JC, U( JC, JC ), U( JC+1, JC ), 1,
766 WORK( 2*N+JC ) = SIGN( ONE, DBLE( U( JC, JC ) ) )
768 CALL ZLARFG( N+1-JC, V( JC, JC ), V( JC+1, JC ), 1,
770 WORK( 3*N+JC ) = SIGN( ONE, DBLE( V( JC, JC ) ) )
773 CTEMP = ZLARND( 3, ISEED )
776 WORK( 3*N ) = CTEMP / ABS( CTEMP )
777 CTEMP = ZLARND( 3, ISEED )
780 WORK( 4*N ) = CTEMP / ABS( CTEMP )
782 * Apply the diagonal matrices
786 A( JR, JC ) = WORK( 2*N+JR )*
787 $ DCONJG( WORK( 3*N+JC ) )*
789 B( JR, JC ) = WORK( 2*N+JR )*
790 $ DCONJG( WORK( 3*N+JC ) )*
794 CALL ZUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, A,
795 $ LDA, WORK( 2*N+1 ), IINFO )
798 CALL ZUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ),
799 $ A, LDA, WORK( 2*N+1 ), IINFO )
802 CALL ZUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, B,
803 $ LDA, WORK( 2*N+1 ), IINFO )
806 CALL ZUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ),
807 $ B, LDA, WORK( 2*N+1 ), IINFO )
817 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
819 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
825 ANORM = ZLANGE( '1', N, N, A, LDA, RWORK )
826 BNORM = ZLANGE( '1', N, N, B, LDA, RWORK )
830 IF( IINFO.NE.0 ) THEN
831 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
839 * Call ZGEQR2, ZUNM2R, and ZGGHRD to compute H, T, U, and V
841 CALL ZLACPY( ' ', N, N, A, LDA, H, LDA )
842 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
846 CALL ZGEQR2( N, N, T, LDA, WORK, WORK( N+1 ), IINFO )
847 IF( IINFO.NE.0 ) THEN
848 WRITE( NOUNIT, FMT = 9999 )'ZGEQR2', IINFO, N, JTYPE,
854 CALL ZUNM2R( 'L', 'C', N, N, N, T, LDA, WORK, H, LDA,
855 $ WORK( N+1 ), IINFO )
856 IF( IINFO.NE.0 ) THEN
857 WRITE( NOUNIT, FMT = 9999 )'ZUNM2R', IINFO, N, JTYPE,
863 CALL ZLASET( 'Full', N, N, CZERO, CONE, U, LDU )
864 CALL ZUNM2R( 'R', 'N', N, N, N, T, LDA, WORK, U, LDU,
865 $ WORK( N+1 ), IINFO )
866 IF( IINFO.NE.0 ) THEN
867 WRITE( NOUNIT, FMT = 9999 )'ZUNM2R', IINFO, N, JTYPE,
873 CALL ZGGHRD( 'V', 'I', N, 1, N, H, LDA, T, LDA, U, LDU, V,
875 IF( IINFO.NE.0 ) THEN
876 WRITE( NOUNIT, FMT = 9999 )'ZGGHRD', IINFO, N, JTYPE,
885 CALL ZGET51( 1, N, A, LDA, H, LDA, U, LDU, V, LDU, WORK,
886 $ RWORK, RESULT( 1 ) )
887 CALL ZGET51( 1, N, B, LDA, T, LDA, U, LDU, V, LDU, WORK,
888 $ RWORK, RESULT( 2 ) )
889 CALL ZGET51( 3, N, B, LDA, T, LDA, U, LDU, U, LDU, WORK,
890 $ RWORK, RESULT( 3 ) )
891 CALL ZGET51( 3, N, B, LDA, T, LDA, V, LDU, V, LDU, WORK,
892 $ RWORK, RESULT( 4 ) )
894 * Call ZHGEQZ to compute S1, P1, S2, P2, Q, and Z, do tests.
900 CALL ZLACPY( ' ', N, N, H, LDA, S2, LDA )
901 CALL ZLACPY( ' ', N, N, T, LDA, P2, LDA )
905 CALL ZHGEQZ( 'E', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
906 $ ALPHA3, BETA3, Q, LDU, Z, LDU, WORK, LWORK,
908 IF( IINFO.NE.0 ) THEN
909 WRITE( NOUNIT, FMT = 9999 )'ZHGEQZ(E)', IINFO, N, JTYPE,
915 * Eigenvalues and Full Schur Form
917 CALL ZLACPY( ' ', N, N, H, LDA, S2, LDA )
918 CALL ZLACPY( ' ', N, N, T, LDA, P2, LDA )
920 CALL ZHGEQZ( 'S', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
921 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
923 IF( IINFO.NE.0 ) THEN
924 WRITE( NOUNIT, FMT = 9999 )'ZHGEQZ(S)', IINFO, N, JTYPE,
930 * Eigenvalues, Schur Form, and Schur Vectors
932 CALL ZLACPY( ' ', N, N, H, LDA, S1, LDA )
933 CALL ZLACPY( ' ', N, N, T, LDA, P1, LDA )
935 CALL ZHGEQZ( 'S', 'I', 'I', N, 1, N, S1, LDA, P1, LDA,
936 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
938 IF( IINFO.NE.0 ) THEN
939 WRITE( NOUNIT, FMT = 9999 )'ZHGEQZ(V)', IINFO, N, JTYPE,
949 CALL ZGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK,
950 $ RWORK, RESULT( 5 ) )
951 CALL ZGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK,
952 $ RWORK, RESULT( 6 ) )
953 CALL ZGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK,
954 $ RWORK, RESULT( 7 ) )
955 CALL ZGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK,
956 $ RWORK, RESULT( 8 ) )
958 * Compute the Left and Right Eigenvectors of (S1,P1)
960 * 9: Compute the left eigenvector Matrix without
966 * To test "SELECT" option, compute half of the eigenvectors
967 * in one call, and half in another
974 LLWORK( J ) = .FALSE.
977 CALL ZTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
978 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
979 IF( IINFO.NE.0 ) THEN
980 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(L,S1)', IINFO, N,
988 LLWORK( J ) = .FALSE.
994 CALL ZTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA,
995 $ EVECTL( 1, I1+1 ), LDU, CDUMMA, LDU, N, IN,
996 $ WORK, RWORK, IINFO )
997 IF( IINFO.NE.0 ) THEN
998 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(L,S2)', IINFO, N,
1004 CALL ZGET52( .TRUE., N, S1, LDA, P1, LDA, EVECTL, LDU,
1005 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
1006 RESULT( 9 ) = DUMMA( 1 )
1007 IF( DUMMA( 2 ).GT.THRSHN ) THEN
1008 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTGEVC(HOWMNY=S)',
1009 $ DUMMA( 2 ), N, JTYPE, IOLDSD
1012 * 10: Compute the left eigenvector Matrix with
1013 * back transforming:
1016 RESULT( 10 ) = ULPINV
1017 CALL ZLACPY( 'F', N, N, Q, LDU, EVECTL, LDU )
1018 CALL ZTGEVC( 'L', 'B', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
1019 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
1020 IF( IINFO.NE.0 ) THEN
1021 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(L,B)', IINFO, N,
1027 CALL ZGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHA1,
1028 $ BETA1, WORK, RWORK, DUMMA( 1 ) )
1029 RESULT( 10 ) = DUMMA( 1 )
1030 IF( DUMMA( 2 ).GT.THRSHN ) THEN
1031 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTGEVC(HOWMNY=B)',
1032 $ DUMMA( 2 ), N, JTYPE, IOLDSD
1035 * 11: Compute the right eigenvector Matrix without
1036 * back transforming:
1039 RESULT( 11 ) = ULPINV
1041 * To test "SELECT" option, compute half of the eigenvectors
1042 * in one call, and half in another
1046 LLWORK( J ) = .TRUE.
1048 DO 170 J = I1 + 1, N
1049 LLWORK( J ) = .FALSE.
1052 CALL ZTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
1053 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
1054 IF( IINFO.NE.0 ) THEN
1055 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(R,S1)', IINFO, N,
1063 LLWORK( J ) = .FALSE.
1065 DO 190 J = I1 + 1, N
1066 LLWORK( J ) = .TRUE.
1069 CALL ZTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
1070 $ LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK,
1072 IF( IINFO.NE.0 ) THEN
1073 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(R,S2)', IINFO, N,
1079 CALL ZGET52( .FALSE., N, S1, LDA, P1, LDA, EVECTR, LDU,
1080 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
1081 RESULT( 11 ) = DUMMA( 1 )
1082 IF( DUMMA( 2 ).GT.THRESH ) THEN
1083 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTGEVC(HOWMNY=S)',
1084 $ DUMMA( 2 ), N, JTYPE, IOLDSD
1087 * 12: Compute the right eigenvector Matrix with
1088 * back transforming:
1091 RESULT( 12 ) = ULPINV
1092 CALL ZLACPY( 'F', N, N, Z, LDU, EVECTR, LDU )
1093 CALL ZTGEVC( 'R', 'B', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
1094 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
1095 IF( IINFO.NE.0 ) THEN
1096 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(R,B)', IINFO, N,
1102 CALL ZGET52( .FALSE., N, H, LDA, T, LDA, EVECTR, LDU,
1103 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
1104 RESULT( 12 ) = DUMMA( 1 )
1105 IF( DUMMA( 2 ).GT.THRESH ) THEN
1106 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTGEVC(HOWMNY=B)',
1107 $ DUMMA( 2 ), N, JTYPE, IOLDSD
1110 * Tests 13--15 are done only on request
1116 CALL ZGET51( 2, N, S1, LDA, S2, LDA, Q, LDU, Z, LDU,
1117 $ WORK, RWORK, RESULT( 13 ) )
1118 CALL ZGET51( 2, N, P1, LDA, P2, LDA, Q, LDU, Z, LDU,
1119 $ WORK, RWORK, RESULT( 14 ) )
1126 TEMP1 = MAX( TEMP1, ABS( ALPHA1( J )-ALPHA3( J ) ) )
1127 TEMP2 = MAX( TEMP2, ABS( BETA1( J )-BETA3( J ) ) )
1130 TEMP1 = TEMP1 / MAX( SAFMIN, ULP*MAX( TEMP1, ANORM ) )
1131 TEMP2 = TEMP2 / MAX( SAFMIN, ULP*MAX( TEMP2, BNORM ) )
1132 RESULT( 15 ) = MAX( TEMP1, TEMP2 )
1141 * End of Loop -- Check for RESULT(j) > THRESH
1145 NTESTT = NTESTT + NTEST
1147 * Print out tests which fail.
1149 DO 220 JR = 1, NTEST
1150 IF( RESULT( JR ).GE.THRESH ) THEN
1152 * If this is the first test to fail,
1153 * print a header to the data file.
1155 IF( NERRS.EQ.0 ) THEN
1156 WRITE( NOUNIT, FMT = 9997 )'ZGG'
1160 WRITE( NOUNIT, FMT = 9996 )
1161 WRITE( NOUNIT, FMT = 9995 )
1162 WRITE( NOUNIT, FMT = 9994 )'Unitary'
1166 WRITE( NOUNIT, FMT = 9993 )'unitary', '*',
1167 $ 'conjugate transpose', ( '*', J = 1, 10 )
1171 IF( RESULT( JR ).LT.10000.0D0 ) THEN
1172 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
1175 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
1186 CALL DLASUM( 'ZGG', NOUNIT, NERRS, NTESTT )
1189 9999 FORMAT( ' ZCHKGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
1190 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
1192 9998 FORMAT( ' ZCHKGG: ', A, ' Eigenvectors from ', A, ' incorrectly ',
1193 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
1194 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
1197 9997 FORMAT( 1X, A3, ' -- Complex Generalized eigenvalue problem' )
1199 9996 FORMAT( ' Matrix types (see ZCHKGG for details): ' )
1201 9995 FORMAT( ' Special Matrices:', 23X,
1202 $ '(J''=transposed Jordan block)',
1203 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
1204 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
1205 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
1206 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
1207 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
1208 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
1209 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
1210 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
1211 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
1212 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
1213 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
1214 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
1215 $ '23=(small,large) 24=(small,small) 25=(large,large)',
1216 $ / ' 26=random O(1) matrices.' )
1218 9993 FORMAT( / ' Tests performed: (H is Hessenberg, S is Schur, B, ',
1219 $ 'T, P are triangular,', / 20X, 'U, V, Q, and Z are ', A,
1220 $ ', l and r are the', / 20X,
1221 $ 'appropriate left and right eigenvectors, resp., a is',
1222 $ / 20X, 'alpha, b is beta, and ', A, ' means ', A, '.)',
1223 $ / ' 1 = | A - U H V', A,
1224 $ ' | / ( |A| n ulp ) 2 = | B - U T V', A,
1225 $ ' | / ( |B| n ulp )', / ' 3 = | I - UU', A,
1226 $ ' | / ( n ulp ) 4 = | I - VV', A,
1227 $ ' | / ( n ulp )', / ' 5 = | H - Q S Z', A,
1228 $ ' | / ( |H| n ulp )', 6X, '6 = | T - Q P Z', A,
1229 $ ' | / ( |T| n ulp )', / ' 7 = | I - QQ', A,
1230 $ ' | / ( n ulp ) 8 = | I - ZZ', A,
1231 $ ' | / ( n ulp )', / ' 9 = max | ( b S - a P )', A,
1232 $ ' l | / const. 10 = max | ( b H - a T )', A,
1233 $ ' l | / const.', /
1234 $ ' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H',
1235 $ ' - a T ) r | / const.', / 1X )
1237 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
1238 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
1239 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
1240 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 )