LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ddrgev3.f
Go to the documentation of this file.
1*> \brief \b DDRGEV3
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 DDRGEV3( 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*> DDRGEV3 checks the nonsymmetric generalized eigenvalue problem driver
38*> routine DGGEV3.
39*>
40*> DGGEV3 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 DDRGEV3 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 DGGEV3:
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*> DDRGEV3 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, DDRGEV3
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 DDRGEV3 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 DGGEV3. 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 DGGEV3.
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 DGGEV3.
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 DGGEV3.
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 DGGEV3.
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 DGGEV3 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 ddrgev3( 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 = 27 )
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, dggev3, 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, 1*4 /
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, 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, 0 /
473 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
474 $ 4*5, 4*3, 1, 1 /
475 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
476 $ 4*6, 4*4, 1, 1 /
477 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
478 $ 2, 1, 3 /
479 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
480 $ 2, 1, 3 /
481 DATA ktrian / 16*0, 11*1 /
482 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
483 $ 5*2, 2*0 /
484 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 10*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( 'DDRGEV3', -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, =4 means random generalized
592* upper Hessenberg.
593* KATYPE: the "type" to be passed to DLATM4 for computing A.
594* KAZERO: the pattern of zeros on the diagonal for A:
595* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
596* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
597* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
598* non-zero entries.)
599* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
600* =2: large, =3: small.
601* IASIGN: 1 if the diagonal elements of A are to be
602* multiplied by a random magnitude 1 number, =2 if
603* randomly chosen diagonal blocks are to be rotated
604* to form 2x2 blocks.
605* KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
606* KTRIAN: =0: don't fill in the upper triangle, =1: do.
607* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
608* RMAGN: used to implement KAMAGN and KBMAGN.
609*
610 IF( mtypes.GT.maxtyp )
611 $ GO TO 100
612 ierr = 0
613 IF( kclass( jtype ).LT.3 ) THEN
614*
615* Generate A (w/o rotation)
616*
617 IF( abs( katype( jtype ) ).EQ.3 ) THEN
618 in = 2*( ( n-1 ) / 2 ) + 1
619 IF( in.NE.n )
620 $ CALL dlaset( 'Full', n, n, zero, zero, a, lda )
621 ELSE
622 in = n
623 END IF
624 CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
625 $ kz2( kazero( jtype ) ), iasign( jtype ),
626 $ rmagn( kamagn( jtype ) ), ulp,
627 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
628 $ iseed, a, lda )
629 iadd = kadd( kazero( jtype ) )
630 IF( iadd.GT.0 .AND. iadd.LE.n )
631 $ a( iadd, iadd ) = one
632*
633* Generate B (w/o rotation)
634*
635 IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
636 in = 2*( ( n-1 ) / 2 ) + 1
637 IF( in.NE.n )
638 $ CALL dlaset( 'Full', n, n, zero, zero, b, lda )
639 ELSE
640 in = n
641 END IF
642 CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
643 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
644 $ rmagn( kbmagn( jtype ) ), one,
645 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
646 $ iseed, b, lda )
647 iadd = kadd( kbzero( jtype ) )
648 IF( iadd.NE.0 .AND. iadd.LE.n )
649 $ b( iadd, iadd ) = one
650*
651 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
652*
653* Include rotations
654*
655* Generate Q, Z as Householder transformations times
656* a diagonal matrix.
657*
658 DO 40 jc = 1, n - 1
659 DO 30 jr = jc, n
660 q( jr, jc ) = dlarnd( 3, iseed )
661 z( jr, jc ) = dlarnd( 3, iseed )
662 30 CONTINUE
663 CALL dlarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
664 $ work( jc ) )
665 work( 2*n+jc ) = sign( one, q( jc, jc ) )
666 q( jc, jc ) = one
667 CALL dlarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
668 $ work( n+jc ) )
669 work( 3*n+jc ) = sign( one, z( jc, jc ) )
670 z( jc, jc ) = one
671 40 CONTINUE
672 q( n, n ) = one
673 work( n ) = zero
674 work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
675 z( n, n ) = one
676 work( 2*n ) = zero
677 work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
678*
679* Apply the diagonal matrices
680*
681 DO 60 jc = 1, n
682 DO 50 jr = 1, n
683 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
684 $ a( jr, jc )
685 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
686 $ b( jr, jc )
687 50 CONTINUE
688 60 CONTINUE
689 CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
690 $ lda, work( 2*n+1 ), ierr )
691 IF( ierr.NE.0 )
692 $ GO TO 90
693 CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
694 $ a, lda, work( 2*n+1 ), ierr )
695 IF( ierr.NE.0 )
696 $ GO TO 90
697 CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
698 $ lda, work( 2*n+1 ), ierr )
699 IF( ierr.NE.0 )
700 $ GO TO 90
701 CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
702 $ b, lda, work( 2*n+1 ), ierr )
703 IF( ierr.NE.0 )
704 $ GO TO 90
705 END IF
706 ELSE IF (kclass( jtype ).EQ.3) THEN
707*
708* Random matrices
709*
710 DO 80 jc = 1, n
711 DO 70 jr = 1, n
712 a( jr, jc ) = rmagn( kamagn( jtype ) )*
713 $ dlarnd( 2, iseed )
714 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
715 $ dlarnd( 2, iseed )
716 70 CONTINUE
717 80 CONTINUE
718 ELSE
719*
720* Random upper Hessenberg pencil with singular B
721*
722 DO 81 jc = 1, n
723 DO 71 jr = 1, min( jc + 1, n)
724 a( jr, jc ) = rmagn( kamagn( jtype ) )*
725 $ dlarnd( 2, iseed )
726 71 CONTINUE
727 DO 72 jr = jc + 2, n
728 a( jr, jc ) = zero
729 72 CONTINUE
730 81 CONTINUE
731 DO 82 jc = 1, n
732 DO 73 jr = 1, jc
733 b( jr, jc ) = rmagn( kamagn( jtype ) )*
734 $ dlarnd( 2, iseed )
735 73 CONTINUE
736 DO 74 jr = jc + 1, n
737 b( jr, jc ) = zero
738 74 CONTINUE
739 82 CONTINUE
740 DO 83 jc = 1, n, 4
741 b( jc, jc ) = zero
742 83 CONTINUE
743
744 END IF
745*
746 90 CONTINUE
747*
748 IF( ierr.NE.0 ) THEN
749 WRITE( nounit, fmt = 9999 )'Generator', ierr, n, jtype,
750 $ ioldsd
751 info = abs( ierr )
752 RETURN
753 END IF
754*
755 100 CONTINUE
756*
757 DO 110 i = 1, 7
758 result( i ) = -one
759 110 CONTINUE
760*
761* Call XLAENV to set the parameters used in DLAQZ0
762*
763 CALL xlaenv( 12, 10 )
764 CALL xlaenv( 13, 12 )
765 CALL xlaenv( 14, 13 )
766 CALL xlaenv( 15, 2 )
767 CALL xlaenv( 17, 10 )
768*
769* Call DGGEV3 to compute eigenvalues and eigenvectors.
770*
771 CALL dlacpy( ' ', n, n, a, lda, s, lda )
772 CALL dlacpy( ' ', n, n, b, lda, t, lda )
773 CALL dggev3( 'V', 'V', n, s, lda, t, lda, alphar, alphai,
774 $ beta, 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 )'DGGEV31', ierr, n, jtype,
778 $ ioldsd
779 info = abs( ierr )
780 GO TO 190
781 END IF
782*
783* Do the tests (1) and (2)
784*
785 CALL dget52( .true., n, a, lda, b, lda, q, ldq, alphar,
786 $ alphai, beta, work, result( 1 ) )
787 IF( result( 2 ).GT.thresh ) THEN
788 WRITE( nounit, fmt = 9998 )'Left', 'DGGEV31',
789 $ result( 2 ), n, jtype, ioldsd
790 END IF
791*
792* Do the tests (3) and (4)
793*
794 CALL dget52( .false., n, a, lda, b, lda, z, ldq, alphar,
795 $ alphai, beta, work, result( 3 ) )
796 IF( result( 4 ).GT.thresh ) THEN
797 WRITE( nounit, fmt = 9998 )'Right', 'DGGEV31',
798 $ result( 4 ), n, jtype, ioldsd
799 END IF
800*
801* Do the test (5)
802*
803 CALL dlacpy( ' ', n, n, a, lda, s, lda )
804 CALL dlacpy( ' ', n, n, b, lda, t, lda )
805 CALL dggev3( 'N', 'N', n, s, lda, t, lda, alphr1, alphi1,
806 $ beta1, q, ldq, z, ldq, work, lwork, ierr )
807 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
808 result( 1 ) = ulpinv
809 WRITE( nounit, fmt = 9999 )'DGGEV32', ierr, n, jtype,
810 $ ioldsd
811 info = abs( ierr )
812 GO TO 190
813 END IF
814*
815 DO 120 j = 1, n
816 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
817 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 5 )
818 $ = ulpinv
819 120 CONTINUE
820*
821* Do the test (6): Compute eigenvalues and left eigenvectors,
822* and test them
823*
824 CALL dlacpy( ' ', n, n, a, lda, s, lda )
825 CALL dlacpy( ' ', n, n, b, lda, t, lda )
826 CALL dggev3( 'V', 'N', n, s, lda, t, lda, alphr1, alphi1,
827 $ beta1, qe, ldqe, z, ldq, work, lwork, ierr )
828 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
829 result( 1 ) = ulpinv
830 WRITE( nounit, fmt = 9999 )'DGGEV33', ierr, n, jtype,
831 $ ioldsd
832 info = abs( ierr )
833 GO TO 190
834 END IF
835*
836 DO 130 j = 1, n
837 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
838 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 6 )
839 $ = ulpinv
840 130 CONTINUE
841*
842 DO 150 j = 1, n
843 DO 140 jc = 1, n
844 IF( q( j, jc ).NE.qe( j, jc ) )
845 $ result( 6 ) = ulpinv
846 140 CONTINUE
847 150 CONTINUE
848*
849* DO the test (7): Compute eigenvalues and right eigenvectors,
850* and test them
851*
852 CALL dlacpy( ' ', n, n, a, lda, s, lda )
853 CALL dlacpy( ' ', n, n, b, lda, t, lda )
854 CALL dggev3( 'N', 'V', n, s, lda, t, lda, alphr1, alphi1,
855 $ beta1, q, ldq, qe, ldqe, work, lwork, ierr )
856 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
857 result( 1 ) = ulpinv
858 WRITE( nounit, fmt = 9999 )'DGGEV34', ierr, n, jtype,
859 $ ioldsd
860 info = abs( ierr )
861 GO TO 190
862 END IF
863*
864 DO 160 j = 1, n
865 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
866 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 7 )
867 $ = ulpinv
868 160 CONTINUE
869*
870 DO 180 j = 1, n
871 DO 170 jc = 1, n
872 IF( z( j, jc ).NE.qe( j, jc ) )
873 $ result( 7 ) = ulpinv
874 170 CONTINUE
875 180 CONTINUE
876*
877* End of Loop -- Check for RESULT(j) > THRESH
878*
879 190 CONTINUE
880*
881 ntestt = ntestt + 7
882*
883* Print out tests which fail.
884*
885 DO 200 jr = 1, 7
886 IF( result( jr ).GE.thresh ) THEN
887*
888* If this is the first test to fail,
889* print a header to the data file.
890*
891 IF( nerrs.EQ.0 ) THEN
892 WRITE( nounit, fmt = 9997 )'DGV'
893*
894* Matrix types
895*
896 WRITE( nounit, fmt = 9996 )
897 WRITE( nounit, fmt = 9995 )
898 WRITE( nounit, fmt = 9994 )'Orthogonal'
899*
900* Tests performed
901*
902 WRITE( nounit, fmt = 9993 )
903*
904 END IF
905 nerrs = nerrs + 1
906 IF( result( jr ).LT.10000.0d0 ) THEN
907 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
908 $ result( jr )
909 ELSE
910 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
911 $ result( jr )
912 END IF
913 END IF
914 200 CONTINUE
915*
916 210 CONTINUE
917 220 CONTINUE
918*
919* Summary
920*
921 CALL alasvm( 'DGV', nounit, nerrs, ntestt, 0 )
922*
923 work( 1 ) = maxwrk
924*
925 RETURN
926*
927 9999 FORMAT( ' DDRGEV3: ', a, ' returned INFO=', i6, '.', / 3x, 'N=',
928 $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
929*
930 9998 FORMAT( ' DDRGEV3: ', a, ' Eigenvectors from ', a,
931 $ ' incorrectly normalized.', / ' Bits of error=', 0p, g10.3,
932 $ ',', 3x, 'N=', i4, ', JTYPE=', i3, ', ISEED=(',
933 $ 4( i4, ',' ), i5, ')' )
934*
935 9997 FORMAT( / 1x, a3, ' -- Real Generalized eigenvalue problem driver'
936 $ )
937*
938 9996 FORMAT( ' Matrix types (see DDRGEV3 for details): ' )
939*
940 9995 FORMAT( ' Special Matrices:', 23x,
941 $ '(J''=transposed Jordan block)',
942 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
943 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
944 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
945 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
946 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
947 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
948 9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
949 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
950 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
951 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
952 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
953 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
954 $ '23=(small,large) 24=(small,small) 25=(large,large)',
955 $ / ' 26=random O(1) matrices.' )
956*
957 9993 FORMAT( / ' Tests performed: ',
958 $ / ' 1 = max | ( b A - a B )''*l | / const.,',
959 $ / ' 2 = | |VR(i)| - 1 | / ulp,',
960 $ / ' 3 = max | ( b A - a B )*r | / const.',
961 $ / ' 4 = | |VL(i)| - 1 | / ulp,',
962 $ / ' 5 = 0 if W same no matter if r or l computed,',
963 $ / ' 6 = 0 if l same no matter if l computed,',
964 $ / ' 7 = 0 if r same no matter if r computed,', / 1x )
965 9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
966 $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
967 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
968 $ 4( i4, ',' ), ' result ', i2, ' is', 1p, d10.3 )
969*
970* End of DDRGEV3
971*
972 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 xlaenv(ISPEC, NVALUE)
XLAENV
Definition: xlaenv.f:81
subroutine dlatm4(ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
DLATM4
Definition: dlatm4.f:175
subroutine dget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
DGET52
Definition: dget52.f:199
subroutine ddrgev3(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)
DDRGEV3
Definition: ddrgev3.f:408
subroutine dggev3(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
DGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (...
Definition: dggev3.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