LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
cdrges3.f
Go to the documentation of this file.
1 *> \brief \b CDRGES3
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 CDRGES3( 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 *> CDRGES3 checks the nonsymmetric generalized eigenvalue (Schur form)
35 *> problem driver CGGES3.
36 *>
37 *> CGGES3 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 CDRGES3 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 *> SDRGES3 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, SDRGES3
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 SDRGES3 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 CGGES3. 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 CGGES3.
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 CGGES3.
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 CGGES3.
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 CGGES3.
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 February 2015
376 *
377 *> \ingroup complex_eig
378 *
379 * =====================================================================
380  SUBROUTINE cdrges3( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
381  $ nounit, a, lda, b, s, t, q, ldq, z, alpha,
382  $ beta, work, lwork, rwork, result, bwork,
383  $ info )
384 *
385 * -- LAPACK test routine (version 3.6.1) --
386 * -- LAPACK is a software package provided by Univ. of Tennessee, --
387 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
388 * February 2015
389 *
390 * .. Scalar Arguments ..
391  INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
392  REAL THRESH
393 * ..
394 * .. Array Arguments ..
395  LOGICAL BWORK( * ), DOTYPE( * )
396  INTEGER ISEED( 4 ), NN( * )
397  REAL RESULT( 13 ), RWORK( * )
398  COMPLEX A( lda, * ), ALPHA( * ), B( lda, * ),
399  $ beta( * ), q( ldq, * ), s( lda, * ),
400  $ t( lda, * ), work( * ), z( ldq, * )
401 * ..
402 *
403 * =====================================================================
404 *
405 * .. Parameters ..
406  REAL ZERO, ONE
407  parameter ( zero = 0.0e+0, one = 1.0e+0 )
408  COMPLEX CZERO, CONE
409  parameter ( czero = ( 0.0e+0, 0.0e+0 ),
410  $ cone = ( 1.0e+0, 0.0e+0 ) )
411  INTEGER MAXTYP
412  parameter ( maxtyp = 26 )
413 * ..
414 * .. Local Scalars ..
415  LOGICAL BADNN, ILABAD
416  CHARACTER SORT
417  INTEGER I, IADD, IINFO, IN, ISORT, J, JC, JR, JSIZE,
418  $ jtype, knteig, maxwrk, minwrk, mtypes, n, n1,
419  $ nb, nerrs, nmats, nmax, ntest, ntestt, rsub,
420  $ sdim
421  REAL SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV
422  COMPLEX CTEMP, X
423 * ..
424 * .. Local Arrays ..
425  LOGICAL LASIGN( maxtyp ), LBSIGN( maxtyp )
426  INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( maxtyp ),
427  $ katype( maxtyp ), kazero( maxtyp ),
428  $ kbmagn( maxtyp ), kbtype( maxtyp ),
429  $ kbzero( maxtyp ), kclass( maxtyp ),
430  $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
431  REAL RMAGN( 0: 3 )
432 * ..
433 * .. External Functions ..
434  LOGICAL CLCTES
435  INTEGER ILAENV
436  REAL SLAMCH
437  COMPLEX CLARND
438  EXTERNAL clctes, ilaenv, slamch, clarnd
439 * ..
440 * .. External Subroutines ..
441  EXTERNAL alasvm, cget51, cget54, cgges3, clacpy, clarfg,
443 * ..
444 * .. Intrinsic Functions ..
445  INTRINSIC abs, aimag, conjg, max, min, REAL, SIGN
446 * ..
447 * .. Statement Functions ..
448  REAL ABS1
449 * ..
450 * .. Statement Function definitions ..
451  abs1( x ) = abs( REAL( X ) ) + abs( AIMAG( x ) )
452 * ..
453 * .. Data statements ..
454  DATA kclass / 15*1, 10*2, 1*3 /
455  DATA kz1 / 0, 1, 2, 1, 3, 3 /
456  DATA kz2 / 0, 0, 1, 2, 1, 1 /
457  DATA kadd / 0, 0, 0, 0, 3, 2 /
458  DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
459  $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
460  DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
461  $ 1, 1, -4, 2, -4, 8*8, 0 /
462  DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
463  $ 4*5, 4*3, 1 /
464  DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
465  $ 4*6, 4*4, 1 /
466  DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
467  $ 2, 1 /
468  DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
469  $ 2, 1 /
470  DATA ktrian / 16*0, 10*1 /
471  DATA lasign / 6*.false., .true., .false., 2*.true.,
472  $ 2*.false., 3*.true., .false., .true.,
473  $ 3*.false., 5*.true., .false. /
474  DATA lbsign / 7*.false., .true., 2*.false.,
475  $ 2*.true., 2*.false., .true., .false., .true.,
476  $ 9*.false. /
477 * ..
478 * .. Executable Statements ..
479 *
480 * Check for errors
481 *
482  info = 0
483 *
484  badnn = .false.
485  nmax = 1
486  DO 10 j = 1, nsizes
487  nmax = max( nmax, nn( j ) )
488  IF( nn( j ).LT.0 )
489  $ badnn = .true.
490  10 CONTINUE
491 *
492  IF( nsizes.LT.0 ) THEN
493  info = -1
494  ELSE IF( badnn ) THEN
495  info = -2
496  ELSE IF( ntypes.LT.0 ) THEN
497  info = -3
498  ELSE IF( thresh.LT.zero ) THEN
499  info = -6
500  ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
501  info = -9
502  ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
503  info = -14
504  END IF
505 *
506 * Compute workspace
507 * (Note: Comments in the code beginning "Workspace:" describe the
508 * minimal amount of workspace needed at that point in the code,
509 * as well as the preferred amount for good performance.
510 * NB refers to the optimal block size for the immediately
511 * following subroutine, as returned by ILAENV.
512 *
513  minwrk = 1
514  IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
515  minwrk = 3*nmax*nmax
516  nb = max( 1, ilaenv( 1, 'CGEQRF', ' ', nmax, nmax, -1, -1 ),
517  $ ilaenv( 1, 'CUNMQR', 'LC', nmax, nmax, nmax, -1 ),
518  $ ilaenv( 1, 'CUNGQR', ' ', nmax, nmax, nmax, -1 ) )
519  maxwrk = max( nmax+nmax*nb, 3*nmax*nmax)
520  work( 1 ) = maxwrk
521  END IF
522 *
523  IF( lwork.LT.minwrk )
524  $ info = -19
525 *
526  IF( info.NE.0 ) THEN
527  CALL xerbla( 'CDRGES3', -info )
528  RETURN
529  END IF
530 *
531 * Quick return if possible
532 *
533  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
534  $ RETURN
535 *
536  ulp = slamch( 'Precision' )
537  safmin = slamch( 'Safe minimum' )
538  safmin = safmin / ulp
539  safmax = one / safmin
540  CALL slabad( safmin, safmax )
541  ulpinv = one / ulp
542 *
543 * The values RMAGN(2:3) depend on N, see below.
544 *
545  rmagn( 0 ) = zero
546  rmagn( 1 ) = one
547 *
548 * Loop over matrix sizes
549 *
550  ntestt = 0
551  nerrs = 0
552  nmats = 0
553 *
554  DO 190 jsize = 1, nsizes
555  n = nn( jsize )
556  n1 = max( 1, n )
557  rmagn( 2 ) = safmax*ulp / REAL( n1 )
558  rmagn( 3 ) = safmin*ulpinv*REAL( n1 )
559 *
560  IF( nsizes.NE.1 ) THEN
561  mtypes = min( maxtyp, ntypes )
562  ELSE
563  mtypes = min( maxtyp+1, ntypes )
564  END IF
565 *
566 * Loop over matrix types
567 *
568  DO 180 jtype = 1, mtypes
569  IF( .NOT.dotype( jtype ) )
570  $ GO TO 180
571  nmats = nmats + 1
572  ntest = 0
573 *
574 * Save ISEED in case of an error.
575 *
576  DO 20 j = 1, 4
577  ioldsd( j ) = iseed( j )
578  20 CONTINUE
579 *
580 * Initialize RESULT
581 *
582  DO 30 j = 1, 13
583  result( j ) = zero
584  30 CONTINUE
585 *
586 * Generate test matrices A and B
587 *
588 * Description of control parameters:
589 *
590 * KCLASS: =1 means w/o rotation, =2 means w/ rotation,
591 * =3 means random.
592 * KATYPE: the "type" to be passed to CLATM4 for computing A.
593 * KAZERO: the pattern of zeros on the diagonal for A:
594 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
595 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
596 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
597 * non-zero entries.)
598 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
599 * =2: large, =3: small.
600 * LASIGN: .TRUE. if the diagonal elements of A are to be
601 * multiplied by a random magnitude 1 number.
602 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
603 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
604 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
605 * RMAGN: used to implement KAMAGN and KBMAGN.
606 *
607  IF( mtypes.GT.maxtyp )
608  $ GO TO 110
609  iinfo = 0
610  IF( kclass( jtype ).LT.3 ) THEN
611 *
612 * Generate A (w/o rotation)
613 *
614  IF( abs( katype( jtype ) ).EQ.3 ) THEN
615  in = 2*( ( n-1 ) / 2 ) + 1
616  IF( in.NE.n )
617  $ CALL claset( 'Full', n, n, czero, czero, a, lda )
618  ELSE
619  in = n
620  END IF
621  CALL clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
622  $ kz2( kazero( jtype ) ), lasign( jtype ),
623  $ rmagn( kamagn( jtype ) ), ulp,
624  $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
625  $ iseed, a, lda )
626  iadd = kadd( kazero( jtype ) )
627  IF( iadd.GT.0 .AND. iadd.LE.n )
628  $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
629 *
630 * Generate B (w/o rotation)
631 *
632  IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
633  in = 2*( ( n-1 ) / 2 ) + 1
634  IF( in.NE.n )
635  $ CALL claset( 'Full', n, n, czero, czero, b, lda )
636  ELSE
637  in = n
638  END IF
639  CALL clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
640  $ kz2( kbzero( jtype ) ), lbsign( jtype ),
641  $ rmagn( kbmagn( jtype ) ), one,
642  $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
643  $ iseed, b, lda )
644  iadd = kadd( kbzero( jtype ) )
645  IF( iadd.NE.0 .AND. iadd.LE.n )
646  $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
647 *
648  IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
649 *
650 * Include rotations
651 *
652 * Generate Q, Z as Householder transformations times
653 * a diagonal matrix.
654 *
655  DO 50 jc = 1, n - 1
656  DO 40 jr = jc, n
657  q( jr, jc ) = clarnd( 3, iseed )
658  z( jr, jc ) = clarnd( 3, iseed )
659  40 CONTINUE
660  CALL clarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
661  $ work( jc ) )
662  work( 2*n+jc ) = sign( one, REAL( Q( JC, JC ) ) )
663  q( jc, jc ) = cone
664  CALL clarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
665  $ work( n+jc ) )
666  work( 3*n+jc ) = sign( one, REAL( Z( JC, JC ) ) )
667  z( jc, jc ) = cone
668  50 CONTINUE
669  ctemp = clarnd( 3, iseed )
670  q( n, n ) = cone
671  work( n ) = czero
672  work( 3*n ) = ctemp / abs( ctemp )
673  ctemp = clarnd( 3, iseed )
674  z( n, n ) = cone
675  work( 2*n ) = czero
676  work( 4*n ) = ctemp / abs( ctemp )
677 *
678 * Apply the diagonal matrices
679 *
680  DO 70 jc = 1, n
681  DO 60 jr = 1, n
682  a( jr, jc ) = work( 2*n+jr )*
683  $ conjg( work( 3*n+jc ) )*
684  $ a( jr, jc )
685  b( jr, jc ) = work( 2*n+jr )*
686  $ conjg( work( 3*n+jc ) )*
687  $ b( jr, jc )
688  60 CONTINUE
689  70 CONTINUE
690  CALL cunm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
691  $ lda, work( 2*n+1 ), iinfo )
692  IF( iinfo.NE.0 )
693  $ GO TO 100
694  CALL cunm2r( 'R', 'C', n, n, n-1, z, ldq, work( n+1 ),
695  $ a, lda, work( 2*n+1 ), iinfo )
696  IF( iinfo.NE.0 )
697  $ GO TO 100
698  CALL cunm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
699  $ lda, work( 2*n+1 ), iinfo )
700  IF( iinfo.NE.0 )
701  $ GO TO 100
702  CALL cunm2r( 'R', 'C', n, n, n-1, z, ldq, work( n+1 ),
703  $ b, lda, work( 2*n+1 ), iinfo )
704  IF( iinfo.NE.0 )
705  $ GO TO 100
706  END IF
707  ELSE
708 *
709 * Random matrices
710 *
711  DO 90 jc = 1, n
712  DO 80 jr = 1, n
713  a( jr, jc ) = rmagn( kamagn( jtype ) )*
714  $ clarnd( 4, iseed )
715  b( jr, jc ) = rmagn( kbmagn( jtype ) )*
716  $ clarnd( 4, iseed )
717  80 CONTINUE
718  90 CONTINUE
719  END IF
720 *
721  100 CONTINUE
722 *
723  IF( iinfo.NE.0 ) THEN
724  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
725  $ ioldsd
726  info = abs( iinfo )
727  RETURN
728  END IF
729 *
730  110 CONTINUE
731 *
732  DO 120 i = 1, 13
733  result( i ) = -one
734  120 CONTINUE
735 *
736 * Test with and without sorting of eigenvalues
737 *
738  DO 150 isort = 0, 1
739  IF( isort.EQ.0 ) THEN
740  sort = 'N'
741  rsub = 0
742  ELSE
743  sort = 'S'
744  rsub = 5
745  END IF
746 *
747 * Call CGGES3 to compute H, T, Q, Z, alpha, and beta.
748 *
749  CALL clacpy( 'Full', n, n, a, lda, s, lda )
750  CALL clacpy( 'Full', n, n, b, lda, t, lda )
751  ntest = 1 + rsub + isort
752  result( 1+rsub+isort ) = ulpinv
753  CALL cgges3( 'V', 'V', sort, clctes, n, s, lda, t, lda,
754  $ sdim, alpha, beta, q, ldq, z, ldq, work,
755  $ lwork, rwork, bwork, iinfo )
756  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
757  result( 1+rsub+isort ) = ulpinv
758  WRITE( nounit, fmt = 9999 )'CGGES3', iinfo, n, jtype,
759  $ ioldsd
760  info = abs( iinfo )
761  GO TO 160
762  END IF
763 *
764  ntest = 4 + rsub
765 *
766 * Do tests 1--4 (or tests 7--9 when reordering )
767 *
768  IF( isort.EQ.0 ) THEN
769  CALL cget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
770  $ work, rwork, result( 1 ) )
771  CALL cget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
772  $ work, rwork, result( 2 ) )
773  ELSE
774  CALL cget54( n, a, lda, b, lda, s, lda, t, lda, q,
775  $ ldq, z, ldq, work, result( 2+rsub ) )
776  END IF
777 *
778  CALL cget51( 3, n, b, lda, t, lda, q, ldq, q, ldq, work,
779  $ rwork, result( 3+rsub ) )
780  CALL cget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
781  $ rwork, result( 4+rsub ) )
782 *
783 * Do test 5 and 6 (or Tests 10 and 11 when reordering):
784 * check Schur form of A and compare eigenvalues with
785 * diagonals.
786 *
787  ntest = 6 + rsub
788  temp1 = zero
789 *
790  DO 130 j = 1, n
791  ilabad = .false.
792  temp2 = ( abs1( alpha( j )-s( j, j ) ) /
793  $ max( safmin, abs1( alpha( j ) ), abs1( s( j,
794  $ j ) ) )+abs1( beta( j )-t( j, j ) ) /
795  $ max( safmin, abs1( beta( j ) ), abs1( t( j,
796  $ j ) ) ) ) / ulp
797 *
798  IF( j.LT.n ) THEN
799  IF( s( j+1, j ).NE.zero ) THEN
800  ilabad = .true.
801  result( 5+rsub ) = ulpinv
802  END IF
803  END IF
804  IF( j.GT.1 ) THEN
805  IF( s( j, j-1 ).NE.zero ) THEN
806  ilabad = .true.
807  result( 5+rsub ) = ulpinv
808  END IF
809  END IF
810  temp1 = max( temp1, temp2 )
811  IF( ilabad ) THEN
812  WRITE( nounit, fmt = 9998 )j, n, jtype, ioldsd
813  END IF
814  130 CONTINUE
815  result( 6+rsub ) = temp1
816 *
817  IF( isort.GE.1 ) THEN
818 *
819 * Do test 12
820 *
821  ntest = 12
822  result( 12 ) = zero
823  knteig = 0
824  DO 140 i = 1, n
825  IF( clctes( alpha( i ), beta( i ) ) )
826  $ knteig = knteig + 1
827  140 CONTINUE
828  IF( sdim.NE.knteig )
829  $ result( 13 ) = ulpinv
830  END IF
831 *
832  150 CONTINUE
833 *
834 * End of Loop -- Check for RESULT(j) > THRESH
835 *
836  160 CONTINUE
837 *
838  ntestt = ntestt + ntest
839 *
840 * Print out tests which fail.
841 *
842  DO 170 jr = 1, ntest
843  IF( result( jr ).GE.thresh ) THEN
844 *
845 * If this is the first test to fail,
846 * print a header to the data file.
847 *
848  IF( nerrs.EQ.0 ) THEN
849  WRITE( nounit, fmt = 9997 )'CGS'
850 *
851 * Matrix types
852 *
853  WRITE( nounit, fmt = 9996 )
854  WRITE( nounit, fmt = 9995 )
855  WRITE( nounit, fmt = 9994 )'Unitary'
856 *
857 * Tests performed
858 *
859  WRITE( nounit, fmt = 9993 )'unitary', '''',
860  $ 'transpose', ( '''', j = 1, 8 )
861 *
862  END IF
863  nerrs = nerrs + 1
864  IF( result( jr ).LT.10000.0 ) THEN
865  WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
866  $ result( jr )
867  ELSE
868  WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
869  $ result( jr )
870  END IF
871  END IF
872  170 CONTINUE
873 *
874  180 CONTINUE
875  190 CONTINUE
876 *
877 * Summary
878 *
879  CALL alasvm( 'CGS', nounit, nerrs, ntestt, 0 )
880 *
881  work( 1 ) = maxwrk
882 *
883  RETURN
884 *
885  9999 FORMAT( ' CDRGES3: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
886  $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
887 *
888  9998 FORMAT( ' CDRGES3: S not in Schur form at eigenvalue ', i6, '.',
889  $ / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
890  $ i5, ')' )
891 *
892  9997 FORMAT( / 1x, a3, ' -- Complex Generalized Schur from problem ',
893  $ 'driver' )
894 *
895  9996 FORMAT( ' Matrix types (see CDRGES3 for details): ' )
896 *
897  9995 FORMAT( ' Special Matrices:', 23x,
898  $ '(J''=transposed Jordan block)',
899  $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
900  $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
901  $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
902  $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
903  $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
904  $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
905  9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
906  $ / ' 16=Transposed Jordan Blocks 19=geometric ',
907  $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
908  $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
909  $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
910  $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
911  $ '23=(small,large) 24=(small,small) 25=(large,large)',
912  $ / ' 26=random O(1) matrices.' )
913 *
914  9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
915  $ 'Q and Z are ', a, ',', / 19x,
916  $ 'l and r are the appropriate left and right', / 19x,
917  $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
918  $ ' means ', a, '.)', / ' Without ordering: ',
919  $ / ' 1 = | A - Q S Z', a,
920  $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
921  $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
922  $ ' | / ( n ulp ) 4 = | I - ZZ', a,
923  $ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
924  $ / ' 6 = difference between (alpha,beta)',
925  $ ' and diagonals of (S,T)', / ' With ordering: ',
926  $ / ' 7 = | (A,B) - Q (S,T) Z', a, ' | / ( |(A,B)| n ulp )',
927  $ / ' 8 = | I - QQ', a,
928  $ ' | / ( n ulp ) 9 = | I - ZZ', a,
929  $ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
930  $ / ' 11 = difference between (alpha,beta) and diagonals',
931  $ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
932  $ 'selected eigenvalues', / )
933  9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
934  $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
935  9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
936  $ 4( i4, ',' ), ' result ', i2, ' is', 1p, e10.3 )
937 *
938 * End of CDRGES3
939 *
940  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 cgges3(JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK, BWORK, INFO)
CGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
Definition: cgges3.f:271
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 cdrges3(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHA, BETA, WORK, LWORK, RWORK, RESULT, BWORK, INFO)
CDRGES3
Definition: cdrges3.f:384
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 clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108