LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ddrges.f
Go to the documentation of this file.
1*> \brief \b DDRGES
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 DDRGES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12* NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHAR,
13* ALPHAI, BETA, WORK, LWORK, RESULT, BWORK,
14* INFO )
15*
16* .. Scalar Arguments ..
17* INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
18* DOUBLE PRECISION THRESH
19* ..
20* .. Array Arguments ..
21* LOGICAL BWORK( * ), DOTYPE( * )
22* INTEGER ISEED( 4 ), NN( * )
23* DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
24* $ B( LDA, * ), BETA( * ), Q( LDQ, * ),
25* $ RESULT( 13 ), S( LDA, * ), T( LDA, * ),
26* $ WORK( * ), Z( LDQ, * )
27* ..
28*
29*
30*> \par Purpose:
31* =============
32*>
33*> \verbatim
34*>
35*> DDRGES checks the nonsymmetric generalized eigenvalue (Schur form)
36*> problem driver DGGES.
37*>
38*> DGGES factors A and B as Q S Z' and Q T Z' , where ' means
39*> transpose, T is upper triangular, S is in generalized Schur form
40*> (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
41*> the 2x2 blocks corresponding to complex conjugate pairs of
42*> generalized eigenvalues), and Q and Z are orthogonal. It also
43*> computes the generalized eigenvalues (alpha(j),beta(j)), j=1,...,n,
44*> Thus, w(j) = alpha(j)/beta(j) is a root of the characteristic
45*> equation
46*> det( A - w(j) B ) = 0
47*> Optionally it also reorder the eigenvalues so that a selected
48*> cluster of eigenvalues appears in the leading diagonal block of the
49*> Schur forms.
50*>
51*> When DDRGES is called, a number of matrix "sizes" ("N's") and a
52*> number of matrix "TYPES" are specified. For each size ("N")
53*> and each TYPE of matrix, a pair of matrices (A, B) will be generated
54*> and used for testing. For each matrix pair, the following 13 tests
55*> will be performed and compared with the threshold THRESH except
56*> the tests (5), (11) and (13).
57*>
58*>
59*> (1) | A - Q S Z' | / ( |A| n ulp ) (no sorting of eigenvalues)
60*>
61*>
62*> (2) | B - Q T Z' | / ( |B| n ulp ) (no sorting of eigenvalues)
63*>
64*>
65*> (3) | I - QQ' | / ( n ulp ) (no sorting of eigenvalues)
66*>
67*>
68*> (4) | I - ZZ' | / ( n ulp ) (no sorting of eigenvalues)
69*>
70*> (5) if A is in Schur form (i.e. quasi-triangular form)
71*> (no sorting of eigenvalues)
72*>
73*> (6) if eigenvalues = diagonal blocks of the Schur form (S, T),
74*> i.e., test the maximum over j of D(j) where:
75*>
76*> if alpha(j) is real:
77*> |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
78*> D(j) = ------------------------ + -----------------------
79*> max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
80*>
81*> if alpha(j) is complex:
82*> | det( s S - w T ) |
83*> D(j) = ---------------------------------------------------
84*> ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
85*>
86*> and S and T are here the 2 x 2 diagonal blocks of S and T
87*> corresponding to the j-th and j+1-th eigenvalues.
88*> (no sorting of eigenvalues)
89*>
90*> (7) | (A,B) - Q (S,T) Z' | / ( | (A,B) | n ulp )
91*> (with sorting of eigenvalues).
92*>
93*> (8) | I - QQ' | / ( n ulp ) (with sorting of eigenvalues).
94*>
95*> (9) | I - ZZ' | / ( n ulp ) (with sorting of eigenvalues).
96*>
97*> (10) if A is in Schur form (i.e. quasi-triangular form)
98*> (with sorting of eigenvalues).
99*>
100*> (11) if eigenvalues = diagonal blocks of the Schur form (S, T),
101*> i.e. test the maximum over j of D(j) where:
102*>
103*> if alpha(j) is real:
104*> |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
105*> D(j) = ------------------------ + -----------------------
106*> max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
107*>
108*> if alpha(j) is complex:
109*> | det( s S - w T ) |
110*> D(j) = ---------------------------------------------------
111*> ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
112*>
113*> and S and T are here the 2 x 2 diagonal blocks of S and T
114*> corresponding to the j-th and j+1-th eigenvalues.
115*> (with sorting of eigenvalues).
116*>
117*> (12) if sorting worked and SDIM is the number of eigenvalues
118*> which were SELECTed.
119*>
120*> Test Matrices
121*> =============
122*>
123*> The sizes of the test matrices are specified by an array
124*> NN(1:NSIZES); the value of each element NN(j) specifies one size.
125*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
126*> DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
127*> Currently, the list of possible types is:
128*>
129*> (1) ( 0, 0 ) (a pair of zero matrices)
130*>
131*> (2) ( I, 0 ) (an identity and a zero matrix)
132*>
133*> (3) ( 0, I ) (an identity and a zero matrix)
134*>
135*> (4) ( I, I ) (a pair of identity matrices)
136*>
137*> t t
138*> (5) ( J , J ) (a pair of transposed Jordan blocks)
139*>
140*> t ( I 0 )
141*> (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
142*> ( 0 I ) ( 0 J )
143*> and I is a k x k identity and J a (k+1)x(k+1)
144*> Jordan block; k=(N-1)/2
145*>
146*> (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal
147*> matrix with those diagonal entries.)
148*> (8) ( I, D )
149*>
150*> (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
151*>
152*> (10) ( small*D, big*I )
153*>
154*> (11) ( big*I, small*D )
155*>
156*> (12) ( small*I, big*D )
157*>
158*> (13) ( big*D, big*I )
159*>
160*> (14) ( small*D, small*I )
161*>
162*> (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
163*> D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
164*> t t
165*> (16) Q ( J , J ) Z where Q and Z are random orthogonal matrices.
166*>
167*> (17) Q ( T1, T2 ) Z where T1 and T2 are upper triangular matrices
168*> with random O(1) entries above the diagonal
169*> and diagonal entries diag(T1) =
170*> ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
171*> ( 0, N-3, N-4,..., 1, 0, 0 )
172*>
173*> (18) Q ( T1, T2 ) Z diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
174*> diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
175*> s = machine precision.
176*>
177*> (19) Q ( T1, T2 ) Z diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
178*> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
179*>
180*> N-5
181*> (20) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
182*> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
183*>
184*> (21) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
185*> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
186*> where r1,..., r(N-4) are random.
187*>
188*> (22) Q ( big*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
189*> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
190*>
191*> (23) Q ( small*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
192*> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
193*>
194*> (24) Q ( small*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
195*> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
196*>
197*> (25) Q ( big*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
198*> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
199*>
200*> (26) Q ( T1, T2 ) Z where T1 and T2 are random upper-triangular
201*> matrices.
202*>
203*> \endverbatim
204*
205* Arguments:
206* ==========
207*
208*> \param[in] NSIZES
209*> \verbatim
210*> NSIZES is INTEGER
211*> The number of sizes of matrices to use. If it is zero,
212*> DDRGES does nothing. NSIZES >= 0.
213*> \endverbatim
214*>
215*> \param[in] NN
216*> \verbatim
217*> NN is INTEGER array, dimension (NSIZES)
218*> An array containing the sizes to be used for the matrices.
219*> Zero values will be skipped. NN >= 0.
220*> \endverbatim
221*>
222*> \param[in] NTYPES
223*> \verbatim
224*> NTYPES is INTEGER
225*> The number of elements in DOTYPE. If it is zero, DDRGES
226*> does nothing. It must be at least zero. If it is MAXTYP+1
227*> and NSIZES is 1, then an additional type, MAXTYP+1 is
228*> defined, which is to use whatever matrix is in A on input.
229*> This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
230*> DOTYPE(MAXTYP+1) is .TRUE. .
231*> \endverbatim
232*>
233*> \param[in] DOTYPE
234*> \verbatim
235*> DOTYPE is LOGICAL array, dimension (NTYPES)
236*> If DOTYPE(j) is .TRUE., then for each size in NN a
237*> matrix of that size and of type j will be generated.
238*> If NTYPES is smaller than the maximum number of types
239*> defined (PARAMETER MAXTYP), then types NTYPES+1 through
240*> MAXTYP will not be generated. If NTYPES is larger
241*> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
242*> will be ignored.
243*> \endverbatim
244*>
245*> \param[in,out] ISEED
246*> \verbatim
247*> ISEED is INTEGER array, dimension (4)
248*> On entry ISEED specifies the seed of the random number
249*> generator. The array elements should be between 0 and 4095;
250*> if not they will be reduced mod 4096. Also, ISEED(4) must
251*> be odd. The random number generator uses a linear
252*> congruential sequence limited to small integers, and so
253*> should produce machine independent random numbers. The
254*> values of ISEED are changed on exit, and can be used in the
255*> next call to DDRGES to continue the same random number
256*> sequence.
257*> \endverbatim
258*>
259*> \param[in] THRESH
260*> \verbatim
261*> THRESH is DOUBLE PRECISION
262*> A test will count as "failed" if the "error", computed as
263*> described above, exceeds THRESH. Note that the error is
264*> scaled to be O(1), so THRESH should be a reasonably small
265*> multiple of 1, e.g., 10 or 100. In particular, it should
266*> not depend on the precision (single vs. double) or the size
267*> of the matrix. THRESH >= 0.
268*> \endverbatim
269*>
270*> \param[in] NOUNIT
271*> \verbatim
272*> NOUNIT is INTEGER
273*> The FORTRAN unit number for printing out error messages
274*> (e.g., if a routine returns IINFO not equal to 0.)
275*> \endverbatim
276*>
277*> \param[in,out] A
278*> \verbatim
279*> A is DOUBLE PRECISION array,
280*> dimension(LDA, max(NN))
281*> Used to hold the original A matrix. Used as input only
282*> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
283*> DOTYPE(MAXTYP+1)=.TRUE.
284*> \endverbatim
285*>
286*> \param[in] LDA
287*> \verbatim
288*> LDA is INTEGER
289*> The leading dimension of A, B, S, and T.
290*> It must be at least 1 and at least max( NN ).
291*> \endverbatim
292*>
293*> \param[in,out] B
294*> \verbatim
295*> B is DOUBLE PRECISION array,
296*> dimension(LDA, max(NN))
297*> Used to hold the original B matrix. Used as input only
298*> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
299*> DOTYPE(MAXTYP+1)=.TRUE.
300*> \endverbatim
301*>
302*> \param[out] S
303*> \verbatim
304*> S is DOUBLE PRECISION array, dimension (LDA, max(NN))
305*> The Schur form matrix computed from A by DGGES. On exit, S
306*> contains the Schur form matrix corresponding to the matrix
307*> in A.
308*> \endverbatim
309*>
310*> \param[out] T
311*> \verbatim
312*> T is DOUBLE PRECISION array, dimension (LDA, max(NN))
313*> The upper triangular matrix computed from B by DGGES.
314*> \endverbatim
315*>
316*> \param[out] Q
317*> \verbatim
318*> Q is DOUBLE PRECISION array, dimension (LDQ, max(NN))
319*> The (left) orthogonal matrix computed by DGGES.
320*> \endverbatim
321*>
322*> \param[in] LDQ
323*> \verbatim
324*> LDQ is INTEGER
325*> The leading dimension of Q and Z. It must
326*> be at least 1 and at least max( NN ).
327*> \endverbatim
328*>
329*> \param[out] Z
330*> \verbatim
331*> Z is DOUBLE PRECISION array, dimension( LDQ, max(NN) )
332*> The (right) orthogonal matrix computed by DGGES.
333*> \endverbatim
334*>
335*> \param[out] ALPHAR
336*> \verbatim
337*> ALPHAR is DOUBLE PRECISION array, dimension (max(NN))
338*> \endverbatim
339*>
340*> \param[out] ALPHAI
341*> \verbatim
342*> ALPHAI is DOUBLE PRECISION array, dimension (max(NN))
343*> \endverbatim
344*>
345*> \param[out] BETA
346*> \verbatim
347*> BETA is DOUBLE PRECISION array, dimension (max(NN))
348*>
349*> The generalized eigenvalues of (A,B) computed by DGGES.
350*> ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
351*> generalized eigenvalue of A and B.
352*> \endverbatim
353*>
354*> \param[out] WORK
355*> \verbatim
356*> WORK is DOUBLE PRECISION array, dimension (LWORK)
357*> \endverbatim
358*>
359*> \param[in] LWORK
360*> \verbatim
361*> LWORK is INTEGER
362*> The dimension of the array WORK.
363*> LWORK >= MAX( 10*(N+1), 3*N*N ), where N is the largest
364*> matrix dimension.
365*> \endverbatim
366*>
367*> \param[out] RESULT
368*> \verbatim
369*> RESULT is DOUBLE PRECISION array, dimension (15)
370*> The values computed by the tests described above.
371*> The values are currently limited to 1/ulp, to avoid overflow.
372*> \endverbatim
373*>
374*> \param[out] BWORK
375*> \verbatim
376*> BWORK is LOGICAL array, dimension (N)
377*> \endverbatim
378*>
379*> \param[out] INFO
380*> \verbatim
381*> INFO is INTEGER
382*> = 0: successful exit
383*> < 0: if INFO = -i, the i-th argument had an illegal value.
384*> > 0: A routine returned an error code. INFO is the
385*> absolute value of the INFO value returned.
386*> \endverbatim
387*
388* Authors:
389* ========
390*
391*> \author Univ. of Tennessee
392*> \author Univ. of California Berkeley
393*> \author Univ. of Colorado Denver
394*> \author NAG Ltd.
395*
396*> \ingroup double_eig
397*
398* =====================================================================
399 SUBROUTINE ddrges( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
400 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHAR,
401 $ ALPHAI, BETA, WORK, LWORK, RESULT, BWORK,
402 $ INFO )
403*
404* -- LAPACK test routine --
405* -- LAPACK is a software package provided by Univ. of Tennessee, --
406* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
407*
408* .. Scalar Arguments ..
409 INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
410 DOUBLE PRECISION THRESH
411* ..
412* .. Array Arguments ..
413 LOGICAL BWORK( * ), DOTYPE( * )
414 INTEGER ISEED( 4 ), NN( * )
415 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
416 $ b( lda, * ), beta( * ), q( ldq, * ),
417 $ result( 13 ), s( lda, * ), t( lda, * ),
418 $ work( * ), z( ldq, * )
419* ..
420*
421* =====================================================================
422*
423* .. Parameters ..
424 DOUBLE PRECISION ZERO, ONE
425 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
426 INTEGER MAXTYP
427 parameter( maxtyp = 26 )
428* ..
429* .. Local Scalars ..
430 LOGICAL BADNN, ILABAD
431 CHARACTER SORT
432 INTEGER I, I1, IADD, IERR, IINFO, IN, ISORT, J, JC, JR,
433 $ jsize, jtype, knteig, maxwrk, minwrk, mtypes,
434 $ n, n1, nb, nerrs, nmats, nmax, ntest, ntestt,
435 $ rsub, sdim
436 DOUBLE PRECISION SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV
437* ..
438* .. Local Arrays ..
439 INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
440 $ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
441 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
442 $ kbmagn( maxtyp ), kbtype( maxtyp ),
443 $ kbzero( maxtyp ), kclass( maxtyp ),
444 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
445 DOUBLE PRECISION RMAGN( 0: 3 )
446* ..
447* .. External Functions ..
448 LOGICAL DLCTES
449 INTEGER ILAENV
450 DOUBLE PRECISION DLAMCH, DLARND
451 EXTERNAL dlctes, ilaenv, dlamch, dlarnd
452* ..
453* .. External Subroutines ..
454 EXTERNAL alasvm, dget51, dget53, dget54, dgges, dlacpy,
456* ..
457* .. Intrinsic Functions ..
458 INTRINSIC abs, dble, max, min, sign
459* ..
460* .. Data statements ..
461 DATA kclass / 15*1, 10*2, 1*3 /
462 DATA kz1 / 0, 1, 2, 1, 3, 3 /
463 DATA kz2 / 0, 0, 1, 2, 1, 1 /
464 DATA kadd / 0, 0, 0, 0, 3, 2 /
465 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
466 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
467 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
468 $ 1, 1, -4, 2, -4, 8*8, 0 /
469 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
470 $ 4*5, 4*3, 1 /
471 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
472 $ 4*6, 4*4, 1 /
473 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
474 $ 2, 1 /
475 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
476 $ 2, 1 /
477 DATA ktrian / 16*0, 10*1 /
478 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
479 $ 5*2, 0 /
480 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
481* ..
482* .. Executable Statements ..
483*
484* Check for errors
485*
486 info = 0
487*
488 badnn = .false.
489 nmax = 1
490 DO 10 j = 1, nsizes
491 nmax = max( nmax, nn( j ) )
492 IF( nn( j ).LT.0 )
493 $ badnn = .true.
494 10 CONTINUE
495*
496 IF( nsizes.LT.0 ) THEN
497 info = -1
498 ELSE IF( badnn ) THEN
499 info = -2
500 ELSE IF( ntypes.LT.0 ) THEN
501 info = -3
502 ELSE IF( thresh.LT.zero ) THEN
503 info = -6
504 ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
505 info = -9
506 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
507 info = -14
508 END IF
509*
510* Compute workspace
511* (Note: Comments in the code beginning "Workspace:" describe the
512* minimal amount of workspace needed at that point in the code,
513* as well as the preferred amount for good performance.
514* NB refers to the optimal block size for the immediately
515* following subroutine, as returned by ILAENV.
516*
517 minwrk = 1
518 IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
519 minwrk = max( 10*( nmax+1 ), 3*nmax*nmax )
520 nb = max( 1, ilaenv( 1, 'DGEQRF', ' ', nmax, nmax, -1, -1 ),
521 $ ilaenv( 1, 'DORMQR', 'LT', nmax, nmax, nmax, -1 ),
522 $ ilaenv( 1, 'DORGQR', ' ', nmax, nmax, nmax, -1 ) )
523 maxwrk = max( 10*( nmax+1 ), 2*nmax+nmax*nb, 3*nmax*nmax )
524 work( 1 ) = maxwrk
525 END IF
526*
527 IF( lwork.LT.minwrk )
528 $ info = -20
529*
530 IF( info.NE.0 ) THEN
531 CALL xerbla( 'DDRGES', -info )
532 RETURN
533 END IF
534*
535* Quick return if possible
536*
537 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
538 $ RETURN
539*
540 safmin = dlamch( 'Safe minimum' )
541 ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
542 safmin = safmin / ulp
543 safmax = one / safmin
544 ulpinv = one / ulp
545*
546* The values RMAGN(2:3) depend on N, see below.
547*
548 rmagn( 0 ) = zero
549 rmagn( 1 ) = one
550*
551* Loop over matrix sizes
552*
553 ntestt = 0
554 nerrs = 0
555 nmats = 0
556*
557 DO 190 jsize = 1, nsizes
558 n = nn( jsize )
559 n1 = max( 1, n )
560 rmagn( 2 ) = safmax*ulp / dble( n1 )
561 rmagn( 3 ) = safmin*ulpinv*dble( n1 )
562*
563 IF( nsizes.NE.1 ) THEN
564 mtypes = min( maxtyp, ntypes )
565 ELSE
566 mtypes = min( maxtyp+1, ntypes )
567 END IF
568*
569* Loop over matrix types
570*
571 DO 180 jtype = 1, mtypes
572 IF( .NOT.dotype( jtype ) )
573 $ GO TO 180
574 nmats = nmats + 1
575 ntest = 0
576*
577* Save ISEED in case of an error.
578*
579 DO 20 j = 1, 4
580 ioldsd( j ) = iseed( j )
581 20 CONTINUE
582*
583* Initialize RESULT
584*
585 DO 30 j = 1, 13
586 result( j ) = zero
587 30 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 110
614 iinfo = 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 50 jc = 1, n - 1
661 DO 40 jr = jc, n
662 q( jr, jc ) = dlarnd( 3, iseed )
663 z( jr, jc ) = dlarnd( 3, iseed )
664 40 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 50 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 70 jc = 1, n
684 DO 60 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 60 CONTINUE
690 70 CONTINUE
691 CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
692 $ lda, work( 2*n+1 ), iinfo )
693 IF( iinfo.NE.0 )
694 $ GO TO 100
695 CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
696 $ a, lda, work( 2*n+1 ), iinfo )
697 IF( iinfo.NE.0 )
698 $ GO TO 100
699 CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
700 $ lda, work( 2*n+1 ), iinfo )
701 IF( iinfo.NE.0 )
702 $ GO TO 100
703 CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
704 $ b, lda, work( 2*n+1 ), iinfo )
705 IF( iinfo.NE.0 )
706 $ GO TO 100
707 END IF
708 ELSE
709*
710* Random matrices
711*
712 DO 90 jc = 1, n
713 DO 80 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 80 CONTINUE
719 90 CONTINUE
720 END IF
721*
722 100 CONTINUE
723*
724 IF( iinfo.NE.0 ) THEN
725 WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
726 $ ioldsd
727 info = abs( iinfo )
728 RETURN
729 END IF
730*
731 110 CONTINUE
732*
733 DO 120 i = 1, 13
734 result( i ) = -one
735 120 CONTINUE
736*
737* Test with and without sorting of eigenvalues
738*
739 DO 150 isort = 0, 1
740 IF( isort.EQ.0 ) THEN
741 sort = 'N'
742 rsub = 0
743 ELSE
744 sort = 'S'
745 rsub = 5
746 END IF
747*
748* Call DGGES to compute H, T, Q, Z, alpha, and beta.
749*
750 CALL dlacpy( 'Full', n, n, a, lda, s, lda )
751 CALL dlacpy( 'Full', n, n, b, lda, t, lda )
752 ntest = 1 + rsub + isort
753 result( 1+rsub+isort ) = ulpinv
754 CALL dgges( 'V', 'V', sort, dlctes, n, s, lda, t, lda,
755 $ sdim, alphar, alphai, beta, q, ldq, z, ldq,
756 $ work, lwork, bwork, iinfo )
757 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
758 result( 1+rsub+isort ) = ulpinv
759 WRITE( nounit, fmt = 9999 )'DGGES', iinfo, n, jtype,
760 $ ioldsd
761 info = abs( iinfo )
762 GO TO 160
763 END IF
764*
765 ntest = 4 + rsub
766*
767* Do tests 1--4 (or tests 7--9 when reordering )
768*
769 IF( isort.EQ.0 ) THEN
770 CALL dget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
771 $ work, result( 1 ) )
772 CALL dget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
773 $ work, result( 2 ) )
774 ELSE
775 CALL dget54( n, a, lda, b, lda, s, lda, t, lda, q,
776 $ ldq, z, ldq, work, result( 7 ) )
777 END IF
778 CALL dget51( 3, n, a, lda, t, lda, q, ldq, q, ldq, work,
779 $ result( 3+rsub ) )
780 CALL dget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
781 $ result( 4+rsub ) )
782*
783* Do test 5 and 6 (or Tests 10 and 11 when reordering):
784* check Schur form of A and compare eigenvalues with
785* diagonals.
786*
787 ntest = 6 + rsub
788 temp1 = zero
789*
790 DO 130 j = 1, n
791 ilabad = .false.
792 IF( alphai( j ).EQ.zero ) THEN
793 temp2 = ( abs( alphar( j )-s( j, j ) ) /
794 $ max( safmin, abs( alphar( j ) ), abs( s( j,
795 $ j ) ) )+abs( beta( j )-t( j, j ) ) /
796 $ max( safmin, abs( beta( j ) ), abs( t( j,
797 $ j ) ) ) ) / ulp
798*
799 IF( j.LT.n ) THEN
800 IF( s( j+1, j ).NE.zero ) THEN
801 ilabad = .true.
802 result( 5+rsub ) = ulpinv
803 END IF
804 END IF
805 IF( j.GT.1 ) THEN
806 IF( s( j, j-1 ).NE.zero ) THEN
807 ilabad = .true.
808 result( 5+rsub ) = ulpinv
809 END IF
810 END IF
811*
812 ELSE
813 IF( alphai( j ).GT.zero ) THEN
814 i1 = j
815 ELSE
816 i1 = j - 1
817 END IF
818 IF( i1.LE.0 .OR. i1.GE.n ) THEN
819 ilabad = .true.
820 ELSE IF( i1.LT.n-1 ) THEN
821 IF( s( i1+2, i1+1 ).NE.zero ) THEN
822 ilabad = .true.
823 result( 5+rsub ) = ulpinv
824 END IF
825 ELSE IF( i1.GT.1 ) THEN
826 IF( s( i1, i1-1 ).NE.zero ) THEN
827 ilabad = .true.
828 result( 5+rsub ) = ulpinv
829 END IF
830 END IF
831 IF( .NOT.ilabad ) THEN
832 CALL dget53( s( i1, i1 ), lda, t( i1, i1 ), lda,
833 $ beta( j ), alphar( j ),
834 $ alphai( j ), temp2, ierr )
835 IF( ierr.GE.3 ) THEN
836 WRITE( nounit, fmt = 9998 )ierr, j, n,
837 $ jtype, ioldsd
838 info = abs( ierr )
839 END IF
840 ELSE
841 temp2 = ulpinv
842 END IF
843*
844 END IF
845 temp1 = max( temp1, temp2 )
846 IF( ilabad ) THEN
847 WRITE( nounit, fmt = 9997 )j, n, jtype, ioldsd
848 END IF
849 130 CONTINUE
850 result( 6+rsub ) = temp1
851*
852 IF( isort.GE.1 ) THEN
853*
854* Do test 12
855*
856 ntest = 12
857 result( 12 ) = zero
858 knteig = 0
859 DO 140 i = 1, n
860 IF( dlctes( alphar( i ), alphai( i ),
861 $ beta( i ) ) .OR. dlctes( alphar( i ),
862 $ -alphai( i ), beta( i ) ) ) THEN
863 knteig = knteig + 1
864 END IF
865 IF( i.LT.n ) THEN
866 IF( ( dlctes( alphar( i+1 ), alphai( i+1 ),
867 $ beta( i+1 ) ) .OR. dlctes( alphar( i+1 ),
868 $ -alphai( i+1 ), beta( i+1 ) ) ) .AND.
869 $ ( .NOT.( dlctes( alphar( i ), alphai( i ),
870 $ beta( i ) ) .OR. dlctes( alphar( i ),
871 $ -alphai( i ), beta( i ) ) ) ) .AND.
872 $ iinfo.NE.n+2 ) THEN
873 result( 12 ) = ulpinv
874 END IF
875 END IF
876 140 CONTINUE
877 IF( sdim.NE.knteig ) THEN
878 result( 12 ) = ulpinv
879 END IF
880 END IF
881*
882 150 CONTINUE
883*
884* End of Loop -- Check for RESULT(j) > THRESH
885*
886 160 CONTINUE
887*
888 ntestt = ntestt + ntest
889*
890* Print out tests which fail.
891*
892 DO 170 jr = 1, ntest
893 IF( result( jr ).GE.thresh ) THEN
894*
895* If this is the first test to fail,
896* print a header to the data file.
897*
898 IF( nerrs.EQ.0 ) THEN
899 WRITE( nounit, fmt = 9996 )'DGS'
900*
901* Matrix types
902*
903 WRITE( nounit, fmt = 9995 )
904 WRITE( nounit, fmt = 9994 )
905 WRITE( nounit, fmt = 9993 )'Orthogonal'
906*
907* Tests performed
908*
909 WRITE( nounit, fmt = 9992 )'orthogonal', '''',
910 $ 'transpose', ( '''', j = 1, 8 )
911*
912 END IF
913 nerrs = nerrs + 1
914 IF( result( jr ).LT.10000.0d0 ) THEN
915 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
916 $ result( jr )
917 ELSE
918 WRITE( nounit, fmt = 9990 )n, jtype, ioldsd, jr,
919 $ result( jr )
920 END IF
921 END IF
922 170 CONTINUE
923*
924 180 CONTINUE
925 190 CONTINUE
926*
927* Summary
928*
929 CALL alasvm( 'DGS', nounit, nerrs, ntestt, 0 )
930*
931 work( 1 ) = maxwrk
932*
933 RETURN
934*
935 9999 FORMAT( ' DDRGES: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
936 $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
937*
938 9998 FORMAT( ' DDRGES: DGET53 returned INFO=', i1, ' for eigenvalue ',
939 $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(',
940 $ 4( i4, ',' ), i5, ')' )
941*
942 9997 FORMAT( ' DDRGES: S not in Schur form at eigenvalue ', i6, '.',
943 $ / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
944 $ i5, ')' )
945*
946 9996 FORMAT( / 1x, a3, ' -- Real Generalized Schur form driver' )
947*
948 9995 FORMAT( ' Matrix types (see DDRGES for details): ' )
949*
950 9994 FORMAT( ' Special Matrices:', 23x,
951 $ '(J''=transposed Jordan block)',
952 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
953 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
954 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
955 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
956 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
957 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
958 9993 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
959 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
960 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
961 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
962 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
963 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
964 $ '23=(small,large) 24=(small,small) 25=(large,large)',
965 $ / ' 26=random O(1) matrices.' )
966*
967 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
968 $ 'Q and Z are ', a, ',', / 19x,
969 $ 'l and r are the appropriate left and right', / 19x,
970 $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
971 $ ' means ', a, '.)', / ' Without ordering: ',
972 $ / ' 1 = | A - Q S Z', a,
973 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
974 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
975 $ ' | / ( n ulp ) 4 = | I - ZZ', a,
976 $ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
977 $ / ' 6 = difference between (alpha,beta)',
978 $ ' and diagonals of (S,T)', / ' With ordering: ',
979 $ / ' 7 = | (A,B) - Q (S,T) Z', a,
980 $ ' | / ( |(A,B)| n ulp ) ', / ' 8 = | I - QQ', a,
981 $ ' | / ( n ulp ) 9 = | I - ZZ', a,
982 $ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
983 $ / ' 11 = difference between (alpha,beta) and diagonals',
984 $ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
985 $ 'selected eigenvalues', / )
986 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
987 $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
988 9990 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
989 $ 4( i4, ',' ), ' result ', i2, ' is', 1p, d10.3 )
990*
991* End of DDRGES
992*
993 END
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ddrges(nsizes, nn, ntypes, dotype, iseed, thresh, nounit, a, lda, b, s, t, q, ldq, z, alphar, alphai, beta, work, lwork, result, bwork, info)
DDRGES
Definition ddrges.f:403
subroutine dget51(itype, n, a, lda, b, ldb, u, ldu, v, ldv, work, result)
DGET51
Definition dget51.f:149
subroutine dget53(a, lda, b, ldb, scale, wr, wi, result, info)
DGET53
Definition dget53.f:126
subroutine dget54(n, a, lda, b, ldb, s, lds, t, ldt, u, ldu, v, ldv, work, result)
DGET54
Definition dget54.f:156
subroutine dlatm4(itype, n, nz1, nz2, isign, amagn, rcond, triang, idist, iseed, a, lda)
DLATM4
Definition dlatm4.f:175
subroutine dgges(jobvsl, jobvsr, sort, selctg, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, bwork, info)
DGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE m...
Definition dgges.f:284
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 dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106
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 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