LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
cdrgev.f
Go to the documentation of this file.
1 *> \brief \b CDRGEV
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 CDRGEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE,
13 * ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK, RWORK,
14 * RESULT, INFO )
15 *
16 * .. Scalar Arguments ..
17 * INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
18 * $ NTYPES
19 * REAL THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL DOTYPE( * )
23 * INTEGER ISEED( 4 ), NN( * )
24 * REAL RESULT( * ), RWORK( * )
25 * COMPLEX A( LDA, * ), ALPHA( * ), ALPHA1( * ),
26 * $ B( LDA, * ), BETA( * ), BETA1( * ),
27 * $ Q( LDQ, * ), QE( LDQE, * ), S( LDA, * ),
28 * $ T( LDA, * ), WORK( * ), Z( LDQ, * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> CDRGEV checks the nonsymmetric generalized eigenvalue problem driver
38 *> routine CGGEV.
39 *>
40 *> CGGEV computes for a pair of n-by-n nonsymmetric matrices (A,B) the
41 *> generalized eigenvalues and, optionally, the left and right
42 *> eigenvectors.
43 *>
44 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
45 *> or a ratio alpha/beta = w, such that A - w*B is singular. It is
46 *> usually represented as the pair (alpha,beta), as there is reasonable
47 *> interpretation for beta=0, and even for both being zero.
48 *>
49 *> A right generalized eigenvector corresponding to a generalized
50 *> eigenvalue w for a pair of matrices (A,B) is a vector r such that
51 *> (A - wB) * r = 0. A left generalized eigenvector is a vector l such
52 *> that l**H * (A - wB) = 0, where l**H is the conjugate-transpose of l.
53 *>
54 *> When CDRGEV is called, a number of matrix "sizes" ("n's") and a
55 *> number of matrix "types" are specified. For each size ("n")
56 *> and each type of matrix, a pair of matrices (A, B) will be generated
57 *> and used for testing. For each matrix pair, the following tests
58 *> will be performed and compared with the threshold THRESH.
59 *>
60 *> Results from CGGEV:
61 *>
62 *> (1) max over all left eigenvalue/-vector pairs (alpha/beta,l) of
63 *>
64 *> | VL**H * (beta A - alpha B) |/( ulp max(|beta A|, |alpha B|) )
65 *>
66 *> where VL**H is the conjugate-transpose of VL.
67 *>
68 *> (2) | |VL(i)| - 1 | / ulp and whether largest component real
69 *>
70 *> VL(i) denotes the i-th column of VL.
71 *>
72 *> (3) max over all left eigenvalue/-vector pairs (alpha/beta,r) of
73 *>
74 *> | (beta A - alpha B) * VR | / ( ulp max(|beta A|, |alpha B|) )
75 *>
76 *> (4) | |VR(i)| - 1 | / ulp and whether largest component real
77 *>
78 *> VR(i) denotes the i-th column of VR.
79 *>
80 *> (5) W(full) = W(partial)
81 *> W(full) denotes the eigenvalues computed when both l and r
82 *> are also computed, and W(partial) denotes the eigenvalues
83 *> computed when only W, only W and r, or only W and l are
84 *> computed.
85 *>
86 *> (6) VL(full) = VL(partial)
87 *> VL(full) denotes the left eigenvectors computed when both l
88 *> and r are computed, and VL(partial) denotes the result
89 *> when only l is computed.
90 *>
91 *> (7) VR(full) = VR(partial)
92 *> VR(full) denotes the right eigenvectors computed when both l
93 *> and r are also computed, and VR(partial) denotes the result
94 *> when only l is computed.
95 *>
96 *>
97 *> Test Matrices
98 *> ---- --------
99 *>
100 *> The sizes of the test matrices are specified by an array
101 *> NN(1:NSIZES); the value of each element NN(j) specifies one size.
102 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
103 *> DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
104 *> Currently, the list of possible types is:
105 *>
106 *> (1) ( 0, 0 ) (a pair of zero matrices)
107 *>
108 *> (2) ( I, 0 ) (an identity and a zero matrix)
109 *>
110 *> (3) ( 0, I ) (an identity and a zero matrix)
111 *>
112 *> (4) ( I, I ) (a pair of identity matrices)
113 *>
114 *> t t
115 *> (5) ( J , J ) (a pair of transposed Jordan blocks)
116 *>
117 *> t ( I 0 )
118 *> (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
119 *> ( 0 I ) ( 0 J )
120 *> and I is a k x k identity and J a (k+1)x(k+1)
121 *> Jordan block; k=(N-1)/2
122 *>
123 *> (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal
124 *> matrix with those diagonal entries.)
125 *> (8) ( I, D )
126 *>
127 *> (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
128 *>
129 *> (10) ( small*D, big*I )
130 *>
131 *> (11) ( big*I, small*D )
132 *>
133 *> (12) ( small*I, big*D )
134 *>
135 *> (13) ( big*D, big*I )
136 *>
137 *> (14) ( small*D, small*I )
138 *>
139 *> (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
140 *> D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
141 *> t t
142 *> (16) Q ( J , J ) Z where Q and Z are random orthogonal matrices.
143 *>
144 *> (17) Q ( T1, T2 ) Z where T1 and T2 are upper triangular matrices
145 *> with random O(1) entries above the diagonal
146 *> and diagonal entries diag(T1) =
147 *> ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
148 *> ( 0, N-3, N-4,..., 1, 0, 0 )
149 *>
150 *> (18) Q ( T1, T2 ) Z diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
151 *> diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
152 *> s = machine precision.
153 *>
154 *> (19) Q ( T1, T2 ) Z diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
155 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
156 *>
157 *> N-5
158 *> (20) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
159 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
160 *>
161 *> (21) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
162 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
163 *> where r1,..., r(N-4) are random.
164 *>
165 *> (22) Q ( big*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
166 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
167 *>
168 *> (23) Q ( small*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
169 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
170 *>
171 *> (24) Q ( small*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
172 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
173 *>
174 *> (25) Q ( big*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
175 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
176 *>
177 *> (26) Q ( T1, T2 ) Z where T1 and T2 are random upper-triangular
178 *> matrices.
179 *>
180 *> \endverbatim
181 *
182 * Arguments:
183 * ==========
184 *
185 *> \param[in] NSIZES
186 *> \verbatim
187 *> NSIZES is INTEGER
188 *> The number of sizes of matrices to use. If it is zero,
189 *> CDRGES does nothing. NSIZES >= 0.
190 *> \endverbatim
191 *>
192 *> \param[in] NN
193 *> \verbatim
194 *> NN is INTEGER array, dimension (NSIZES)
195 *> An array containing the sizes to be used for the matrices.
196 *> Zero values will be skipped. NN >= 0.
197 *> \endverbatim
198 *>
199 *> \param[in] NTYPES
200 *> \verbatim
201 *> NTYPES is INTEGER
202 *> The number of elements in DOTYPE. If it is zero, CDRGEV
203 *> does nothing. It must be at least zero. If it is MAXTYP+1
204 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
205 *> defined, which is to use whatever matrix is in A. This
206 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
207 *> DOTYPE(MAXTYP+1) is .TRUE. .
208 *> \endverbatim
209 *>
210 *> \param[in] DOTYPE
211 *> \verbatim
212 *> DOTYPE is LOGICAL array, dimension (NTYPES)
213 *> If DOTYPE(j) is .TRUE., then for each size in NN a
214 *> matrix of that size and of type j will be generated.
215 *> If NTYPES is smaller than the maximum number of types
216 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
217 *> MAXTYP will not be generated. If NTYPES is larger
218 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
219 *> will be ignored.
220 *> \endverbatim
221 *>
222 *> \param[in,out] ISEED
223 *> \verbatim
224 *> ISEED is INTEGER array, dimension (4)
225 *> On entry ISEED specifies the seed of the random number
226 *> generator. The array elements should be between 0 and 4095;
227 *> if not they will be reduced mod 4096. Also, ISEED(4) must
228 *> be odd. The random number generator uses a linear
229 *> congruential sequence limited to small integers, and so
230 *> should produce machine independent random numbers. The
231 *> values of ISEED are changed on exit, and can be used in the
232 *> next call to CDRGES to continue the same random number
233 *> sequence.
234 *> \endverbatim
235 *>
236 *> \param[in] THRESH
237 *> \verbatim
238 *> THRESH is REAL
239 *> A test will count as "failed" if the "error", computed as
240 *> described above, exceeds THRESH. Note that the error is
241 *> scaled to be O(1), so THRESH should be a reasonably small
242 *> multiple of 1, e.g., 10 or 100. In particular, it should
243 *> not depend on the precision (single vs. double) or the size
244 *> of the matrix. It must be at least zero.
245 *> \endverbatim
246 *>
247 *> \param[in] NOUNIT
248 *> \verbatim
249 *> NOUNIT is INTEGER
250 *> The FORTRAN unit number for printing out error messages
251 *> (e.g., if a routine returns IERR not equal to 0.)
252 *> \endverbatim
253 *>
254 *> \param[in,out] A
255 *> \verbatim
256 *> A is COMPLEX array, dimension(LDA, max(NN))
257 *> Used to hold the original A matrix. Used as input only
258 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
259 *> DOTYPE(MAXTYP+1)=.TRUE.
260 *> \endverbatim
261 *>
262 *> \param[in] LDA
263 *> \verbatim
264 *> LDA is INTEGER
265 *> The leading dimension of A, B, S, and T.
266 *> It must be at least 1 and at least max( NN ).
267 *> \endverbatim
268 *>
269 *> \param[in,out] B
270 *> \verbatim
271 *> B is COMPLEX array, dimension(LDA, max(NN))
272 *> Used to hold the original B matrix. Used as input only
273 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
274 *> DOTYPE(MAXTYP+1)=.TRUE.
275 *> \endverbatim
276 *>
277 *> \param[out] S
278 *> \verbatim
279 *> S is COMPLEX array, dimension (LDA, max(NN))
280 *> The Schur form matrix computed from A by CGGEV. On exit, S
281 *> contains the Schur form matrix corresponding to the matrix
282 *> in A.
283 *> \endverbatim
284 *>
285 *> \param[out] T
286 *> \verbatim
287 *> T is COMPLEX array, dimension (LDA, max(NN))
288 *> The upper triangular matrix computed from B by CGGEV.
289 *> \endverbatim
290 *>
291 *> \param[out] Q
292 *> \verbatim
293 *> Q is COMPLEX array, dimension (LDQ, max(NN))
294 *> The (left) eigenvectors matrix computed by CGGEV.
295 *> \endverbatim
296 *>
297 *> \param[in] LDQ
298 *> \verbatim
299 *> LDQ is INTEGER
300 *> The leading dimension of Q and Z. It must
301 *> be at least 1 and at least max( NN ).
302 *> \endverbatim
303 *>
304 *> \param[out] Z
305 *> \verbatim
306 *> Z is COMPLEX array, dimension( LDQ, max(NN) )
307 *> The (right) orthogonal matrix computed by CGGEV.
308 *> \endverbatim
309 *>
310 *> \param[out] QE
311 *> \verbatim
312 *> QE is COMPLEX array, dimension( LDQ, max(NN) )
313 *> QE holds the computed right or left eigenvectors.
314 *> \endverbatim
315 *>
316 *> \param[in] LDQE
317 *> \verbatim
318 *> LDQE is INTEGER
319 *> The leading dimension of QE. LDQE >= max(1,max(NN)).
320 *> \endverbatim
321 *>
322 *> \param[out] ALPHA
323 *> \verbatim
324 *> ALPHA is COMPLEX array, dimension (max(NN))
325 *> \endverbatim
326 *>
327 *> \param[out] BETA
328 *> \verbatim
329 *> BETA is COMPLEX array, dimension (max(NN))
330 *>
331 *> The generalized eigenvalues of (A,B) computed by CGGEV.
332 *> ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
333 *> generalized eigenvalue of A and B.
334 *> \endverbatim
335 *>
336 *> \param[out] ALPHA1
337 *> \verbatim
338 *> ALPHA1 is COMPLEX array, dimension (max(NN))
339 *> \endverbatim
340 *>
341 *> \param[out] BETA1
342 *> \verbatim
343 *> BETA1 is COMPLEX array, dimension (max(NN))
344 *>
345 *> Like ALPHAR, ALPHAI, BETA, these arrays contain the
346 *> eigenvalues of A and B, but those computed when CGGEV only
347 *> computes a partial eigendecomposition, i.e. not the
348 *> eigenvalues and left and right eigenvectors.
349 *> \endverbatim
350 *>
351 *> \param[out] WORK
352 *> \verbatim
353 *> WORK is COMPLEX array, dimension (LWORK)
354 *> \endverbatim
355 *>
356 *> \param[in] LWORK
357 *> \verbatim
358 *> LWORK is INTEGER
359 *> The number of entries in WORK. LWORK >= N*(N+1)
360 *> \endverbatim
361 *>
362 *> \param[out] RWORK
363 *> \verbatim
364 *> RWORK is REAL array, dimension (8*N)
365 *> Real workspace.
366 *> \endverbatim
367 *>
368 *> \param[out] RESULT
369 *> \verbatim
370 *> RESULT is REAL array, dimension (2)
371 *> The values computed by the tests described above.
372 *> The values are currently limited to 1/ulp, to avoid overflow.
373 *> \endverbatim
374 *>
375 *> \param[out] INFO
376 *> \verbatim
377 *> INFO is INTEGER
378 *> = 0: successful exit
379 *> < 0: if INFO = -i, the i-th argument had an illegal value.
380 *> > 0: A routine returned an error code. INFO is the
381 *> absolute value of the INFO value returned.
382 *> \endverbatim
383 *
384 * Authors:
385 * ========
386 *
387 *> \author Univ. of Tennessee
388 *> \author Univ. of California Berkeley
389 *> \author Univ. of Colorado Denver
390 *> \author NAG Ltd.
391 *
392 *> \date June 2016
393 *
394 *> \ingroup complex_eig
395 *
396 * =====================================================================
397  SUBROUTINE cdrgev( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
398  $ nounit, a, lda, b, s, t, q, ldq, z, qe, ldqe,
399  $ alpha, beta, alpha1, beta1, work, lwork, rwork,
400  $ result, info )
401 *
402 * -- LAPACK test routine (version 3.6.1) --
403 * -- LAPACK is a software package provided by Univ. of Tennessee, --
404 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
405 * June 2016
406 *
407 * .. Scalar Arguments ..
408  INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
409  $ ntypes
410  REAL THRESH
411 * ..
412 * .. Array Arguments ..
413  LOGICAL DOTYPE( * )
414  INTEGER ISEED( 4 ), NN( * )
415  REAL RESULT( * ), RWORK( * )
416  COMPLEX A( lda, * ), ALPHA( * ), ALPHA1( * ),
417  $ b( lda, * ), beta( * ), beta1( * ),
418  $ q( ldq, * ), qe( ldqe, * ), s( lda, * ),
419  $ t( lda, * ), work( * ), z( ldq, * )
420 * ..
421 *
422 * =====================================================================
423 *
424 * .. Parameters ..
425  REAL ZERO, ONE
426  parameter ( zero = 0.0e+0, one = 1.0e+0 )
427  COMPLEX CZERO, CONE
428  parameter ( czero = ( 0.0e+0, 0.0e+0 ),
429  $ cone = ( 1.0e+0, 0.0e+0 ) )
430  INTEGER MAXTYP
431  parameter ( maxtyp = 26 )
432 * ..
433 * .. Local Scalars ..
434  LOGICAL BADNN
435  INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
436  $ maxwrk, minwrk, mtypes, n, n1, nb, nerrs,
437  $ nmats, nmax, ntestt
438  REAL SAFMAX, SAFMIN, ULP, ULPINV
439  COMPLEX CTEMP
440 * ..
441 * .. Local Arrays ..
442  LOGICAL LASIGN( maxtyp ), LBSIGN( maxtyp )
443  INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( maxtyp ),
444  $ katype( maxtyp ), kazero( maxtyp ),
445  $ kbmagn( maxtyp ), kbtype( maxtyp ),
446  $ kbzero( maxtyp ), kclass( maxtyp ),
447  $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
448  REAL RMAGN( 0: 3 )
449 * ..
450 * .. External Functions ..
451  INTEGER ILAENV
452  REAL SLAMCH
453  COMPLEX CLARND
454  EXTERNAL ilaenv, slamch, clarnd
455 * ..
456 * .. External Subroutines ..
457  EXTERNAL alasvm, cget52, cggev, clacpy, clarfg, claset,
459 * ..
460 * .. Intrinsic Functions ..
461  INTRINSIC abs, conjg, max, min, REAL, SIGN
462 * ..
463 * .. Data statements ..
464  DATA kclass / 15*1, 10*2, 1*3 /
465  DATA kz1 / 0, 1, 2, 1, 3, 3 /
466  DATA kz2 / 0, 0, 1, 2, 1, 1 /
467  DATA kadd / 0, 0, 0, 0, 3, 2 /
468  DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
469  $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
470  DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
471  $ 1, 1, -4, 2, -4, 8*8, 0 /
472  DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
473  $ 4*5, 4*3, 1 /
474  DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
475  $ 4*6, 4*4, 1 /
476  DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
477  $ 2, 1 /
478  DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
479  $ 2, 1 /
480  DATA ktrian / 16*0, 10*1 /
481  DATA lasign / 6*.false., .true., .false., 2*.true.,
482  $ 2*.false., 3*.true., .false., .true.,
483  $ 3*.false., 5*.true., .false. /
484  DATA lbsign / 7*.false., .true., 2*.false.,
485  $ 2*.true., 2*.false., .true., .false., .true.,
486  $ 9*.false. /
487 * ..
488 * .. Executable Statements ..
489 *
490 * Check for errors
491 *
492  info = 0
493 *
494  badnn = .false.
495  nmax = 1
496  DO 10 j = 1, nsizes
497  nmax = max( nmax, nn( j ) )
498  IF( nn( j ).LT.0 )
499  $ badnn = .true.
500  10 CONTINUE
501 *
502  IF( nsizes.LT.0 ) THEN
503  info = -1
504  ELSE IF( badnn ) THEN
505  info = -2
506  ELSE IF( ntypes.LT.0 ) THEN
507  info = -3
508  ELSE IF( thresh.LT.zero ) THEN
509  info = -6
510  ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
511  info = -9
512  ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
513  info = -14
514  ELSE IF( ldqe.LE.1 .OR. ldqe.LT.nmax ) THEN
515  info = -17
516  END IF
517 *
518 * Compute workspace
519 * (Note: Comments in the code beginning "Workspace:" describe the
520 * minimal amount of workspace needed at that point in the code,
521 * as well as the preferred amount for good performance.
522 * NB refers to the optimal block size for the immediately
523 * following subroutine, as returned by ILAENV.
524 *
525  minwrk = 1
526  IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
527  minwrk = nmax*( nmax+1 )
528  nb = max( 1, ilaenv( 1, 'CGEQRF', ' ', nmax, nmax, -1, -1 ),
529  $ ilaenv( 1, 'CUNMQR', 'LC', nmax, nmax, nmax, -1 ),
530  $ ilaenv( 1, 'CUNGQR', ' ', nmax, nmax, nmax, -1 ) )
531  maxwrk = max( 2*nmax, nmax*( nb+1 ), nmax*( nmax+1 ) )
532  work( 1 ) = maxwrk
533  END IF
534 *
535  IF( lwork.LT.minwrk )
536  $ info = -23
537 *
538  IF( info.NE.0 ) THEN
539  CALL xerbla( 'CDRGEV', -info )
540  RETURN
541  END IF
542 *
543 * Quick return if possible
544 *
545  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
546  $ RETURN
547 *
548  ulp = slamch( 'Precision' )
549  safmin = slamch( 'Safe minimum' )
550  safmin = safmin / ulp
551  safmax = one / safmin
552  CALL slabad( safmin, safmax )
553  ulpinv = one / ulp
554 *
555 * The values RMAGN(2:3) depend on N, see below.
556 *
557  rmagn( 0 ) = zero
558  rmagn( 1 ) = one
559 *
560 * Loop over sizes, types
561 *
562  ntestt = 0
563  nerrs = 0
564  nmats = 0
565 *
566  DO 220 jsize = 1, nsizes
567  n = nn( jsize )
568  n1 = max( 1, n )
569  rmagn( 2 ) = safmax*ulp / REAL( n1 )
570  rmagn( 3 ) = safmin*ulpinv*n1
571 *
572  IF( nsizes.NE.1 ) THEN
573  mtypes = min( maxtyp, ntypes )
574  ELSE
575  mtypes = min( maxtyp+1, ntypes )
576  END IF
577 *
578  DO 210 jtype = 1, mtypes
579  IF( .NOT.dotype( jtype ) )
580  $ GO TO 210
581  nmats = nmats + 1
582 *
583 * Save ISEED in case of an error.
584 *
585  DO 20 j = 1, 4
586  ioldsd( j ) = iseed( j )
587  20 CONTINUE
588 *
589 * Generate test matrices A and B
590 *
591 * Description of control parameters:
592 *
593 * KCLASS: =1 means w/o rotation, =2 means w/ rotation,
594 * =3 means random.
595 * KATYPE: the "type" to be passed to CLATM4 for computing A.
596 * KAZERO: the pattern of zeros on the diagonal for A:
597 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
598 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
599 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
600 * non-zero entries.)
601 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
602 * =2: large, =3: small.
603 * LASIGN: .TRUE. if the diagonal elements of A are to be
604 * multiplied by a random magnitude 1 number.
605 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
606 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
607 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
608 * RMAGN: used to implement KAMAGN and KBMAGN.
609 *
610  IF( mtypes.GT.maxtyp )
611  $ GO TO 100
612  ierr = 0
613  IF( kclass( jtype ).LT.3 ) THEN
614 *
615 * Generate A (w/o rotation)
616 *
617  IF( abs( katype( jtype ) ).EQ.3 ) THEN
618  in = 2*( ( n-1 ) / 2 ) + 1
619  IF( in.NE.n )
620  $ CALL claset( 'Full', n, n, czero, czero, a, lda )
621  ELSE
622  in = n
623  END IF
624  CALL clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
625  $ kz2( kazero( jtype ) ), lasign( jtype ),
626  $ rmagn( kamagn( jtype ) ), ulp,
627  $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
628  $ iseed, a, lda )
629  iadd = kadd( kazero( jtype ) )
630  IF( iadd.GT.0 .AND. iadd.LE.n )
631  $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
632 *
633 * Generate B (w/o rotation)
634 *
635  IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
636  in = 2*( ( n-1 ) / 2 ) + 1
637  IF( in.NE.n )
638  $ CALL claset( 'Full', n, n, czero, czero, b, lda )
639  ELSE
640  in = n
641  END IF
642  CALL clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
643  $ kz2( kbzero( jtype ) ), lbsign( jtype ),
644  $ rmagn( kbmagn( jtype ) ), one,
645  $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
646  $ iseed, b, lda )
647  iadd = kadd( kbzero( jtype ) )
648  IF( iadd.NE.0 .AND. iadd.LE.n )
649  $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
650 *
651  IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
652 *
653 * Include rotations
654 *
655 * Generate Q, Z as Householder transformations times
656 * a diagonal matrix.
657 *
658  DO 40 jc = 1, n - 1
659  DO 30 jr = jc, n
660  q( jr, jc ) = clarnd( 3, iseed )
661  z( jr, jc ) = clarnd( 3, iseed )
662  30 CONTINUE
663  CALL clarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
664  $ work( jc ) )
665  work( 2*n+jc ) = sign( one, REAL( Q( JC, JC ) ) )
666  q( jc, jc ) = cone
667  CALL clarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
668  $ work( n+jc ) )
669  work( 3*n+jc ) = sign( one, REAL( Z( JC, JC ) ) )
670  z( jc, jc ) = cone
671  40 CONTINUE
672  ctemp = clarnd( 3, iseed )
673  q( n, n ) = cone
674  work( n ) = czero
675  work( 3*n ) = ctemp / abs( ctemp )
676  ctemp = clarnd( 3, iseed )
677  z( n, n ) = cone
678  work( 2*n ) = czero
679  work( 4*n ) = ctemp / abs( ctemp )
680 *
681 * Apply the diagonal matrices
682 *
683  DO 60 jc = 1, n
684  DO 50 jr = 1, n
685  a( jr, jc ) = work( 2*n+jr )*
686  $ conjg( work( 3*n+jc ) )*
687  $ a( jr, jc )
688  b( jr, jc ) = work( 2*n+jr )*
689  $ conjg( work( 3*n+jc ) )*
690  $ b( jr, jc )
691  50 CONTINUE
692  60 CONTINUE
693  CALL cunm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
694  $ lda, work( 2*n+1 ), ierr )
695  IF( ierr.NE.0 )
696  $ GO TO 90
697  CALL cunm2r( 'R', 'C', n, n, n-1, z, ldq, work( n+1 ),
698  $ a, lda, work( 2*n+1 ), ierr )
699  IF( ierr.NE.0 )
700  $ GO TO 90
701  CALL cunm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
702  $ lda, work( 2*n+1 ), ierr )
703  IF( ierr.NE.0 )
704  $ GO TO 90
705  CALL cunm2r( 'R', 'C', n, n, n-1, z, ldq, work( n+1 ),
706  $ b, lda, work( 2*n+1 ), ierr )
707  IF( ierr.NE.0 )
708  $ GO TO 90
709  END IF
710  ELSE
711 *
712 * Random matrices
713 *
714  DO 80 jc = 1, n
715  DO 70 jr = 1, n
716  a( jr, jc ) = rmagn( kamagn( jtype ) )*
717  $ clarnd( 4, iseed )
718  b( jr, jc ) = rmagn( kbmagn( jtype ) )*
719  $ clarnd( 4, iseed )
720  70 CONTINUE
721  80 CONTINUE
722  END IF
723 *
724  90 CONTINUE
725 *
726  IF( ierr.NE.0 ) THEN
727  WRITE( nounit, fmt = 9999 )'Generator', ierr, n, jtype,
728  $ ioldsd
729  info = abs( ierr )
730  RETURN
731  END IF
732 *
733  100 CONTINUE
734 *
735  DO 110 i = 1, 7
736  result( i ) = -one
737  110 CONTINUE
738 *
739 * Call CGGEV to compute eigenvalues and eigenvectors.
740 *
741  CALL clacpy( ' ', n, n, a, lda, s, lda )
742  CALL clacpy( ' ', n, n, b, lda, t, lda )
743  CALL cggev( 'V', 'V', n, s, lda, t, lda, alpha, beta, q,
744  $ ldq, z, ldq, work, lwork, rwork, ierr )
745  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
746  result( 1 ) = ulpinv
747  WRITE( nounit, fmt = 9999 )'CGGEV1', ierr, n, jtype,
748  $ ioldsd
749  info = abs( ierr )
750  GO TO 190
751  END IF
752 *
753 * Do the tests (1) and (2)
754 *
755  CALL cget52( .true., n, a, lda, b, lda, q, ldq, alpha, beta,
756  $ work, rwork, result( 1 ) )
757  IF( result( 2 ).GT.thresh ) THEN
758  WRITE( nounit, fmt = 9998 )'Left', 'CGGEV1',
759  $ result( 2 ), n, jtype, ioldsd
760  END IF
761 *
762 * Do the tests (3) and (4)
763 *
764  CALL cget52( .false., n, a, lda, b, lda, z, ldq, alpha,
765  $ beta, work, rwork, result( 3 ) )
766  IF( result( 4 ).GT.thresh ) THEN
767  WRITE( nounit, fmt = 9998 )'Right', 'CGGEV1',
768  $ result( 4 ), n, jtype, ioldsd
769  END IF
770 *
771 * Do test (5)
772 *
773  CALL clacpy( ' ', n, n, a, lda, s, lda )
774  CALL clacpy( ' ', n, n, b, lda, t, lda )
775  CALL cggev( 'N', 'N', n, s, lda, t, lda, alpha1, beta1, q,
776  $ ldq, z, ldq, work, lwork, rwork, ierr )
777  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
778  result( 1 ) = ulpinv
779  WRITE( nounit, fmt = 9999 )'CGGEV2', ierr, n, jtype,
780  $ ioldsd
781  info = abs( ierr )
782  GO TO 190
783  END IF
784 *
785  DO 120 j = 1, n
786  IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
787  $ beta1( j ) )result( 5 ) = ulpinv
788  120 CONTINUE
789 *
790 * Do test (6): Compute eigenvalues and left eigenvectors,
791 * and test them
792 *
793  CALL clacpy( ' ', n, n, a, lda, s, lda )
794  CALL clacpy( ' ', n, n, b, lda, t, lda )
795  CALL cggev( 'V', 'N', n, s, lda, t, lda, alpha1, beta1, qe,
796  $ ldqe, z, ldq, work, lwork, rwork, ierr )
797  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
798  result( 1 ) = ulpinv
799  WRITE( nounit, fmt = 9999 )'CGGEV3', ierr, n, jtype,
800  $ ioldsd
801  info = abs( ierr )
802  GO TO 190
803  END IF
804 *
805  DO 130 j = 1, n
806  IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
807  $ beta1( j ) )result( 6 ) = ulpinv
808  130 CONTINUE
809 *
810  DO 150 j = 1, n
811  DO 140 jc = 1, n
812  IF( q( j, jc ).NE.qe( j, jc ) )
813  $ result( 6 ) = ulpinv
814  140 CONTINUE
815  150 CONTINUE
816 *
817 * Do test (7): Compute eigenvalues and right eigenvectors,
818 * and test them
819 *
820  CALL clacpy( ' ', n, n, a, lda, s, lda )
821  CALL clacpy( ' ', n, n, b, lda, t, lda )
822  CALL cggev( 'N', 'V', n, s, lda, t, lda, alpha1, beta1, q,
823  $ ldq, qe, ldqe, work, lwork, rwork, ierr )
824  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
825  result( 1 ) = ulpinv
826  WRITE( nounit, fmt = 9999 )'CGGEV4', ierr, n, jtype,
827  $ ioldsd
828  info = abs( ierr )
829  GO TO 190
830  END IF
831 *
832  DO 160 j = 1, n
833  IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
834  $ beta1( j ) )result( 7 ) = ulpinv
835  160 CONTINUE
836 *
837  DO 180 j = 1, n
838  DO 170 jc = 1, n
839  IF( z( j, jc ).NE.qe( j, jc ) )
840  $ result( 7 ) = ulpinv
841  170 CONTINUE
842  180 CONTINUE
843 *
844 * End of Loop -- Check for RESULT(j) > THRESH
845 *
846  190 CONTINUE
847 *
848  ntestt = ntestt + 7
849 *
850 * Print out tests which fail.
851 *
852  DO 200 jr = 1, 7
853  IF( result( jr ).GE.thresh ) THEN
854 *
855 * If this is the first test to fail,
856 * print a header to the data file.
857 *
858  IF( nerrs.EQ.0 ) THEN
859  WRITE( nounit, fmt = 9997 )'CGV'
860 *
861 * Matrix types
862 *
863  WRITE( nounit, fmt = 9996 )
864  WRITE( nounit, fmt = 9995 )
865  WRITE( nounit, fmt = 9994 )'Orthogonal'
866 *
867 * Tests performed
868 *
869  WRITE( nounit, fmt = 9993 )
870 *
871  END IF
872  nerrs = nerrs + 1
873  IF( result( jr ).LT.10000.0 ) THEN
874  WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
875  $ result( jr )
876  ELSE
877  WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
878  $ result( jr )
879  END IF
880  END IF
881  200 CONTINUE
882 *
883  210 CONTINUE
884  220 CONTINUE
885 *
886 * Summary
887 *
888  CALL alasvm( 'CGV', nounit, nerrs, ntestt, 0 )
889 *
890  work( 1 ) = maxwrk
891 *
892  RETURN
893 *
894  9999 FORMAT( ' CDRGEV: ', a, ' returned INFO=', i6, '.', / 3x, 'N=',
895  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
896 *
897  9998 FORMAT( ' CDRGEV: ', a, ' Eigenvectors from ', a, ' incorrectly ',
898  $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 3x,
899  $ 'N=', i4, ', JTYPE=', i3, ', ISEED=(', 3( i4, ',' ), i5,
900  $ ')' )
901 *
902  9997 FORMAT( / 1x, a3, ' -- Complex Generalized eigenvalue problem ',
903  $ 'driver' )
904 *
905  9996 FORMAT( ' Matrix types (see CDRGEV for details): ' )
906 *
907  9995 FORMAT( ' Special Matrices:', 23x,
908  $ '(J''=transposed Jordan block)',
909  $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
910  $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
911  $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
912  $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
913  $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
914  $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
915  9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
916  $ / ' 16=Transposed Jordan Blocks 19=geometric ',
917  $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
918  $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
919  $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
920  $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
921  $ '23=(small,large) 24=(small,small) 25=(large,large)',
922  $ / ' 26=random O(1) matrices.' )
923 *
924  9993 FORMAT( / ' Tests performed: ',
925  $ / ' 1 = max | ( b A - a B )''*l | / const.,',
926  $ / ' 2 = | |VR(i)| - 1 | / ulp,',
927  $ / ' 3 = max | ( b A - a B )*r | / const.',
928  $ / ' 4 = | |VL(i)| - 1 | / ulp,',
929  $ / ' 5 = 0 if W same no matter if r or l computed,',
930  $ / ' 6 = 0 if l same no matter if l computed,',
931  $ / ' 7 = 0 if r same no matter if r computed,', / 1x )
932  9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
933  $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
934  9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
935  $ 4( i4, ',' ), ' result ', i2, ' is', 1p, e10.3 )
936 *
937 * End of CDRGEV
938 *
939  END
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
subroutine cunm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
CUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition: cunm2r.f:161
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
subroutine cget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA, WORK, RWORK, RESULT)
CGET52
Definition: cget52.f:163
subroutine clatm4(ITYPE, N, NZ1, NZ2, RSIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
CLATM4
Definition: clatm4.f:173
subroutine cdrgev(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK, RWORK, RESULT, INFO)
CDRGEV
Definition: cdrgev.f:401
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine cggev(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO)
CGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
Definition: cggev.f:219
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108