LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
sdrvgg.f
Go to the documentation of this file.
1 *> \brief \b SDRVGG
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 SDRVGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * THRSHN, NOUNIT, A, LDA, B, S, T, S2, T2, Q,
13 * LDQ, Z, ALPHR1, ALPHI1, BETA1, ALPHR2, ALPHI2,
14 * BETA2, VL, VR, WORK, LWORK, RESULT, INFO )
15 *
16 * .. Scalar Arguments ..
17 * INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
18 * REAL THRESH, THRSHN
19 * ..
20 * .. Array Arguments ..
21 * LOGICAL DOTYPE( * )
22 * INTEGER ISEED( 4 ), NN( * )
23 * REAL A( LDA, * ), ALPHI1( * ), ALPHI2( * ),
24 * $ ALPHR1( * ), ALPHR2( * ), B( LDA, * ),
25 * $ BETA1( * ), BETA2( * ), Q( LDQ, * ),
26 * $ RESULT( * ), S( LDA, * ), S2( LDA, * ),
27 * $ T( LDA, * ), T2( LDA, * ), VL( LDQ, * ),
28 * $ VR( LDQ, * ), WORK( * ), Z( LDQ, * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> SDRVGG checks the nonsymmetric generalized eigenvalue driver
38 *> routines.
39 *> T T T
40 *> SGEGS factors A and B as Q S Z and Q T Z , where means
41 *> transpose, T is upper triangular, S is in generalized Schur form
42 *> (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
43 *> the 2x2 blocks corresponding to complex conjugate pairs of
44 *> generalized eigenvalues), and Q and Z are orthogonal. It also
45 *> computes the generalized eigenvalues (alpha(1),beta(1)), ...,
46 *> (alpha(n),beta(n)), where alpha(j)=S(j,j) and beta(j)=P(j,j) --
47 *> thus, w(j) = alpha(j)/beta(j) is a root of the generalized
48 *> eigenvalue problem
49 *>
50 *> det( A - w(j) B ) = 0
51 *>
52 *> and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent
53 *> problem
54 *>
55 *> det( m(j) A - B ) = 0
56 *>
57 *> SGEGV computes the generalized eigenvalues (alpha(1),beta(1)), ...,
58 *> (alpha(n),beta(n)), the matrix L whose columns contain the
59 *> generalized left eigenvectors l, and the matrix R whose columns
60 *> contain the generalized right eigenvectors r for the pair (A,B).
61 *>
62 *> When SDRVGG is called, a number of matrix "sizes" ("n's") and a
63 *> number of matrix "types" are specified. For each size ("n")
64 *> and each type of matrix, one matrix will be generated and used
65 *> to test the nonsymmetric eigenroutines. For each matrix, 7
66 *> tests will be performed and compared with the threshhold THRESH:
67 *>
68 *> Results from SGEGS:
69 *>
70 *> T
71 *> (1) | A - Q S Z | / ( |A| n ulp )
72 *>
73 *> T
74 *> (2) | B - Q T Z | / ( |B| n ulp )
75 *>
76 *> T
77 *> (3) | I - QQ | / ( n ulp )
78 *>
79 *> T
80 *> (4) | I - ZZ | / ( n ulp )
81 *>
82 *> (5) maximum over j of D(j) where:
83 *>
84 *> if alpha(j) is real:
85 *> |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
86 *> D(j) = ------------------------ + -----------------------
87 *> max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
88 *>
89 *> if alpha(j) is complex:
90 *> | det( s S - w T ) |
91 *> D(j) = ---------------------------------------------------
92 *> ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
93 *>
94 *> and S and T are here the 2 x 2 diagonal blocks of S and T
95 *> corresponding to the j-th eigenvalue.
96 *>
97 *> Results from SGEGV:
98 *>
99 *> (6) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
100 *>
101 *> | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) )
102 *>
103 *> where l**H is the conjugate tranpose of l.
104 *>
105 *> (7) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
106 *>
107 *> | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
108 *>
109 *> Test Matrices
110 *> ---- --------
111 *>
112 *> The sizes of the test matrices are specified by an array
113 *> NN(1:NSIZES); the value of each element NN(j) specifies one size.
114 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
115 *> DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
116 *> Currently, the list of possible types is:
117 *>
118 *> (1) ( 0, 0 ) (a pair of zero matrices)
119 *>
120 *> (2) ( I, 0 ) (an identity and a zero matrix)
121 *>
122 *> (3) ( 0, I ) (an identity and a zero matrix)
123 *>
124 *> (4) ( I, I ) (a pair of identity matrices)
125 *>
126 *> t t
127 *> (5) ( J , J ) (a pair of transposed Jordan blocks)
128 *>
129 *> t ( I 0 )
130 *> (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
131 *> ( 0 I ) ( 0 J )
132 *> and I is a k x k identity and J a (k+1)x(k+1)
133 *> Jordan block; k=(N-1)/2
134 *>
135 *> (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal
136 *> matrix with those diagonal entries.)
137 *> (8) ( I, D )
138 *>
139 *> (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
140 *>
141 *> (10) ( small*D, big*I )
142 *>
143 *> (11) ( big*I, small*D )
144 *>
145 *> (12) ( small*I, big*D )
146 *>
147 *> (13) ( big*D, big*I )
148 *>
149 *> (14) ( small*D, small*I )
150 *>
151 *> (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
152 *> D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
153 *> t t
154 *> (16) Q ( J , J ) Z where Q and Z are random orthogonal matrices.
155 *>
156 *> (17) Q ( T1, T2 ) Z where T1 and T2 are upper triangular matrices
157 *> with random O(1) entries above the diagonal
158 *> and diagonal entries diag(T1) =
159 *> ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
160 *> ( 0, N-3, N-4,..., 1, 0, 0 )
161 *>
162 *> (18) Q ( T1, T2 ) Z diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
163 *> diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
164 *> s = machine precision.
165 *>
166 *> (19) Q ( T1, T2 ) Z diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
167 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
168 *>
169 *> N-5
170 *> (20) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
171 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
172 *>
173 *> (21) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
174 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
175 *> where r1,..., r(N-4) are random.
176 *>
177 *> (22) Q ( big*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
178 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
179 *>
180 *> (23) Q ( small*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
181 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
182 *>
183 *> (24) Q ( small*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
184 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
185 *>
186 *> (25) Q ( big*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
187 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
188 *>
189 *> (26) Q ( T1, T2 ) Z where T1 and T2 are random upper-triangular
190 *> matrices.
191 *> \endverbatim
192 *
193 * Arguments:
194 * ==========
195 *
196 *> \param[in] NSIZES
197 *> \verbatim
198 *> NSIZES is INTEGER
199 *> The number of sizes of matrices to use. If it is zero,
200 *> SDRVGG does nothing. It must be at least zero.
201 *> \endverbatim
202 *>
203 *> \param[in] NN
204 *> \verbatim
205 *> NN is INTEGER array, dimension (NSIZES)
206 *> An array containing the sizes to be used for the matrices.
207 *> Zero values will be skipped. The values must be at least
208 *> zero.
209 *> \endverbatim
210 *>
211 *> \param[in] NTYPES
212 *> \verbatim
213 *> NTYPES is INTEGER
214 *> The number of elements in DOTYPE. If it is zero, SDRVGG
215 *> does nothing. It must be at least zero. If it is MAXTYP+1
216 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
217 *> defined, which is to use whatever matrix is in A. This
218 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
219 *> DOTYPE(MAXTYP+1) is .TRUE. .
220 *> \endverbatim
221 *>
222 *> \param[in] DOTYPE
223 *> \verbatim
224 *> DOTYPE is LOGICAL array, dimension (NTYPES)
225 *> If DOTYPE(j) is .TRUE., then for each size in NN a
226 *> matrix of that size and of type j will be generated.
227 *> If NTYPES is smaller than the maximum number of types
228 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
229 *> MAXTYP will not be generated. If NTYPES is larger
230 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
231 *> will be ignored.
232 *> \endverbatim
233 *>
234 *> \param[in,out] ISEED
235 *> \verbatim
236 *> ISEED is INTEGER array, dimension (4)
237 *> On entry ISEED specifies the seed of the random number
238 *> generator. The array elements should be between 0 and 4095;
239 *> if not they will be reduced mod 4096. Also, ISEED(4) must
240 *> be odd. The random number generator uses a linear
241 *> congruential sequence limited to small integers, and so
242 *> should produce machine independent random numbers. The
243 *> values of ISEED are changed on exit, and can be used in the
244 *> next call to SDRVGG to continue the same random number
245 *> sequence.
246 *> \endverbatim
247 *>
248 *> \param[in] THRESH
249 *> \verbatim
250 *> THRESH is REAL
251 *> A test will count as "failed" if the "error", computed as
252 *> described above, exceeds THRESH. Note that the error is
253 *> scaled to be O(1), so THRESH should be a reasonably small
254 *> multiple of 1, e.g., 10 or 100. In particular, it should
255 *> not depend on the precision (single vs. double) or the size
256 *> of the matrix. It must be at least zero.
257 *> \endverbatim
258 *>
259 *> \param[in] THRSHN
260 *> \verbatim
261 *> THRSHN is REAL
262 *> Threshhold for reporting eigenvector normalization error.
263 *> If the normalization of any eigenvector differs from 1 by
264 *> more than THRSHN*ulp, then a special error message will be
265 *> printed. (This is handled separately from the other tests,
266 *> since only a compiler or programming error should cause an
267 *> error message, at least if THRSHN is at least 5--10.)
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, dimension
280 *> (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, T, S2, and T2.
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, dimension
296 *> (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 SGEGS. 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 SGEGS.
314 *> \endverbatim
315 *>
316 *> \param[out] S2
317 *> \verbatim
318 *> S2 is REAL array, dimension (LDA, max(NN))
319 *> The matrix computed from A by SGEGV. This will be the
320 *> Schur form of some matrix related to A, but will not, in
321 *> general, be the same as S.
322 *> \endverbatim
323 *>
324 *> \param[out] T2
325 *> \verbatim
326 *> T2 is REAL array, dimension (LDA, max(NN))
327 *> The matrix computed from B by SGEGV. This will be the
328 *> Schur form of some matrix related to B, but will not, in
329 *> general, be the same as T.
330 *> \endverbatim
331 *>
332 *> \param[out] Q
333 *> \verbatim
334 *> Q is REAL array, dimension (LDQ, max(NN))
335 *> The (left) orthogonal matrix computed by SGEGS.
336 *> \endverbatim
337 *>
338 *> \param[in] LDQ
339 *> \verbatim
340 *> LDQ is INTEGER
341 *> The leading dimension of Q, Z, VL, and VR. It must
342 *> be at least 1 and at least max( NN ).
343 *> \endverbatim
344 *>
345 *> \param[out] Z
346 *> \verbatim
347 *> Z is REAL array of
348 *> dimension( LDQ, max(NN) )
349 *> The (right) orthogonal matrix computed by SGEGS.
350 *> \endverbatim
351 *>
352 *> \param[out] ALPHR1
353 *> \verbatim
354 *> ALPHR1 is REAL array, dimension (max(NN))
355 *> \endverbatim
356 *>
357 *> \param[out] ALPHI1
358 *> \verbatim
359 *> ALPHI1 is REAL array, dimension (max(NN))
360 *> \endverbatim
361 *>
362 *> \param[out] BETA1
363 *> \verbatim
364 *> BETA1 is REAL array, dimension (max(NN))
365 *>
366 *> The generalized eigenvalues of (A,B) computed by SGEGS.
367 *> ( ALPHR1(k)+ALPHI1(k)*i ) / BETA1(k) is the k-th
368 *> generalized eigenvalue of the matrices in A and B.
369 *> \endverbatim
370 *>
371 *> \param[out] ALPHR2
372 *> \verbatim
373 *> ALPHR2 is REAL array, dimension (max(NN))
374 *> \endverbatim
375 *>
376 *> \param[out] ALPHI2
377 *> \verbatim
378 *> ALPHI2 is REAL array, dimension (max(NN))
379 *> \endverbatim
380 *>
381 *> \param[out] BETA2
382 *> \verbatim
383 *> BETA2 is REAL array, dimension (max(NN))
384 *>
385 *> The generalized eigenvalues of (A,B) computed by SGEGV.
386 *> ( ALPHR2(k)+ALPHI2(k)*i ) / BETA2(k) is the k-th
387 *> generalized eigenvalue of the matrices in A and B.
388 *> \endverbatim
389 *>
390 *> \param[out] VL
391 *> \verbatim
392 *> VL is REAL array, dimension (LDQ, max(NN))
393 *> The (block lower triangular) left eigenvector matrix for
394 *> the matrices in A and B. (See STGEVC for the format.)
395 *> \endverbatim
396 *>
397 *> \param[out] VR
398 *> \verbatim
399 *> VR is REAL array, dimension (LDQ, max(NN))
400 *> The (block upper triangular) right eigenvector matrix for
401 *> the matrices in A and B. (See STGEVC for the format.)
402 *> \endverbatim
403 *>
404 *> \param[out] WORK
405 *> \verbatim
406 *> WORK is REAL array, dimension (LWORK)
407 *> \endverbatim
408 *>
409 *> \param[in] LWORK
410 *> \verbatim
411 *> LWORK is INTEGER
412 *> The number of entries in WORK. This must be at least
413 *> 2*N + MAX( 6*N, N*(NB+1), (k+1)*(2*k+N+1) ), where
414 *> "k" is the sum of the blocksize and number-of-shifts for
415 *> SHGEQZ, and NB is the greatest of the blocksizes for
416 *> SGEQRF, SORMQR, and SORGQR. (The blocksizes and the
417 *> number-of-shifts are retrieved through calls to ILAENV.)
418 *> \endverbatim
419 *>
420 *> \param[out] RESULT
421 *> \verbatim
422 *> RESULT is REAL array, dimension (15)
423 *> The values computed by the tests described above.
424 *> The values are currently limited to 1/ulp, to avoid
425 *> overflow.
426 *> \endverbatim
427 *>
428 *> \param[out] INFO
429 *> \verbatim
430 *> INFO is INTEGER
431 *> = 0: successful exit
432 *> < 0: if INFO = -i, the i-th argument had an illegal value.
433 *> > 0: A routine returned an error code. INFO is the
434 *> absolute value of the INFO value returned.
435 *> \endverbatim
436 *
437 * Authors:
438 * ========
439 *
440 *> \author Univ. of Tennessee
441 *> \author Univ. of California Berkeley
442 *> \author Univ. of Colorado Denver
443 *> \author NAG Ltd.
444 *
445 *> \date November 2011
446 *
447 *> \ingroup single_eig
448 *
449 * =====================================================================
450  SUBROUTINE sdrvgg( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
451  $ thrshn, nounit, a, lda, b, s, t, s2, t2, q,
452  $ ldq, z, alphr1, alphi1, beta1, alphr2, alphi2,
453  $ beta2, vl, vr, work, lwork, result, info )
454 *
455 * -- LAPACK test routine (version 3.4.0) --
456 * -- LAPACK is a software package provided by Univ. of Tennessee, --
457 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
458 * November 2011
459 *
460 * .. Scalar Arguments ..
461  INTEGER info, lda, ldq, lwork, nounit, nsizes, ntypes
462  REAL thresh, thrshn
463 * ..
464 * .. Array Arguments ..
465  LOGICAL dotype( * )
466  INTEGER iseed( 4 ), nn( * )
467  REAL a( lda, * ), alphi1( * ), alphi2( * ),
468  $ alphr1( * ), alphr2( * ), b( lda, * ),
469  $ beta1( * ), beta2( * ), q( ldq, * ),
470  $ result( * ), s( lda, * ), s2( lda, * ),
471  $ t( lda, * ), t2( lda, * ), vl( ldq, * ),
472  $ vr( ldq, * ), work( * ), z( ldq, * )
473 * ..
474 *
475 * =====================================================================
476 *
477 * .. Parameters ..
478  REAL zero, one
479  parameter( zero = 0.0, one = 1.0 )
480  INTEGER maxtyp
481  parameter( maxtyp = 26 )
482 * ..
483 * .. Local Scalars ..
484  LOGICAL badnn, ilabad
485  INTEGER i1, iadd, iinfo, in, j, jc, jr, jsize, jtype,
486  $ lwkopt, mtypes, n, n1, nb, nbz, nerrs, nmats,
487  $ nmax, ns, ntest, ntestt
488  REAL safmax, safmin, temp1, temp2, ulp, ulpinv
489 * ..
490 * .. Local Arrays ..
491  INTEGER iasign( maxtyp ), ibsign( maxtyp ),
492  $ ioldsd( 4 ), kadd( 6 ), kamagn( maxtyp ),
493  $ katype( maxtyp ), kazero( maxtyp ),
494  $ kbmagn( maxtyp ), kbtype( maxtyp ),
495  $ kbzero( maxtyp ), kclass( maxtyp ),
496  $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
497  REAL dumma( 4 ), rmagn( 0: 3 )
498 * ..
499 * .. External Functions ..
500  INTEGER ilaenv
501  REAL slamch, slarnd
502  EXTERNAL ilaenv, slamch, slarnd
503 * ..
504 * .. External Subroutines ..
505  EXTERNAL alasvm, sgegs, sgegv, sget51, sget52, sget53,
507  $ xerbla
508 * ..
509 * .. Intrinsic Functions ..
510  INTRINSIC abs, max, min, REAL, sign
511 * ..
512 * .. Data statements ..
513  DATA kclass / 15*1, 10*2, 1*3 /
514  DATA kz1 / 0, 1, 2, 1, 3, 3 /
515  DATA kz2 / 0, 0, 1, 2, 1, 1 /
516  DATA kadd / 0, 0, 0, 0, 3, 2 /
517  DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
518  $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
519  DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
520  $ 1, 1, -4, 2, -4, 8*8, 0 /
521  DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
522  $ 4*5, 4*3, 1 /
523  DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
524  $ 4*6, 4*4, 1 /
525  DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
526  $ 2, 1 /
527  DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
528  $ 2, 1 /
529  DATA ktrian / 16*0, 10*1 /
530  DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
531  $ 5*2, 0 /
532  DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
533 * ..
534 * .. Executable Statements ..
535 *
536 * Check for errors
537 *
538  info = 0
539 *
540  badnn = .false.
541  nmax = 1
542  DO 10 j = 1, nsizes
543  nmax = max( nmax, nn( j ) )
544  IF( nn( j ).LT.0 )
545  $ badnn = .true.
546  10 CONTINUE
547 *
548 * Maximum blocksize and shift -- we assume that blocksize and number
549 * of shifts are monotone increasing functions of N.
550 *
551  nb = max( 1, ilaenv( 1, 'SGEQRF', ' ', nmax, nmax, -1, -1 ),
552  $ ilaenv( 1, 'SORMQR', 'LT', nmax, nmax, nmax, -1 ),
553  $ ilaenv( 1, 'SORGQR', ' ', nmax, nmax, nmax, -1 ) )
554  nbz = ilaenv( 1, 'SHGEQZ', 'SII', nmax, 1, nmax, 0 )
555  ns = ilaenv( 4, 'SHGEQZ', 'SII', nmax, 1, nmax, 0 )
556  i1 = nbz + ns
557  lwkopt = 2*nmax + max( 6*nmax, nmax*( nb+1 ),
558  $ ( 2*i1+nmax+1 )*( i1+1 ) )
559 *
560 * Check for errors
561 *
562  IF( nsizes.LT.0 ) THEN
563  info = -1
564  ELSE IF( badnn ) THEN
565  info = -2
566  ELSE IF( ntypes.LT.0 ) THEN
567  info = -3
568  ELSE IF( thresh.LT.zero ) THEN
569  info = -6
570  ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
571  info = -10
572  ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
573  info = -19
574  ELSE IF( lwkopt.GT.lwork ) THEN
575  info = -30
576  END IF
577 *
578  IF( info.NE.0 ) THEN
579  CALL xerbla( 'SDRVGG', -info )
580  RETURN
581  END IF
582 *
583 * Quick return if possible
584 *
585  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
586  $ RETURN
587 *
588  safmin = slamch( 'Safe minimum' )
589  ulp = slamch( 'Epsilon' )*slamch( 'Base' )
590  safmin = safmin / ulp
591  safmax = one / safmin
592  CALL slabad( safmin, safmax )
593  ulpinv = one / ulp
594 *
595 * The values RMAGN(2:3) depend on N, see below.
596 *
597  rmagn( 0 ) = zero
598  rmagn( 1 ) = one
599 *
600 * Loop over sizes, types
601 *
602  ntestt = 0
603  nerrs = 0
604  nmats = 0
605 *
606  DO 170 jsize = 1, nsizes
607  n = nn( jsize )
608  n1 = max( 1, n )
609  rmagn( 2 ) = safmax*ulp / REAL( n1 )
610  rmagn( 3 ) = safmin*ulpinv*n1
611 *
612  IF( nsizes.NE.1 ) THEN
613  mtypes = min( maxtyp, ntypes )
614  ELSE
615  mtypes = min( maxtyp+1, ntypes )
616  END IF
617 *
618  DO 160 jtype = 1, mtypes
619  IF( .NOT.dotype( jtype ) )
620  $ go to 160
621  nmats = nmats + 1
622  ntest = 0
623 *
624 * Save ISEED in case of an error.
625 *
626  DO 20 j = 1, 4
627  ioldsd( j ) = iseed( j )
628  20 CONTINUE
629 *
630 * Initialize RESULT
631 *
632  DO 30 j = 1, 15
633  result( j ) = zero
634  30 CONTINUE
635 *
636 * Compute A and B
637 *
638 * Description of control parameters:
639 *
640 * KCLASS: =1 means w/o rotation, =2 means w/ rotation,
641 * =3 means random.
642 * KATYPE: the "type" to be passed to SLATM4 for computing A.
643 * KAZERO: the pattern of zeros on the diagonal for A:
644 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
645 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
646 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
647 * non-zero entries.)
648 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
649 * =2: large, =3: small.
650 * IASIGN: 1 if the diagonal elements of A are to be
651 * multiplied by a random magnitude 1 number, =2 if
652 * randomly chosen diagonal blocks are to be rotated
653 * to form 2x2 blocks.
654 * KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
655 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
656 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
657 * RMAGN: used to implement KAMAGN and KBMAGN.
658 *
659  IF( mtypes.GT.maxtyp )
660  $ go to 110
661  iinfo = 0
662  IF( kclass( jtype ).LT.3 ) THEN
663 *
664 * Generate A (w/o rotation)
665 *
666  IF( abs( katype( jtype ) ).EQ.3 ) THEN
667  in = 2*( ( n-1 ) / 2 ) + 1
668  IF( in.NE.n )
669  $ CALL slaset( 'Full', n, n, zero, zero, a, lda )
670  ELSE
671  in = n
672  END IF
673  CALL slatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
674  $ kz2( kazero( jtype ) ), iasign( jtype ),
675  $ rmagn( kamagn( jtype ) ), ulp,
676  $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
677  $ iseed, a, lda )
678  iadd = kadd( kazero( jtype ) )
679  IF( iadd.GT.0 .AND. iadd.LE.n )
680  $ a( iadd, iadd ) = one
681 *
682 * Generate B (w/o rotation)
683 *
684  IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
685  in = 2*( ( n-1 ) / 2 ) + 1
686  IF( in.NE.n )
687  $ CALL slaset( 'Full', n, n, zero, zero, b, lda )
688  ELSE
689  in = n
690  END IF
691  CALL slatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
692  $ kz2( kbzero( jtype ) ), ibsign( jtype ),
693  $ rmagn( kbmagn( jtype ) ), one,
694  $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
695  $ iseed, b, lda )
696  iadd = kadd( kbzero( jtype ) )
697  IF( iadd.NE.0 .AND. iadd.LE.n )
698  $ b( iadd, iadd ) = one
699 *
700  IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
701 *
702 * Include rotations
703 *
704 * Generate Q, Z as Householder transformations times
705 * a diagonal matrix.
706 *
707  DO 50 jc = 1, n - 1
708  DO 40 jr = jc, n
709  q( jr, jc ) = slarnd( 3, iseed )
710  z( jr, jc ) = slarnd( 3, iseed )
711  40 CONTINUE
712  CALL slarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
713  $ work( jc ) )
714  work( 2*n+jc ) = sign( one, q( jc, jc ) )
715  q( jc, jc ) = one
716  CALL slarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
717  $ work( n+jc ) )
718  work( 3*n+jc ) = sign( one, z( jc, jc ) )
719  z( jc, jc ) = one
720  50 CONTINUE
721  q( n, n ) = one
722  work( n ) = zero
723  work( 3*n ) = sign( one, slarnd( 2, iseed ) )
724  z( n, n ) = one
725  work( 2*n ) = zero
726  work( 4*n ) = sign( one, slarnd( 2, iseed ) )
727 *
728 * Apply the diagonal matrices
729 *
730  DO 70 jc = 1, n
731  DO 60 jr = 1, n
732  a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
733  $ a( jr, jc )
734  b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
735  $ b( jr, jc )
736  60 CONTINUE
737  70 CONTINUE
738  CALL sorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
739  $ lda, work( 2*n+1 ), iinfo )
740  IF( iinfo.NE.0 )
741  $ go to 100
742  CALL sorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
743  $ a, lda, work( 2*n+1 ), iinfo )
744  IF( iinfo.NE.0 )
745  $ go to 100
746  CALL sorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
747  $ lda, work( 2*n+1 ), iinfo )
748  IF( iinfo.NE.0 )
749  $ go to 100
750  CALL sorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
751  $ b, lda, work( 2*n+1 ), iinfo )
752  IF( iinfo.NE.0 )
753  $ go to 100
754  END IF
755  ELSE
756 *
757 * Random matrices
758 *
759  DO 90 jc = 1, n
760  DO 80 jr = 1, n
761  a( jr, jc ) = rmagn( kamagn( jtype ) )*
762  $ slarnd( 2, iseed )
763  b( jr, jc ) = rmagn( kbmagn( jtype ) )*
764  $ slarnd( 2, iseed )
765  80 CONTINUE
766  90 CONTINUE
767  END IF
768 *
769  100 CONTINUE
770 *
771  IF( iinfo.NE.0 ) THEN
772  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
773  $ ioldsd
774  info = abs( iinfo )
775  RETURN
776  END IF
777 *
778  110 CONTINUE
779 *
780 * Call SGEGS to compute H, T, Q, Z, alpha, and beta.
781 *
782  CALL slacpy( ' ', n, n, a, lda, s, lda )
783  CALL slacpy( ' ', n, n, b, lda, t, lda )
784  ntest = 1
785  result( 1 ) = ulpinv
786 *
787  CALL sgegs( 'V', 'V', n, s, lda, t, lda, alphr1, alphi1,
788  $ beta1, q, ldq, z, ldq, work, lwork, iinfo )
789  IF( iinfo.NE.0 ) THEN
790  WRITE( nounit, fmt = 9999 )'SGEGS', iinfo, n, jtype,
791  $ ioldsd
792  info = abs( iinfo )
793  go to 140
794  END IF
795 *
796  ntest = 4
797 *
798 * Do tests 1--4
799 *
800  CALL sget51( 1, n, a, lda, s, lda, q, ldq, z, ldq, work,
801  $ result( 1 ) )
802  CALL sget51( 1, n, b, lda, t, lda, q, ldq, z, ldq, work,
803  $ result( 2 ) )
804  CALL sget51( 3, n, b, lda, t, lda, q, ldq, q, ldq, work,
805  $ result( 3 ) )
806  CALL sget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
807  $ result( 4 ) )
808 *
809 * Do test 5: compare eigenvalues with diagonals.
810 * Also check Schur form of A.
811 *
812  temp1 = zero
813 *
814  DO 120 j = 1, n
815  ilabad = .false.
816  IF( alphi1( j ).EQ.zero ) THEN
817  temp2 = ( abs( alphr1( j )-s( j, j ) ) /
818  $ max( safmin, abs( alphr1( j ) ), abs( s( j,
819  $ j ) ) )+abs( beta1( j )-t( j, j ) ) /
820  $ max( safmin, abs( beta1( j ) ), abs( t( j,
821  $ j ) ) ) ) / ulp
822  IF( j.LT.n ) THEN
823  IF( s( j+1, j ).NE.zero )
824  $ ilabad = .true.
825  END IF
826  IF( j.GT.1 ) THEN
827  IF( s( j, j-1 ).NE.zero )
828  $ ilabad = .true.
829  END IF
830  ELSE
831  IF( alphi1( j ).GT.zero ) THEN
832  i1 = j
833  ELSE
834  i1 = j - 1
835  END IF
836  IF( i1.LE.0 .OR. i1.GE.n ) THEN
837  ilabad = .true.
838  ELSE IF( i1.LT.n-1 ) THEN
839  IF( s( i1+2, i1+1 ).NE.zero )
840  $ ilabad = .true.
841  ELSE IF( i1.GT.1 ) THEN
842  IF( s( i1, i1-1 ).NE.zero )
843  $ ilabad = .true.
844  END IF
845  IF( .NOT.ilabad ) THEN
846  CALL sget53( s( i1, i1 ), lda, t( i1, i1 ), lda,
847  $ beta1( j ), alphr1( j ), alphi1( j ),
848  $ temp2, iinfo )
849  IF( iinfo.GE.3 ) THEN
850  WRITE( nounit, fmt = 9997 )iinfo, j, n, jtype,
851  $ ioldsd
852  info = abs( iinfo )
853  END IF
854  ELSE
855  temp2 = ulpinv
856  END IF
857  END IF
858  temp1 = max( temp1, temp2 )
859  IF( ilabad ) THEN
860  WRITE( nounit, fmt = 9996 )j, n, jtype, ioldsd
861  END IF
862  120 CONTINUE
863  result( 5 ) = temp1
864 *
865 * Call SGEGV to compute S2, T2, VL, and VR, do tests.
866 *
867 * Eigenvalues and Eigenvectors
868 *
869  CALL slacpy( ' ', n, n, a, lda, s2, lda )
870  CALL slacpy( ' ', n, n, b, lda, t2, lda )
871  ntest = 6
872  result( 6 ) = ulpinv
873 *
874  CALL sgegv( 'V', 'V', n, s2, lda, t2, lda, alphr2, alphi2,
875  $ beta2, vl, ldq, vr, ldq, work, lwork, iinfo )
876  IF( iinfo.NE.0 ) THEN
877  WRITE( nounit, fmt = 9999 )'SGEGV', iinfo, n, jtype,
878  $ ioldsd
879  info = abs( iinfo )
880  go to 140
881  END IF
882 *
883  ntest = 7
884 *
885 * Do Tests 6 and 7
886 *
887  CALL sget52( .true., n, a, lda, b, lda, vl, ldq, alphr2,
888  $ alphi2, beta2, work, dumma( 1 ) )
889  result( 6 ) = dumma( 1 )
890  IF( dumma( 2 ).GT.thrshn ) THEN
891  WRITE( nounit, fmt = 9998 )'Left', 'SGEGV', dumma( 2 ),
892  $ n, jtype, ioldsd
893  END IF
894 *
895  CALL sget52( .false., n, a, lda, b, lda, vr, ldq, alphr2,
896  $ alphi2, beta2, work, dumma( 1 ) )
897  result( 7 ) = dumma( 1 )
898  IF( dumma( 2 ).GT.thresh ) THEN
899  WRITE( nounit, fmt = 9998 )'Right', 'SGEGV', dumma( 2 ),
900  $ n, jtype, ioldsd
901  END IF
902 *
903 * Check form of Complex eigenvalues.
904 *
905  DO 130 j = 1, n
906  ilabad = .false.
907  IF( alphi2( j ).GT.zero ) THEN
908  IF( j.EQ.n ) THEN
909  ilabad = .true.
910  ELSE IF( alphi2( j+1 ).GE.zero ) THEN
911  ilabad = .true.
912  END IF
913  ELSE IF( alphi2( j ).LT.zero ) THEN
914  IF( j.EQ.1 ) THEN
915  ilabad = .true.
916  ELSE IF( alphi2( j-1 ).LE.zero ) THEN
917  ilabad = .true.
918  END IF
919  END IF
920  IF( ilabad ) THEN
921  WRITE( nounit, fmt = 9996 )j, n, jtype, ioldsd
922  END IF
923  130 CONTINUE
924 *
925 * End of Loop -- Check for RESULT(j) > THRESH
926 *
927  140 CONTINUE
928 *
929  ntestt = ntestt + ntest
930 *
931 * Print out tests which fail.
932 *
933  DO 150 jr = 1, ntest
934  IF( result( jr ).GE.thresh ) THEN
935 *
936 * If this is the first test to fail,
937 * print a header to the data file.
938 *
939  IF( nerrs.EQ.0 ) THEN
940  WRITE( nounit, fmt = 9995 )'SGG'
941 *
942 * Matrix types
943 *
944  WRITE( nounit, fmt = 9994 )
945  WRITE( nounit, fmt = 9993 )
946  WRITE( nounit, fmt = 9992 )'Orthogonal'
947 *
948 * Tests performed
949 *
950  WRITE( nounit, fmt = 9991 )'orthogonal', '''',
951  $ 'transpose', ( '''', j = 1, 5 )
952 *
953  END IF
954  nerrs = nerrs + 1
955  IF( result( jr ).LT.10000.0 ) THEN
956  WRITE( nounit, fmt = 9990 )n, jtype, ioldsd, jr,
957  $ result( jr )
958  ELSE
959  WRITE( nounit, fmt = 9989 )n, jtype, ioldsd, jr,
960  $ result( jr )
961  END IF
962  END IF
963  150 CONTINUE
964 *
965  160 CONTINUE
966  170 CONTINUE
967 *
968 * Summary
969 *
970  CALL alasvm( 'SGG', nounit, nerrs, ntestt, 0 )
971  RETURN
972 *
973  9999 FORMAT( ' SDRVGG: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
974  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
975 *
976  9998 FORMAT( ' SDRVGG: ', a, ' Eigenvectors from ', a, ' incorrectly ',
977  $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
978  $ 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5,
979  $ ')' )
980 *
981  9997 FORMAT( ' SDRVGG: SGET53 returned INFO=', i1, ' for eigenvalue ',
982  $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(',
983  $ 3( i5, ',' ), i5, ')' )
984 *
985  9996 FORMAT( ' SDRVGG: S not in Schur form at eigenvalue ', i6, '.',
986  $ / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
987  $ i5, ')' )
988 *
989  9995 FORMAT( / 1x, a3, ' -- Real Generalized eigenvalue problem driver'
990  $ )
991 *
992  9994 FORMAT( ' Matrix types (see SDRVGG for details): ' )
993 *
994  9993 FORMAT( ' Special Matrices:', 23x,
995  $ '(J''=transposed Jordan block)',
996  $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
997  $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
998  $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
999  $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
1000  $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
1001  $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
1002  9992 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
1003  $ / ' 16=Transposed Jordan Blocks 19=geometric ',
1004  $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
1005  $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
1006  $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
1007  $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
1008  $ '23=(small,large) 24=(small,small) 25=(large,large)',
1009  $ / ' 26=random O(1) matrices.' )
1010 *
1011  9991 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
1012  $ 'Q and Z are ', a, ',', / 20x,
1013  $ 'l and r are the appropriate left and right', / 19x,
1014  $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
1015  $ ' means ', a, '.)', / ' 1 = | A - Q S Z', a,
1016  $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
1017  $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
1018  $ ' | / ( n ulp ) 4 = | I - ZZ', a,
1019  $ ' | / ( n ulp )', /
1020  $ ' 5 = difference between (alpha,beta) and diagonals of',
1021  $ ' (S,T)', / ' 6 = max | ( b A - a B )', a,
1022  $ ' l | / const. 7 = max | ( b A - a B ) r | / const.',
1023  $ / 1x )
1024  9990 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
1025  $ 4( i4, ',' ), ' result ', i3, ' is', 0p, f8.2 )
1026  9989 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
1027  $ 4( i4, ',' ), ' result ', i3, ' is', 1p, e10.3 )
1028 *
1029 * End of SDRVGG
1030 *
1031  END
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:111
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:74
subroutine sgegv(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
SGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
Definition: sgegv.f:306
subroutine sget53(A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO)
SGET53
Definition: sget53.f:127
subroutine sget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
SGET52
Definition: sget52.f:199
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:104
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
real function slarnd(IDIST, ISEED)
SLARND
Definition: slarnd.f:74
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:107
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
Definition: xerbla-fortran:9
subroutine slatm4(ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
SLATM4
Definition: slatm4.f:175
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:75
subroutine sgegs(JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, INFO)
SGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
Definition: sgegs.f:226
subroutine sget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
SGET51
Definition: sget51.f:149
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:81
subroutine sdrvgg(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, THRSHN, NOUNIT, A, LDA, B, S, T, S2, T2, Q, LDQ, Z, ALPHR1, ALPHI1, BETA1, ALPHR2, ALPHI2, BETA2, VL, VR, WORK, LWORK, RESULT, INFO)
SDRVGG
Definition: sdrvgg.f:450