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