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