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