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