LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
ddrgev.f
Go to the documentation of this file.
1 *> \brief \b DDRGEV
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 DDRGEV( 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 * DOUBLE PRECISION THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL DOTYPE( * )
23 * INTEGER ISEED( 4 ), NN( * )
24 * DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
25 * $ ALPHI1( * ), 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 *> DDRGEV checks the nonsymmetric generalized eigenvalue problem driver
38 *> routine DGGEV.
39 *>
40 *> DGGEV 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 reasonalbe
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 DDRGEV 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 threshhold THRESH.
59 *>
60 *> Results from DGGEV:
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 *> DDRGES 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, DDRGES
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 DDRGES to continue the same random number
233 *> sequence.
234 *> \endverbatim
235 *>
236 *> \param[in] THRESH
237 *> \verbatim
238 *> THRESH is DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array,
282 *> dimension (LDA, max(NN))
283 *> The Schur form matrix computed from A by DGGES. 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 DOUBLE PRECISION array,
291 *> dimension (LDA, max(NN))
292 *> The upper triangular matrix computed from B by DGGES.
293 *> \endverbatim
294 *>
295 *> \param[out] Q
296 *> \verbatim
297 *> Q is DOUBLE PRECISION array,
298 *> dimension (LDQ, max(NN))
299 *> The (left) eigenvectors matrix computed by DGGEV.
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 DOUBLE PRECISION array, dimension( LDQ, max(NN) )
312 *> The (right) orthogonal matrix computed by DGGES.
313 *> \endverbatim
314 *>
315 *> \param[out] QE
316 *> \verbatim
317 *> QE is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (max(NN))
330 *> \endverbatim
331 *>
332 *> \param[out] ALPHAI
333 *> \verbatim
334 *> ALPHAI is DOUBLE PRECISION array, dimension (max(NN))
335 *> \endverbatim
336 *>
337 *> \param[out] BETA
338 *> \verbatim
339 *> BETA is DOUBLE PRECISION array, dimension (max(NN))
340 *>
341 *> The generalized eigenvalues of (A,B) computed by DGGEV.
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 DOUBLE PRECISION array, dimension (max(NN))
349 *> \endverbatim
350 *>
351 *> \param[out] ALPHI1
352 *> \verbatim
353 *> ALPHI1 is DOUBLE PRECISION array, dimension (max(NN))
354 *> \endverbatim
355 *>
356 *> \param[out] BETA1
357 *> \verbatim
358 *> BETA1 is DOUBLE PRECISION array, dimension (max(NN))
359 *>
360 *> Like ALPHAR, ALPHAI, BETA, these arrays contain the
361 *> eigenvalues of A and B, but those computed when DGGEV 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 November 2011
402 *
403 *> \ingroup double_eig
404 *
405 * =====================================================================
406  SUBROUTINE ddrgev( 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.4.0) --
412 * -- LAPACK is a software package provided by Univ. of Tennessee, --
413 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
414 * November 2011
415 *
416 * .. Scalar Arguments ..
417  INTEGER info, lda, ldq, ldqe, lwork, nounit, nsizes,
418  $ ntypes
419  DOUBLE PRECISION thresh
420 * ..
421 * .. Array Arguments ..
422  LOGICAL dotype( * )
423  INTEGER iseed( 4 ), nn( * )
424  DOUBLE PRECISION a( lda, * ), alphai( * ), alphar( * ),
425  $ alphi1( * ), 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  DOUBLE PRECISION zero, one
435  parameter( zero = 0.0d+0, one = 1.0d+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  DOUBLE PRECISION 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  DOUBLE PRECISION rmagn( 0: 3 )
454 * ..
455 * .. External Functions ..
456  INTEGER ilaenv
457  DOUBLE PRECISION dlamch, dlarnd
458  EXTERNAL ilaenv, dlamch, dlarnd
459 * ..
460 * .. External Subroutines ..
461  EXTERNAL alasvm, dget52, dggev, dlabad, dlacpy, dlarfg,
463 * ..
464 * .. Intrinsic Functions ..
465  INTRINSIC abs, dble, max, min, 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, 'DGEQRF', ' ', 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( 'DDRGEV', -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 = dlamch( 'Safe minimum' )
549  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
550  safmin = safmin / ulp
551  safmax = one / safmin
552  CALL dlabad( 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 / dble( 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 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
594 * =3 means random.
595 * KATYPE: the "type" to be passed to DLATM4 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 dlaset( 'Full', n, n, zero, zero, a, lda )
623  ELSE
624  in = n
625  END IF
626  CALL dlatm4( 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 dlaset( 'Full', n, n, zero, zero, b, lda )
641  ELSE
642  in = n
643  END IF
644  CALL dlatm4( 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 ) = dlarnd( 3, iseed )
663  z( jr, jc ) = dlarnd( 3, iseed )
664  30 CONTINUE
665  CALL dlarfg( 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 dlarfg( 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, dlarnd( 2, iseed ) )
677  z( n, n ) = one
678  work( 2*n ) = zero
679  work( 4*n ) = sign( one, dlarnd( 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 dorm2r( '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 dorm2r( '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 dorm2r( '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 dorm2r( '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  $ dlarnd( 2, iseed )
716  b( jr, jc ) = rmagn( kbmagn( jtype ) )*
717  $ dlarnd( 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 DGGEV to compute eigenvalues and eigenvectors.
738 *
739  CALL dlacpy( ' ', n, n, a, lda, s, lda )
740  CALL dlacpy( ' ', n, n, b, lda, t, lda )
741  CALL dggev( '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 )'DGGEV1', 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 dget52( .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', 'DGGEV1',
757  $ result( 2 ), n, jtype, ioldsd
758  END IF
759 *
760 * Do the tests (3) and (4)
761 *
762  CALL dget52( .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', 'DGGEV1',
766  $ result( 4 ), n, jtype, ioldsd
767  END IF
768 *
769 * Do the test (5)
770 *
771  CALL dlacpy( ' ', n, n, a, lda, s, lda )
772  CALL dlacpy( ' ', n, n, b, lda, t, lda )
773  CALL dggev( '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 )'DGGEV2', 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. alphai( j ).NE.
785  $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 5 )
786  $ = ulpinv
787  120 CONTINUE
788 *
789 * Do the test (6): Compute eigenvalues and left eigenvectors,
790 * and test them
791 *
792  CALL dlacpy( ' ', n, n, a, lda, s, lda )
793  CALL dlacpy( ' ', n, n, b, lda, t, lda )
794  CALL dggev( 'V', 'N', n, s, lda, t, lda, alphr1, alphi1,
795  $ beta1, qe, ldqe, z, ldq, work, lwork, ierr )
796  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
797  result( 1 ) = ulpinv
798  WRITE( nounit, fmt = 9999 )'DGGEV3', ierr, n, jtype,
799  $ ioldsd
800  info = abs( ierr )
801  go to 190
802  END IF
803 *
804  DO 130 j = 1, n
805  IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
806  $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 6 )
807  $ = ulpinv
808  130 CONTINUE
809 *
810  DO 150 j = 1, n
811  DO 140 jc = 1, n
812  IF( q( j, jc ).NE.qe( j, jc ) )
813  $ result( 6 ) = ulpinv
814  140 CONTINUE
815  150 CONTINUE
816 *
817 * DO the test (7): Compute eigenvalues and right eigenvectors,
818 * and test them
819 *
820  CALL dlacpy( ' ', n, n, a, lda, s, lda )
821  CALL dlacpy( ' ', n, n, b, lda, t, lda )
822  CALL dggev( 'N', 'V', n, s, lda, t, lda, alphr1, alphi1,
823  $ beta1, q, ldq, qe, ldqe, work, lwork, ierr )
824  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
825  result( 1 ) = ulpinv
826  WRITE( nounit, fmt = 9999 )'DGGEV4', ierr, n, jtype,
827  $ ioldsd
828  info = abs( ierr )
829  go to 190
830  END IF
831 *
832  DO 160 j = 1, n
833  IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
834  $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 7 )
835  $ = ulpinv
836  160 CONTINUE
837 *
838  DO 180 j = 1, n
839  DO 170 jc = 1, n
840  IF( z( j, jc ).NE.qe( j, jc ) )
841  $ result( 7 ) = ulpinv
842  170 CONTINUE
843  180 CONTINUE
844 *
845 * End of Loop -- Check for RESULT(j) > THRESH
846 *
847  190 CONTINUE
848 *
849  ntestt = ntestt + 7
850 *
851 * Print out tests which fail.
852 *
853  DO 200 jr = 1, 7
854  IF( result( jr ).GE.thresh ) THEN
855 *
856 * If this is the first test to fail,
857 * print a header to the data file.
858 *
859  IF( nerrs.EQ.0 ) THEN
860  WRITE( nounit, fmt = 9997 )'DGV'
861 *
862 * Matrix types
863 *
864  WRITE( nounit, fmt = 9996 )
865  WRITE( nounit, fmt = 9995 )
866  WRITE( nounit, fmt = 9994 )'Orthogonal'
867 *
868 * Tests performed
869 *
870  WRITE( nounit, fmt = 9993 )
871 *
872  END IF
873  nerrs = nerrs + 1
874  IF( result( jr ).LT.10000.0d0 ) THEN
875  WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
876  $ result( jr )
877  ELSE
878  WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
879  $ result( jr )
880  END IF
881  END IF
882  200 CONTINUE
883 *
884  210 CONTINUE
885  220 CONTINUE
886 *
887 * Summary
888 *
889  CALL alasvm( 'DGV', nounit, nerrs, ntestt, 0 )
890 *
891  work( 1 ) = maxwrk
892 *
893  RETURN
894 *
895  9999 FORMAT( ' DDRGEV: ', a, ' returned INFO=', i6, '.', / 3x, 'N=',
896  $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
897 *
898  9998 FORMAT( ' DDRGEV: ', a, ' Eigenvectors from ', a, ' incorrectly ',
899  $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 3x,
900  $ 'N=', i4, ', JTYPE=', i3, ', ISEED=(', 4( i4, ',' ), i5,
901  $ ')' )
902 *
903  9997 FORMAT( / 1x, a3, ' -- Real Generalized eigenvalue problem driver'
904  $ )
905 *
906  9996 FORMAT( ' Matrix types (see DDRGEV for details): ' )
907 *
908  9995 FORMAT( ' Special Matrices:', 23x,
909  $ '(J''=transposed Jordan block)',
910  $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
911  $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
912  $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
913  $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
914  $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
915  $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
916  9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
917  $ / ' 16=Transposed Jordan Blocks 19=geometric ',
918  $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
919  $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
920  $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
921  $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
922  $ '23=(small,large) 24=(small,small) 25=(large,large)',
923  $ / ' 26=random O(1) matrices.' )
924 *
925  9993 FORMAT( / ' Tests performed: ',
926  $ / ' 1 = max | ( b A - a B )''*l | / const.,',
927  $ / ' 2 = | |VR(i)| - 1 | / ulp,',
928  $ / ' 3 = max | ( b A - a B )*r | / const.',
929  $ / ' 4 = | |VL(i)| - 1 | / ulp,',
930  $ / ' 5 = 0 if W same no matter if r or l computed,',
931  $ / ' 6 = 0 if l same no matter if l computed,',
932  $ / ' 7 = 0 if r same no matter if r computed,', / 1x )
933  9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
934  $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
935  9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
936  $ 4( i4, ',' ), ' result ', i2, ' is', 1p, d10.3 )
937 *
938 * End of DDRGEV
939 *
940  END