LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 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 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 threshold 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 *> \ingroup double_eig
402 *
403 * =====================================================================
404  SUBROUTINE ddrgev( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
405  $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE,
406  $ ALPHAR, ALPHAI, BETA, ALPHR1, ALPHI1, BETA1,
407  $ WORK, LWORK, RESULT, INFO )
408 *
409 * -- LAPACK test routine --
410 * -- LAPACK is a software package provided by Univ. of Tennessee, --
411 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
412 *
413 * .. Scalar Arguments ..
414  INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
415  $ NTYPES
416  DOUBLE PRECISION THRESH
417 * ..
418 * .. Array Arguments ..
419  LOGICAL DOTYPE( * )
420  INTEGER ISEED( 4 ), NN( * )
421  DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
422  $ alphi1( * ), alphr1( * ), b( lda, * ),
423  $ beta( * ), beta1( * ), q( ldq, * ),
424  $ qe( ldqe, * ), result( * ), s( lda, * ),
425  $ t( lda, * ), work( * ), z( ldq, * )
426 * ..
427 *
428 * =====================================================================
429 *
430 * .. Parameters ..
431  DOUBLE PRECISION ZERO, ONE
432  PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
433  INTEGER MAXTYP
434  parameter( maxtyp = 26 )
435 * ..
436 * .. Local Scalars ..
437  LOGICAL BADNN
438  INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
439  $ MAXWRK, MINWRK, MTYPES, N, N1, NERRS, NMATS,
440  $ nmax, ntestt
441  DOUBLE PRECISION SAFMAX, SAFMIN, ULP, ULPINV
442 * ..
443 * .. Local Arrays ..
444  INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
445  $ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
446  $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
447  $ kbmagn( maxtyp ), kbtype( maxtyp ),
448  $ kbzero( maxtyp ), kclass( maxtyp ),
449  $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
450  DOUBLE PRECISION RMAGN( 0: 3 )
451 * ..
452 * .. External Functions ..
453  INTEGER ILAENV
454  DOUBLE PRECISION DLAMCH, DLARND
455  EXTERNAL ILAENV, DLAMCH, DLARND
456 * ..
457 * .. External Subroutines ..
458  EXTERNAL alasvm, dget52, dggev, dlabad, dlacpy, dlarfg,
460 * ..
461 * .. Intrinsic Functions ..
462  INTRINSIC abs, dble, max, min, sign
463 * ..
464 * .. Data statements ..
465  DATA kclass / 15*1, 10*2, 1*3 /
466  DATA kz1 / 0, 1, 2, 1, 3, 3 /
467  DATA kz2 / 0, 0, 1, 2, 1, 1 /
468  DATA kadd / 0, 0, 0, 0, 3, 2 /
469  DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
470  $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
471  DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
472  $ 1, 1, -4, 2, -4, 8*8, 0 /
473  DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
474  $ 4*5, 4*3, 1 /
475  DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
476  $ 4*6, 4*4, 1 /
477  DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
478  $ 2, 1 /
479  DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
480  $ 2, 1 /
481  DATA ktrian / 16*0, 10*1 /
482  DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
483  $ 5*2, 0 /
484  DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
485 * ..
486 * .. Executable Statements ..
487 *
488 * Check for errors
489 *
490  info = 0
491 *
492  badnn = .false.
493  nmax = 1
494  DO 10 j = 1, nsizes
495  nmax = max( nmax, nn( j ) )
496  IF( nn( j ).LT.0 )
497  $ badnn = .true.
498  10 CONTINUE
499 *
500  IF( nsizes.LT.0 ) THEN
501  info = -1
502  ELSE IF( badnn ) THEN
503  info = -2
504  ELSE IF( ntypes.LT.0 ) THEN
505  info = -3
506  ELSE IF( thresh.LT.zero ) THEN
507  info = -6
508  ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
509  info = -9
510  ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
511  info = -14
512  ELSE IF( ldqe.LE.1 .OR. ldqe.LT.nmax ) THEN
513  info = -17
514  END IF
515 *
516 * Compute workspace
517 * (Note: Comments in the code beginning "Workspace:" describe the
518 * minimal amount of workspace needed at that point in the code,
519 * as well as the preferred amount for good performance.
520 * NB refers to the optimal block size for the immediately
521 * following subroutine, as returned by ILAENV.
522 *
523  minwrk = 1
524  IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
525  minwrk = max( 1, 8*nmax, nmax*( nmax+1 ) )
526  maxwrk = 7*nmax + nmax*ilaenv( 1, 'DGEQRF', ' ', nmax, 1, nmax,
527  $ 0 )
528  maxwrk = max( maxwrk, nmax*( nmax+1 ) )
529  work( 1 ) = maxwrk
530  END IF
531 *
532  IF( lwork.LT.minwrk )
533  $ info = -25
534 *
535  IF( info.NE.0 ) THEN
536  CALL xerbla( 'DDRGEV', -info )
537  RETURN
538  END IF
539 *
540 * Quick return if possible
541 *
542  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
543  $ RETURN
544 *
545  safmin = dlamch( 'Safe minimum' )
546  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
547  safmin = safmin / ulp
548  safmax = one / safmin
549  CALL dlabad( safmin, safmax )
550  ulpinv = one / ulp
551 *
552 * The values RMAGN(2:3) depend on N, see below.
553 *
554  rmagn( 0 ) = zero
555  rmagn( 1 ) = one
556 *
557 * Loop over sizes, types
558 *
559  ntestt = 0
560  nerrs = 0
561  nmats = 0
562 *
563  DO 220 jsize = 1, nsizes
564  n = nn( jsize )
565  n1 = max( 1, n )
566  rmagn( 2 ) = safmax*ulp / dble( n1 )
567  rmagn( 3 ) = safmin*ulpinv*n1
568 *
569  IF( nsizes.NE.1 ) THEN
570  mtypes = min( maxtyp, ntypes )
571  ELSE
572  mtypes = min( maxtyp+1, ntypes )
573  END IF
574 *
575  DO 210 jtype = 1, mtypes
576  IF( .NOT.dotype( jtype ) )
577  $ GO TO 210
578  nmats = nmats + 1
579 *
580 * Save ISEED in case of an error.
581 *
582  DO 20 j = 1, 4
583  ioldsd( j ) = iseed( j )
584  20 CONTINUE
585 *
586 * Generate test matrices A and B
587 *
588 * Description of control parameters:
589 *
590 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
591 * =3 means random.
592 * KATYPE: the "type" to be passed to DLATM4 for computing A.
593 * KAZERO: the pattern of zeros on the diagonal for A:
594 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
595 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
596 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
597 * non-zero entries.)
598 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
599 * =2: large, =3: small.
600 * IASIGN: 1 if the diagonal elements of A are to be
601 * multiplied by a random magnitude 1 number, =2 if
602 * randomly chosen diagonal blocks are to be rotated
603 * to form 2x2 blocks.
604 * KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
605 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
606 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
607 * RMAGN: used to implement KAMAGN and KBMAGN.
608 *
609  IF( mtypes.GT.maxtyp )
610  $ GO TO 100
611  ierr = 0
612  IF( kclass( jtype ).LT.3 ) THEN
613 *
614 * Generate A (w/o rotation)
615 *
616  IF( abs( katype( jtype ) ).EQ.3 ) THEN
617  in = 2*( ( n-1 ) / 2 ) + 1
618  IF( in.NE.n )
619  $ CALL dlaset( 'Full', n, n, zero, zero, a, lda )
620  ELSE
621  in = n
622  END IF
623  CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
624  $ kz2( kazero( jtype ) ), iasign( jtype ),
625  $ rmagn( kamagn( jtype ) ), ulp,
626  $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
627  $ iseed, a, lda )
628  iadd = kadd( kazero( jtype ) )
629  IF( iadd.GT.0 .AND. iadd.LE.n )
630  $ a( iadd, iadd ) = one
631 *
632 * Generate B (w/o rotation)
633 *
634  IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
635  in = 2*( ( n-1 ) / 2 ) + 1
636  IF( in.NE.n )
637  $ CALL dlaset( 'Full', n, n, zero, zero, b, lda )
638  ELSE
639  in = n
640  END IF
641  CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
642  $ kz2( kbzero( jtype ) ), ibsign( jtype ),
643  $ rmagn( kbmagn( jtype ) ), one,
644  $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
645  $ iseed, b, lda )
646  iadd = kadd( kbzero( jtype ) )
647  IF( iadd.NE.0 .AND. iadd.LE.n )
648  $ b( iadd, iadd ) = one
649 *
650  IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
651 *
652 * Include rotations
653 *
654 * Generate Q, Z as Householder transformations times
655 * a diagonal matrix.
656 *
657  DO 40 jc = 1, n - 1
658  DO 30 jr = jc, n
659  q( jr, jc ) = dlarnd( 3, iseed )
660  z( jr, jc ) = dlarnd( 3, iseed )
661  30 CONTINUE
662  CALL dlarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
663  $ work( jc ) )
664  work( 2*n+jc ) = sign( one, q( jc, jc ) )
665  q( jc, jc ) = one
666  CALL dlarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
667  $ work( n+jc ) )
668  work( 3*n+jc ) = sign( one, z( jc, jc ) )
669  z( jc, jc ) = one
670  40 CONTINUE
671  q( n, n ) = one
672  work( n ) = zero
673  work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
674  z( n, n ) = one
675  work( 2*n ) = zero
676  work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
677 *
678 * Apply the diagonal matrices
679 *
680  DO 60 jc = 1, n
681  DO 50 jr = 1, n
682  a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
683  $ a( jr, jc )
684  b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
685  $ b( jr, jc )
686  50 CONTINUE
687  60 CONTINUE
688  CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
689  $ lda, work( 2*n+1 ), ierr )
690  IF( ierr.NE.0 )
691  $ GO TO 90
692  CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
693  $ a, lda, work( 2*n+1 ), ierr )
694  IF( ierr.NE.0 )
695  $ GO TO 90
696  CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
697  $ lda, work( 2*n+1 ), ierr )
698  IF( ierr.NE.0 )
699  $ GO TO 90
700  CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
701  $ b, lda, work( 2*n+1 ), ierr )
702  IF( ierr.NE.0 )
703  $ GO TO 90
704  END IF
705  ELSE
706 *
707 * Random matrices
708 *
709  DO 80 jc = 1, n
710  DO 70 jr = 1, n
711  a( jr, jc ) = rmagn( kamagn( jtype ) )*
712  $ dlarnd( 2, iseed )
713  b( jr, jc ) = rmagn( kbmagn( jtype ) )*
714  $ dlarnd( 2, iseed )
715  70 CONTINUE
716  80 CONTINUE
717  END IF
718 *
719  90 CONTINUE
720 *
721  IF( ierr.NE.0 ) THEN
722  WRITE( nounit, fmt = 9999 )'Generator', ierr, n, jtype,
723  $ ioldsd
724  info = abs( ierr )
725  RETURN
726  END IF
727 *
728  100 CONTINUE
729 *
730  DO 110 i = 1, 7
731  result( i ) = -one
732  110 CONTINUE
733 *
734 * Call DGGEV to compute eigenvalues and eigenvectors.
735 *
736  CALL dlacpy( ' ', n, n, a, lda, s, lda )
737  CALL dlacpy( ' ', n, n, b, lda, t, lda )
738  CALL dggev( 'V', 'V', n, s, lda, t, lda, alphar, alphai,
739  $ beta, q, ldq, z, ldq, work, lwork, ierr )
740  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
741  result( 1 ) = ulpinv
742  WRITE( nounit, fmt = 9999 )'DGGEV1', ierr, n, jtype,
743  $ ioldsd
744  info = abs( ierr )
745  GO TO 190
746  END IF
747 *
748 * Do the tests (1) and (2)
749 *
750  CALL dget52( .true., n, a, lda, b, lda, q, ldq, alphar,
751  $ alphai, beta, work, result( 1 ) )
752  IF( result( 2 ).GT.thresh ) THEN
753  WRITE( nounit, fmt = 9998 )'Left', 'DGGEV1',
754  $ result( 2 ), n, jtype, ioldsd
755  END IF
756 *
757 * Do the tests (3) and (4)
758 *
759  CALL dget52( .false., n, a, lda, b, lda, z, ldq, alphar,
760  $ alphai, beta, work, result( 3 ) )
761  IF( result( 4 ).GT.thresh ) THEN
762  WRITE( nounit, fmt = 9998 )'Right', 'DGGEV1',
763  $ result( 4 ), n, jtype, ioldsd
764  END IF
765 *
766 * Do the test (5)
767 *
768  CALL dlacpy( ' ', n, n, a, lda, s, lda )
769  CALL dlacpy( ' ', n, n, b, lda, t, lda )
770  CALL dggev( 'N', 'N', n, s, lda, t, lda, alphr1, alphi1,
771  $ beta1, q, ldq, z, ldq, work, lwork, ierr )
772  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
773  result( 1 ) = ulpinv
774  WRITE( nounit, fmt = 9999 )'DGGEV2', ierr, n, jtype,
775  $ ioldsd
776  info = abs( ierr )
777  GO TO 190
778  END IF
779 *
780  DO 120 j = 1, n
781  IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
782  $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 5 )
783  $ = ulpinv
784  120 CONTINUE
785 *
786 * Do the test (6): Compute eigenvalues and left eigenvectors,
787 * and test them
788 *
789  CALL dlacpy( ' ', n, n, a, lda, s, lda )
790  CALL dlacpy( ' ', n, n, b, lda, t, lda )
791  CALL dggev( 'V', 'N', n, s, lda, t, lda, alphr1, alphi1,
792  $ beta1, qe, ldqe, z, ldq, work, lwork, ierr )
793  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
794  result( 1 ) = ulpinv
795  WRITE( nounit, fmt = 9999 )'DGGEV3', ierr, n, jtype,
796  $ ioldsd
797  info = abs( ierr )
798  GO TO 190
799  END IF
800 *
801  DO 130 j = 1, n
802  IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
803  $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 6 )
804  $ = ulpinv
805  130 CONTINUE
806 *
807  DO 150 j = 1, n
808  DO 140 jc = 1, n
809  IF( q( j, jc ).NE.qe( j, jc ) )
810  $ result( 6 ) = ulpinv
811  140 CONTINUE
812  150 CONTINUE
813 *
814 * DO the test (7): Compute eigenvalues and right eigenvectors,
815 * and test them
816 *
817  CALL dlacpy( ' ', n, n, a, lda, s, lda )
818  CALL dlacpy( ' ', n, n, b, lda, t, lda )
819  CALL dggev( 'N', 'V', n, s, lda, t, lda, alphr1, alphi1,
820  $ beta1, q, ldq, qe, ldqe, work, lwork, ierr )
821  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
822  result( 1 ) = ulpinv
823  WRITE( nounit, fmt = 9999 )'DGGEV4', ierr, n, jtype,
824  $ ioldsd
825  info = abs( ierr )
826  GO TO 190
827  END IF
828 *
829  DO 160 j = 1, n
830  IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
831  $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 7 )
832  $ = ulpinv
833  160 CONTINUE
834 *
835  DO 180 j = 1, n
836  DO 170 jc = 1, n
837  IF( z( j, jc ).NE.qe( j, jc ) )
838  $ result( 7 ) = ulpinv
839  170 CONTINUE
840  180 CONTINUE
841 *
842 * End of Loop -- Check for RESULT(j) > THRESH
843 *
844  190 CONTINUE
845 *
846  ntestt = ntestt + 7
847 *
848 * Print out tests which fail.
849 *
850  DO 200 jr = 1, 7
851  IF( result( jr ).GE.thresh ) THEN
852 *
853 * If this is the first test to fail,
854 * print a header to the data file.
855 *
856  IF( nerrs.EQ.0 ) THEN
857  WRITE( nounit, fmt = 9997 )'DGV'
858 *
859 * Matrix types
860 *
861  WRITE( nounit, fmt = 9996 )
862  WRITE( nounit, fmt = 9995 )
863  WRITE( nounit, fmt = 9994 )'Orthogonal'
864 *
865 * Tests performed
866 *
867  WRITE( nounit, fmt = 9993 )
868 *
869  END IF
870  nerrs = nerrs + 1
871  IF( result( jr ).LT.10000.0d0 ) THEN
872  WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
873  $ result( jr )
874  ELSE
875  WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
876  $ result( jr )
877  END IF
878  END IF
879  200 CONTINUE
880 *
881  210 CONTINUE
882  220 CONTINUE
883 *
884 * Summary
885 *
886  CALL alasvm( 'DGV', nounit, nerrs, ntestt, 0 )
887 *
888  work( 1 ) = maxwrk
889 *
890  RETURN
891 *
892  9999 FORMAT( ' DDRGEV: ', a, ' returned INFO=', i6, '.', / 3x, 'N=',
893  $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
894 *
895  9998 FORMAT( ' DDRGEV: ', a, ' Eigenvectors from ', a, ' incorrectly ',
896  $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 3x,
897  $ 'N=', i4, ', JTYPE=', i3, ', ISEED=(', 4( i4, ',' ), i5,
898  $ ')' )
899 *
900  9997 FORMAT( / 1x, a3, ' -- Real Generalized eigenvalue problem driver'
901  $ )
902 *
903  9996 FORMAT( ' Matrix types (see DDRGEV for details): ' )
904 *
905  9995 FORMAT( ' Special Matrices:', 23x,
906  $ '(J''=transposed Jordan block)',
907  $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
908  $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
909  $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
910  $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
911  $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
912  $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
913  9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
914  $ / ' 16=Transposed Jordan Blocks 19=geometric ',
915  $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
916  $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
917  $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
918  $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
919  $ '23=(small,large) 24=(small,small) 25=(large,large)',
920  $ / ' 26=random O(1) matrices.' )
921 *
922  9993 FORMAT( / ' Tests performed: ',
923  $ / ' 1 = max | ( b A - a B )''*l | / const.,',
924  $ / ' 2 = | |VR(i)| - 1 | / ulp,',
925  $ / ' 3 = max | ( b A - a B )*r | / const.',
926  $ / ' 4 = | |VL(i)| - 1 | / ulp,',
927  $ / ' 5 = 0 if W same no matter if r or l computed,',
928  $ / ' 6 = 0 if l same no matter if l computed,',
929  $ / ' 7 = 0 if r same no matter if r computed,', / 1x )
930  9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
931  $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
932  9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
933  $ 4( i4, ',' ), ' result ', i2, ' is', 1p, d10.3 )
934 *
935 * End of DDRGEV
936 *
937  END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:73
subroutine dlatm4(ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
DLATM4
Definition: dlatm4.f:175
subroutine ddrgev(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)
DDRGEV
Definition: ddrgev.f:408
subroutine dget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
DGET52
Definition: dget52.f:199
subroutine dggev(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
DGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition: dggev.f:226
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:106
subroutine dorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: dorm2r.f:159