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