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