LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 threshhold 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 November 2011
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.4.0) --
385 * -- LAPACK is a software package provided by Univ. of Tennessee, --
386 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
387 * November 2011
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