LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
cdrgev3.f
Go to the documentation of this file.
1 *> \brief \b CDRGEV3
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 CDRGEV3( 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 *> CDRGEV3 checks the nonsymmetric generalized eigenvalue problem driver
38 *> routine CGGEV3.
39 *>
40 *> CGGEV3 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 CDRGEV3 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 CGGEV3:
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 *> CDRGEV3 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, CDRGEV3
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 CDRGEV3 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 CGGEV3. 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 CGGEV3.
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 CGGEV3.
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 CGGEV3.
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 CGGEV3.
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 CGGEV3 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 January 2015
393 *
394 *> \ingroup complex_eig
395 *
396 * =====================================================================
397  SUBROUTINE cdrgev3( 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,
400  $ RWORK, 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 * January 2015
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, cggev3, clacpy, clarfg, claset,
458  $ clatm4, cunm2r, slabad, xerbla
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( 'CDRGEV3', -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 CGGEV3 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 cggev3( '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 )'CGGEV31', 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', 'CGGEV31',
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', 'CGGEV31',
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 cggev3( '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 )'CGGEV32', 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 the 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 cggev3( '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 )'CGGEV33', ierr, n, jtype,
800  $ ioldsd
801  info = abs( ierr )
802  GO TO 190
803  END IF
804 
805 *
806  DO 130 j = 1, n
807  IF( alpha( j ).NE.alpha1( j ) .OR.
808  $ beta( j ).NE.beta1( j ) ) THEN
809  result( 6 ) = ulpinv
810  ENDIF
811  130 CONTINUE
812 *
813  DO 150 j = 1, n
814  DO 140 jc = 1, n
815  IF( q( j, jc ).NE.qe( j, jc ) ) THEN
816  result( 6 ) = ulpinv
817  END IF
818  140 CONTINUE
819  150 CONTINUE
820 *
821 * DO the test (7): Compute eigenvalues and right eigenvectors,
822 * and test them
823 *
824  CALL clacpy( ' ', n, n, a, lda, s, lda )
825  CALL clacpy( ' ', n, n, b, lda, t, lda )
826  CALL cggev3( 'N', 'V', n, s, lda, t, lda, alpha1, beta1, q,
827  $ ldq, qe, ldqe, work, lwork, rwork, ierr )
828  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
829  result( 1 ) = ulpinv
830  WRITE( nounit, fmt = 9999 )'CGGEV34', ierr, n, jtype,
831  $ ioldsd
832  info = abs( ierr )
833  GO TO 190
834  END IF
835 *
836  DO 160 j = 1, n
837  IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
838  $ beta1( j ) )result( 7 ) = ulpinv
839  160 CONTINUE
840 *
841  DO 180 j = 1, n
842  DO 170 jc = 1, n
843  IF( z( j, jc ).NE.qe( j, jc ) )
844  $ result( 7 ) = ulpinv
845  170 CONTINUE
846  180 CONTINUE
847 *
848 * End of Loop -- Check for RESULT(j) > THRESH
849 *
850  190 CONTINUE
851 *
852  ntestt = ntestt + 7
853 *
854 * Print out tests which fail.
855 *
856  DO 200 jr = 1, 7
857  IF( result( jr ).GE.thresh ) THEN
858 *
859 * If this is the first test to fail,
860 * print a header to the data file.
861 *
862  IF( nerrs.EQ.0 ) THEN
863  WRITE( nounit, fmt = 9997 )'CGV'
864 *
865 * Matrix types
866 *
867  WRITE( nounit, fmt = 9996 )
868  WRITE( nounit, fmt = 9995 )
869  WRITE( nounit, fmt = 9994 )'Orthogonal'
870 *
871 * Tests performed
872 *
873  WRITE( nounit, fmt = 9993 )
874 *
875  END IF
876  nerrs = nerrs + 1
877  IF( result( jr ).LT.10000.0 ) THEN
878  WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
879  $ result( jr )
880  ELSE
881  WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
882  $ result( jr )
883  END IF
884  END IF
885  200 CONTINUE
886 *
887  210 CONTINUE
888  220 CONTINUE
889 *
890 * Summary
891 *
892  CALL alasvm( 'CGV3', nounit, nerrs, ntestt, 0 )
893 *
894  work( 1 ) = maxwrk
895 *
896  RETURN
897 *
898  9999 FORMAT( ' CDRGEV3: ', a, ' returned INFO=', i6, '.', / 3x, 'N=',
899  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
900 *
901  9998 FORMAT( ' CDRGEV3: ', a, ' Eigenvectors from ', a,
902  $ ' incorrectly normalized.', / ' Bits of error=', 0p, g10.3,
903  $ ',', 3x, 'N=', i4, ', JTYPE=', i3, ', ISEED=(',
904  $ 3( i4, ',' ), i5, ')' )
905 *
906  9997 FORMAT( / 1x, a3, ' -- Complex Generalized eigenvalue problem ',
907  $ 'driver' )
908 *
909  9996 FORMAT( ' Matrix types (see CDRGEV3 for details): ' )
910 *
911  9995 FORMAT( ' Special Matrices:', 23x,
912  $ '(J''=transposed Jordan block)',
913  $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
914  $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
915  $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
916  $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
917  $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
918  $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
919  9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
920  $ / ' 16=Transposed Jordan Blocks 19=geometric ',
921  $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
922  $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
923  $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
924  $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
925  $ '23=(small,large) 24=(small,small) 25=(large,large)',
926  $ / ' 26=random O(1) matrices.' )
927 *
928  9993 FORMAT( / ' Tests performed: ',
929  $ / ' 1 = max | ( b A - a B )''*l | / const.,',
930  $ / ' 2 = | |VR(i)| - 1 | / ulp,',
931  $ / ' 3 = max | ( b A - a B )*r | / const.',
932  $ / ' 4 = | |VL(i)| - 1 | / ulp,',
933  $ / ' 5 = 0 if W same no matter if r or l computed,',
934  $ / ' 6 = 0 if l same no matter if l computed,',
935  $ / ' 7 = 0 if r same no matter if r computed,', / 1x )
936  9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
937  $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
938  9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
939  $ 4( i4, ',' ), ' result ', i2, ' is', 1p, e10.3 )
940 *
941 * End of CDRGEV3
942 *
943  END
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
clarfg
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108
alasvm
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
clacpy
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
claset
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
clatm4
subroutine clatm4(ITYPE, N, NZ1, NZ2, RSIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
CLATM4
Definition: clatm4.f:173
cdrgev3
subroutine cdrgev3(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)
CDRGEV3
Definition: cdrgev3.f:401
cunm2r
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
cget52
subroutine cget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA, WORK, RWORK, RESULT)
CGET52
Definition: cget52.f:163
cggev3
subroutine cggev3(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO)
CGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (...
Definition: cggev3.f:218