LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sdrges3.f
Go to the documentation of this file.
1*> \brief \b SDRGES3
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 SDRGES3( 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* REAL THRESH
19* ..
20* .. Array Arguments ..
21* LOGICAL BWORK( * ), DOTYPE( * )
22* INTEGER ISEED( 4 ), NN( * )
23* REAL 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*> SDRGES3 checks the nonsymmetric generalized eigenvalue (Schur form)
36*> problem driver SGGES3.
37*>
38*> SGGES3 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 SDRGES3 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*> SDRGES3 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, SDRGES3
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 SDRGES3 to continue the same random number
256*> sequence.
257*> \endverbatim
258*>
259*> \param[in] THRESH
260*> \verbatim
261*> THRESH is REAL
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 REAL 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 REAL 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 REAL array, dimension (LDA, max(NN))
305*> The Schur form matrix computed from A by SGGES3. 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 REAL array, dimension (LDA, max(NN))
313*> The upper triangular matrix computed from B by SGGES3.
314*> \endverbatim
315*>
316*> \param[out] Q
317*> \verbatim
318*> Q is REAL array, dimension (LDQ, max(NN))
319*> The (left) orthogonal matrix computed by SGGES3.
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 REAL array, dimension( LDQ, max(NN) )
332*> The (right) orthogonal matrix computed by SGGES3.
333*> \endverbatim
334*>
335*> \param[out] ALPHAR
336*> \verbatim
337*> ALPHAR is REAL array, dimension (max(NN))
338*> \endverbatim
339*>
340*> \param[out] ALPHAI
341*> \verbatim
342*> ALPHAI is REAL array, dimension (max(NN))
343*> \endverbatim
344*>
345*> \param[out] BETA
346*> \verbatim
347*> BETA is REAL array, dimension (max(NN))
348*>
349*> The generalized eigenvalues of (A,B) computed by SGGES3.
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 REAL 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 REAL 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 single_eig
397*
398* =====================================================================
399 SUBROUTINE sdrges3( 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 REAL THRESH
411* ..
412* .. Array Arguments ..
413 LOGICAL BWORK( * ), DOTYPE( * )
414 INTEGER ISEED( 4 ), NN( * )
415 REAL 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 REAL ZERO, ONE
425 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+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 REAL 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 REAL RMAGN( 0: 3 )
446* ..
447* .. External Functions ..
448 LOGICAL SLCTES
449 INTEGER ILAENV
450 REAL SLAMCH, SLARND
451 EXTERNAL slctes, ilaenv, slamch, slarnd
452* ..
453* .. External Subroutines ..
454 EXTERNAL alasvm, sget51, sget53, sget54, sgges3, slabad,
456* ..
457* .. Intrinsic Functions ..
458 INTRINSIC abs, max, min, real, 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, 'SGEQRF', ' ', nmax, nmax, -1, -1 ),
521 $ ilaenv( 1, 'SORMQR', 'LT', nmax, nmax, nmax, -1 ),
522 $ ilaenv( 1, 'SORGQR', ' ', 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( 'SDRGES3', -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 = slamch( 'Safe minimum' )
541 ulp = slamch( 'Epsilon' )*slamch( 'Base' )
542 safmin = safmin / ulp
543 safmax = one / safmin
544 CALL slabad( safmin, safmax )
545 ulpinv = one / ulp
546*
547* The values RMAGN(2:3) depend on N, see below.
548*
549 rmagn( 0 ) = zero
550 rmagn( 1 ) = one
551*
552* Loop over matrix sizes
553*
554 ntestt = 0
555 nerrs = 0
556 nmats = 0
557*
558 DO 190 jsize = 1, nsizes
559 n = nn( jsize )
560 n1 = max( 1, n )
561 rmagn( 2 ) = safmax*ulp / real( n1 )
562 rmagn( 3 ) = safmin*ulpinv*real( n1 )
563*
564 IF( nsizes.NE.1 ) THEN
565 mtypes = min( maxtyp, ntypes )
566 ELSE
567 mtypes = min( maxtyp+1, ntypes )
568 END IF
569*
570* Loop over matrix types
571*
572 DO 180 jtype = 1, mtypes
573 IF( .NOT.dotype( jtype ) )
574 $ GO TO 180
575 nmats = nmats + 1
576 ntest = 0
577*
578* Save ISEED in case of an error.
579*
580 DO 20 j = 1, 4
581 ioldsd( j ) = iseed( j )
582 20 CONTINUE
583*
584* Initialize RESULT
585*
586 DO 30 j = 1, 13
587 result( j ) = zero
588 30 CONTINUE
589*
590* Generate test matrices A and B
591*
592* Description of control parameters:
593*
594* KCLASS: =1 means w/o rotation, =2 means w/ rotation,
595* =3 means random.
596* KATYPE: the "type" to be passed to SLATM4 for computing A.
597* KAZERO: the pattern of zeros on the diagonal for A:
598* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
599* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
600* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
601* non-zero entries.)
602* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
603* =2: large, =3: small.
604* IASIGN: 1 if the diagonal elements of A are to be
605* multiplied by a random magnitude 1 number, =2 if
606* randomly chosen diagonal blocks are to be rotated
607* to form 2x2 blocks.
608* KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
609* KTRIAN: =0: don't fill in the upper triangle, =1: do.
610* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
611* RMAGN: used to implement KAMAGN and KBMAGN.
612*
613 IF( mtypes.GT.maxtyp )
614 $ GO TO 110
615 iinfo = 0
616 IF( kclass( jtype ).LT.3 ) THEN
617*
618* Generate A (w/o rotation)
619*
620 IF( abs( katype( jtype ) ).EQ.3 ) THEN
621 in = 2*( ( n-1 ) / 2 ) + 1
622 IF( in.NE.n )
623 $ CALL slaset( 'Full', n, n, zero, zero, a, lda )
624 ELSE
625 in = n
626 END IF
627 CALL slatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
628 $ kz2( kazero( jtype ) ), iasign( jtype ),
629 $ rmagn( kamagn( jtype ) ), ulp,
630 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
631 $ iseed, a, lda )
632 iadd = kadd( kazero( jtype ) )
633 IF( iadd.GT.0 .AND. iadd.LE.n )
634 $ a( iadd, iadd ) = one
635*
636* Generate B (w/o rotation)
637*
638 IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
639 in = 2*( ( n-1 ) / 2 ) + 1
640 IF( in.NE.n )
641 $ CALL slaset( 'Full', n, n, zero, zero, b, lda )
642 ELSE
643 in = n
644 END IF
645 CALL slatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
646 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
647 $ rmagn( kbmagn( jtype ) ), one,
648 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
649 $ iseed, b, lda )
650 iadd = kadd( kbzero( jtype ) )
651 IF( iadd.NE.0 .AND. iadd.LE.n )
652 $ b( iadd, iadd ) = one
653*
654 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
655*
656* Include rotations
657*
658* Generate Q, Z as Householder transformations times
659* a diagonal matrix.
660*
661 DO 50 jc = 1, n - 1
662 DO 40 jr = jc, n
663 q( jr, jc ) = slarnd( 3, iseed )
664 z( jr, jc ) = slarnd( 3, iseed )
665 40 CONTINUE
666 CALL slarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
667 $ work( jc ) )
668 work( 2*n+jc ) = sign( one, q( jc, jc ) )
669 q( jc, jc ) = one
670 CALL slarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
671 $ work( n+jc ) )
672 work( 3*n+jc ) = sign( one, z( jc, jc ) )
673 z( jc, jc ) = one
674 50 CONTINUE
675 q( n, n ) = one
676 work( n ) = zero
677 work( 3*n ) = sign( one, slarnd( 2, iseed ) )
678 z( n, n ) = one
679 work( 2*n ) = zero
680 work( 4*n ) = sign( one, slarnd( 2, iseed ) )
681*
682* Apply the diagonal matrices
683*
684 DO 70 jc = 1, n
685 DO 60 jr = 1, n
686 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
687 $ a( jr, jc )
688 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
689 $ b( jr, jc )
690 60 CONTINUE
691 70 CONTINUE
692 CALL sorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
693 $ lda, work( 2*n+1 ), iinfo )
694 IF( iinfo.NE.0 )
695 $ GO TO 100
696 CALL sorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
697 $ a, lda, work( 2*n+1 ), iinfo )
698 IF( iinfo.NE.0 )
699 $ GO TO 100
700 CALL sorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
701 $ lda, work( 2*n+1 ), iinfo )
702 IF( iinfo.NE.0 )
703 $ GO TO 100
704 CALL sorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
705 $ b, lda, work( 2*n+1 ), iinfo )
706 IF( iinfo.NE.0 )
707 $ GO TO 100
708 END IF
709 ELSE
710*
711* Random matrices
712*
713 DO 90 jc = 1, n
714 DO 80 jr = 1, n
715 a( jr, jc ) = rmagn( kamagn( jtype ) )*
716 $ slarnd( 2, iseed )
717 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
718 $ slarnd( 2, iseed )
719 80 CONTINUE
720 90 CONTINUE
721 END IF
722*
723 100 CONTINUE
724*
725 IF( iinfo.NE.0 ) THEN
726 WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
727 $ ioldsd
728 info = abs( iinfo )
729 RETURN
730 END IF
731*
732 110 CONTINUE
733*
734 DO 120 i = 1, 13
735 result( i ) = -one
736 120 CONTINUE
737*
738* Test with and without sorting of eigenvalues
739*
740 DO 150 isort = 0, 1
741 IF( isort.EQ.0 ) THEN
742 sort = 'N'
743 rsub = 0
744 ELSE
745 sort = 'S'
746 rsub = 5
747 END IF
748*
749* Call XLAENV to set the parameters used in SLAQZ0
750*
751 CALL xlaenv( 12, 10 )
752 CALL xlaenv( 13, 12 )
753 CALL xlaenv( 14, 13 )
754 CALL xlaenv( 15, 2 )
755 CALL xlaenv( 17, 10 )
756*
757* Call SGGES3 to compute H, T, Q, Z, alpha, and beta.
758*
759 CALL slacpy( 'Full', n, n, a, lda, s, lda )
760 CALL slacpy( 'Full', n, n, b, lda, t, lda )
761 ntest = 1 + rsub + isort
762 result( 1+rsub+isort ) = ulpinv
763 CALL sgges3( 'V', 'V', sort, slctes, n, s, lda, t, lda,
764 $ sdim, alphar, alphai, beta, q, ldq, z, ldq,
765 $ work, lwork, bwork, iinfo )
766 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
767 result( 1+rsub+isort ) = ulpinv
768 WRITE( nounit, fmt = 9999 )'SGGES3', iinfo, n, jtype,
769 $ ioldsd
770 info = abs( iinfo )
771 GO TO 160
772 END IF
773*
774 ntest = 4 + rsub
775*
776* Do tests 1--4 (or tests 7--9 when reordering )
777*
778 IF( isort.EQ.0 ) THEN
779 CALL sget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
780 $ work, result( 1 ) )
781 CALL sget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
782 $ work, result( 2 ) )
783 ELSE
784 CALL sget54( n, a, lda, b, lda, s, lda, t, lda, q,
785 $ ldq, z, ldq, work, result( 7 ) )
786 END IF
787 CALL sget51( 3, n, a, lda, t, lda, q, ldq, q, ldq, work,
788 $ result( 3+rsub ) )
789 CALL sget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
790 $ result( 4+rsub ) )
791*
792* Do test 5 and 6 (or Tests 10 and 11 when reordering):
793* check Schur form of A and compare eigenvalues with
794* diagonals.
795*
796 ntest = 6 + rsub
797 temp1 = zero
798*
799 DO 130 j = 1, n
800 ilabad = .false.
801 IF( alphai( j ).EQ.zero ) THEN
802 temp2 = ( abs( alphar( j )-s( j, j ) ) /
803 $ max( safmin, abs( alphar( j ) ), abs( s( j,
804 $ j ) ) )+abs( beta( j )-t( j, j ) ) /
805 $ max( safmin, abs( beta( j ) ), abs( t( j,
806 $ j ) ) ) ) / ulp
807*
808 IF( j.LT.n ) THEN
809 IF( s( j+1, j ).NE.zero ) THEN
810 ilabad = .true.
811 result( 5+rsub ) = ulpinv
812 END IF
813 END IF
814 IF( j.GT.1 ) THEN
815 IF( s( j, j-1 ).NE.zero ) THEN
816 ilabad = .true.
817 result( 5+rsub ) = ulpinv
818 END IF
819 END IF
820*
821 ELSE
822 IF( alphai( j ).GT.zero ) THEN
823 i1 = j
824 ELSE
825 i1 = j - 1
826 END IF
827 IF( i1.LE.0 .OR. i1.GE.n ) THEN
828 ilabad = .true.
829 ELSE IF( i1.LT.n-1 ) THEN
830 IF( s( i1+2, i1+1 ).NE.zero ) THEN
831 ilabad = .true.
832 result( 5+rsub ) = ulpinv
833 END IF
834 ELSE IF( i1.GT.1 ) THEN
835 IF( s( i1, i1-1 ).NE.zero ) THEN
836 ilabad = .true.
837 result( 5+rsub ) = ulpinv
838 END IF
839 END IF
840 IF( .NOT.ilabad ) THEN
841 CALL sget53( s( i1, i1 ), lda, t( i1, i1 ), lda,
842 $ beta( j ), alphar( j ),
843 $ alphai( j ), temp2, ierr )
844 IF( ierr.GE.3 ) THEN
845 WRITE( nounit, fmt = 9998 )ierr, j, n,
846 $ jtype, ioldsd
847 info = abs( ierr )
848 END IF
849 ELSE
850 temp2 = ulpinv
851 END IF
852*
853 END IF
854 temp1 = max( temp1, temp2 )
855 IF( ilabad ) THEN
856 WRITE( nounit, fmt = 9997 )j, n, jtype, ioldsd
857 END IF
858 130 CONTINUE
859 result( 6+rsub ) = temp1
860*
861 IF( isort.GE.1 ) THEN
862*
863* Do test 12
864*
865 ntest = 12
866 result( 12 ) = zero
867 knteig = 0
868 DO 140 i = 1, n
869 IF( slctes( alphar( i ), alphai( i ),
870 $ beta( i ) ) .OR. slctes( alphar( i ),
871 $ -alphai( i ), beta( i ) ) ) THEN
872 knteig = knteig + 1
873 END IF
874 IF( i.LT.n ) THEN
875 IF( ( slctes( alphar( i+1 ), alphai( i+1 ),
876 $ beta( i+1 ) ) .OR. slctes( alphar( i+1 ),
877 $ -alphai( i+1 ), beta( i+1 ) ) ) .AND.
878 $ ( .NOT.( slctes( alphar( i ), alphai( i ),
879 $ beta( i ) ) .OR. slctes( alphar( i ),
880 $ -alphai( i ), beta( i ) ) ) ) .AND.
881 $ iinfo.NE.n+2 ) THEN
882 result( 12 ) = ulpinv
883 END IF
884 END IF
885 140 CONTINUE
886 IF( sdim.NE.knteig ) THEN
887 result( 12 ) = ulpinv
888 END IF
889 END IF
890*
891 150 CONTINUE
892*
893* End of Loop -- Check for RESULT(j) > THRESH
894*
895 160 CONTINUE
896*
897 ntestt = ntestt + ntest
898*
899* Print out tests which fail.
900*
901 DO 170 jr = 1, ntest
902 IF( result( jr ).GE.thresh ) THEN
903*
904* If this is the first test to fail,
905* print a header to the data file.
906*
907 IF( nerrs.EQ.0 ) THEN
908 WRITE( nounit, fmt = 9996 )'SGS'
909*
910* Matrix types
911*
912 WRITE( nounit, fmt = 9995 )
913 WRITE( nounit, fmt = 9994 )
914 WRITE( nounit, fmt = 9993 )'Orthogonal'
915*
916* Tests performed
917*
918 WRITE( nounit, fmt = 9992 )'orthogonal', '''',
919 $ 'transpose', ( '''', j = 1, 8 )
920*
921 END IF
922 nerrs = nerrs + 1
923 IF( result( jr ).LT.10000.0 ) THEN
924 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
925 $ result( jr )
926 ELSE
927 WRITE( nounit, fmt = 9990 )n, jtype, ioldsd, jr,
928 $ result( jr )
929 END IF
930 END IF
931 170 CONTINUE
932*
933 180 CONTINUE
934 190 CONTINUE
935*
936* Summary
937*
938 CALL alasvm( 'SGS', nounit, nerrs, ntestt, 0 )
939*
940 work( 1 ) = maxwrk
941*
942 RETURN
943*
944 9999 FORMAT( ' SDRGES3: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
945 $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
946*
947 9998 FORMAT( ' SDRGES3: SGET53 returned INFO=', i1, ' for eigenvalue ',
948 $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(',
949 $ 4( i4, ',' ), i5, ')' )
950*
951 9997 FORMAT( ' SDRGES3: S not in Schur form at eigenvalue ', i6, '.',
952 $ / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
953 $ i5, ')' )
954*
955 9996 FORMAT( / 1x, a3, ' -- Real Generalized Schur form driver' )
956*
957 9995 FORMAT( ' Matrix types (see SDRGES3 for details): ' )
958*
959 9994 FORMAT( ' Special Matrices:', 23x,
960 $ '(J''=transposed Jordan block)',
961 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
962 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
963 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
964 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
965 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
966 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
967 9993 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
968 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
969 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
970 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
971 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
972 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
973 $ '23=(small,large) 24=(small,small) 25=(large,large)',
974 $ / ' 26=random O(1) matrices.' )
975*
976 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
977 $ 'Q and Z are ', a, ',', / 19x,
978 $ 'l and r are the appropriate left and right', / 19x,
979 $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
980 $ ' means ', a, '.)', / ' Without ordering: ',
981 $ / ' 1 = | A - Q S Z', a,
982 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
983 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
984 $ ' | / ( n ulp ) 4 = | I - ZZ', a,
985 $ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
986 $ / ' 6 = difference between (alpha,beta)',
987 $ ' and diagonals of (S,T)', / ' With ordering: ',
988 $ / ' 7 = | (A,B) - Q (S,T) Z', a,
989 $ ' | / ( |(A,B)| n ulp ) ', / ' 8 = | I - QQ', a,
990 $ ' | / ( n ulp ) 9 = | I - ZZ', a,
991 $ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
992 $ / ' 11 = difference between (alpha,beta) and diagonals',
993 $ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
994 $ 'selected eigenvalues', / )
995 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
996 $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
997 9990 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
998 $ 4( i4, ',' ), ' result ', i2, ' is', 1p, e10.3 )
999*
1000* End of SDRGES3
1001*
1002 END
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
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 sgges3(JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, BWORK, INFO)
SGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: sgges3.f:282
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:106
subroutine sorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: sorm2r.f:159
subroutine sget54(N, A, LDA, B, LDB, S, LDS, T, LDT, U, LDU, V, LDV, WORK, RESULT)
SGET54
Definition: sget54.f:156
subroutine sget53(A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO)
SGET53
Definition: sget53.f:126
subroutine sdrges3(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHAR, ALPHAI, BETA, WORK, LWORK, RESULT, BWORK, INFO)
SDRGES3
Definition: sdrges3.f:403
subroutine sget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
SGET51
Definition: sget51.f:149
subroutine slatm4(ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
SLATM4
Definition: slatm4.f:175