LAPACK 3.12.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, slacpy,
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 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 / real( n1 )
561 rmagn( 3 ) = safmin*ulpinv*real( 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* KCLASS: =1 means w/o rotation, =2 means w/ rotation,
594* =3 means random.
595* KATYPE: the "type" to be passed to SLATM4 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 slaset( 'Full', n, n, zero, zero, a, lda )
623 ELSE
624 in = n
625 END IF
626 CALL slatm4( 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 slaset( 'Full', n, n, zero, zero, b, lda )
641 ELSE
642 in = n
643 END IF
644 CALL slatm4( 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 ) = slarnd( 3, iseed )
663 z( jr, jc ) = slarnd( 3, iseed )
664 40 CONTINUE
665 CALL slarfg( 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 slarfg( 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, slarnd( 2, iseed ) )
677 z( n, n ) = one
678 work( 2*n ) = zero
679 work( 4*n ) = sign( one, slarnd( 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 sorm2r( '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 sorm2r( '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 sorm2r( '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 sorm2r( '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 $ slarnd( 2, iseed )
716 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
717 $ slarnd( 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 XLAENV to set the parameters used in SLAQZ0
749*
750 CALL xlaenv( 12, 10 )
751 CALL xlaenv( 13, 12 )
752 CALL xlaenv( 14, 13 )
753 CALL xlaenv( 15, 2 )
754 CALL xlaenv( 17, 10 )
755*
756* Call SGGES3 to compute H, T, Q, Z, alpha, and beta.
757*
758 CALL slacpy( 'Full', n, n, a, lda, s, lda )
759 CALL slacpy( 'Full', n, n, b, lda, t, lda )
760 ntest = 1 + rsub + isort
761 result( 1+rsub+isort ) = ulpinv
762 CALL sgges3( 'V', 'V', sort, slctes, n, s, lda, t, lda,
763 $ sdim, alphar, alphai, beta, q, ldq, z, ldq,
764 $ work, lwork, bwork, iinfo )
765 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
766 result( 1+rsub+isort ) = ulpinv
767 WRITE( nounit, fmt = 9999 )'SGGES3', iinfo, n, jtype,
768 $ ioldsd
769 info = abs( iinfo )
770 GO TO 160
771 END IF
772*
773 ntest = 4 + rsub
774*
775* Do tests 1--4 (or tests 7--9 when reordering )
776*
777 IF( isort.EQ.0 ) THEN
778 CALL sget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
779 $ work, result( 1 ) )
780 CALL sget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
781 $ work, result( 2 ) )
782 ELSE
783 CALL sget54( n, a, lda, b, lda, s, lda, t, lda, q,
784 $ ldq, z, ldq, work, result( 7 ) )
785 END IF
786 CALL sget51( 3, n, a, lda, t, lda, q, ldq, q, ldq, work,
787 $ result( 3+rsub ) )
788 CALL sget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
789 $ result( 4+rsub ) )
790*
791* Do test 5 and 6 (or Tests 10 and 11 when reordering):
792* check Schur form of A and compare eigenvalues with
793* diagonals.
794*
795 ntest = 6 + rsub
796 temp1 = zero
797*
798 DO 130 j = 1, n
799 ilabad = .false.
800 IF( alphai( j ).EQ.zero ) THEN
801 temp2 = ( abs( alphar( j )-s( j, j ) ) /
802 $ max( safmin, abs( alphar( j ) ), abs( s( j,
803 $ j ) ) )+abs( beta( j )-t( j, j ) ) /
804 $ max( safmin, abs( beta( j ) ), abs( t( j,
805 $ j ) ) ) ) / ulp
806*
807 IF( j.LT.n ) THEN
808 IF( s( j+1, j ).NE.zero ) THEN
809 ilabad = .true.
810 result( 5+rsub ) = ulpinv
811 END IF
812 END IF
813 IF( j.GT.1 ) THEN
814 IF( s( j, j-1 ).NE.zero ) THEN
815 ilabad = .true.
816 result( 5+rsub ) = ulpinv
817 END IF
818 END IF
819*
820 ELSE
821 IF( alphai( j ).GT.zero ) THEN
822 i1 = j
823 ELSE
824 i1 = j - 1
825 END IF
826 IF( i1.LE.0 .OR. i1.GE.n ) THEN
827 ilabad = .true.
828 ELSE IF( i1.LT.n-1 ) THEN
829 IF( s( i1+2, i1+1 ).NE.zero ) THEN
830 ilabad = .true.
831 result( 5+rsub ) = ulpinv
832 END IF
833 ELSE IF( i1.GT.1 ) THEN
834 IF( s( i1, i1-1 ).NE.zero ) THEN
835 ilabad = .true.
836 result( 5+rsub ) = ulpinv
837 END IF
838 END IF
839 IF( .NOT.ilabad ) THEN
840 CALL sget53( s( i1, i1 ), lda, t( i1, i1 ), lda,
841 $ beta( j ), alphar( j ),
842 $ alphai( j ), temp2, ierr )
843 IF( ierr.GE.3 ) THEN
844 WRITE( nounit, fmt = 9998 )ierr, j, n,
845 $ jtype, ioldsd
846 info = abs( ierr )
847 END IF
848 ELSE
849 temp2 = ulpinv
850 END IF
851*
852 END IF
853 temp1 = max( temp1, temp2 )
854 IF( ilabad ) THEN
855 WRITE( nounit, fmt = 9997 )j, n, jtype, ioldsd
856 END IF
857 130 CONTINUE
858 result( 6+rsub ) = temp1
859*
860 IF( isort.GE.1 ) THEN
861*
862* Do test 12
863*
864 ntest = 12
865 result( 12 ) = zero
866 knteig = 0
867 DO 140 i = 1, n
868 IF( slctes( alphar( i ), alphai( i ),
869 $ beta( i ) ) .OR. slctes( alphar( i ),
870 $ -alphai( i ), beta( i ) ) ) THEN
871 knteig = knteig + 1
872 END IF
873 IF( i.LT.n ) THEN
874 IF( ( slctes( alphar( i+1 ), alphai( i+1 ),
875 $ beta( i+1 ) ) .OR. slctes( alphar( i+1 ),
876 $ -alphai( i+1 ), beta( i+1 ) ) ) .AND.
877 $ ( .NOT.( slctes( alphar( i ), alphai( i ),
878 $ beta( i ) ) .OR. slctes( alphar( i ),
879 $ -alphai( i ), beta( i ) ) ) ) .AND.
880 $ iinfo.NE.n+2 ) THEN
881 result( 12 ) = ulpinv
882 END IF
883 END IF
884 140 CONTINUE
885 IF( sdim.NE.knteig ) THEN
886 result( 12 ) = ulpinv
887 END IF
888 END IF
889*
890 150 CONTINUE
891*
892* End of Loop -- Check for RESULT(j) > THRESH
893*
894 160 CONTINUE
895*
896 ntestt = ntestt + ntest
897*
898* Print out tests which fail.
899*
900 DO 170 jr = 1, ntest
901 IF( result( jr ).GE.thresh ) THEN
902*
903* If this is the first test to fail,
904* print a header to the data file.
905*
906 IF( nerrs.EQ.0 ) THEN
907 WRITE( nounit, fmt = 9996 )'SGS'
908*
909* Matrix types
910*
911 WRITE( nounit, fmt = 9995 )
912 WRITE( nounit, fmt = 9994 )
913 WRITE( nounit, fmt = 9993 )'Orthogonal'
914*
915* Tests performed
916*
917 WRITE( nounit, fmt = 9992 )'orthogonal', '''',
918 $ 'transpose', ( '''', j = 1, 8 )
919*
920 END IF
921 nerrs = nerrs + 1
922 IF( result( jr ).LT.10000.0 ) THEN
923 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
924 $ result( jr )
925 ELSE
926 WRITE( nounit, fmt = 9990 )n, jtype, ioldsd, jr,
927 $ result( jr )
928 END IF
929 END IF
930 170 CONTINUE
931*
932 180 CONTINUE
933 190 CONTINUE
934*
935* Summary
936*
937 CALL alasvm( 'SGS', nounit, nerrs, ntestt, 0 )
938*
939 work( 1 ) = maxwrk
940*
941 RETURN
942*
943 9999 FORMAT( ' SDRGES3: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
944 $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
945*
946 9998 FORMAT( ' SDRGES3: SGET53 returned INFO=', i1, ' for eigenvalue ',
947 $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(',
948 $ 4( i4, ',' ), i5, ')' )
949*
950 9997 FORMAT( ' SDRGES3: S not in Schur form at eigenvalue ', i6, '.',
951 $ / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
952 $ i5, ')' )
953*
954 9996 FORMAT( / 1x, a3, ' -- Real Generalized Schur form driver' )
955*
956 9995 FORMAT( ' Matrix types (see SDRGES3 for details): ' )
957*
958 9994 FORMAT( ' Special Matrices:', 23x,
959 $ '(J''=transposed Jordan block)',
960 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
961 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
962 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
963 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
964 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
965 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
966 9993 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
967 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
968 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
969 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
970 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
971 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
972 $ '23=(small,large) 24=(small,small) 25=(large,large)',
973 $ / ' 26=random O(1) matrices.' )
974*
975 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
976 $ 'Q and Z are ', a, ',', / 19x,
977 $ 'l and r are the appropriate left and right', / 19x,
978 $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
979 $ ' means ', a, '.)', / ' Without ordering: ',
980 $ / ' 1 = | A - Q S Z', a,
981 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
982 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
983 $ ' | / ( n ulp ) 4 = | I - ZZ', a,
984 $ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
985 $ / ' 6 = difference between (alpha,beta)',
986 $ ' and diagonals of (S,T)', / ' With ordering: ',
987 $ / ' 7 = | (A,B) - Q (S,T) Z', a,
988 $ ' | / ( |(A,B)| n ulp ) ', / ' 8 = | I - QQ', a,
989 $ ' | / ( n ulp ) 9 = | I - ZZ', a,
990 $ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
991 $ / ' 11 = difference between (alpha,beta) and diagonals',
992 $ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
993 $ 'selected eigenvalues', / )
994 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
995 $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
996 9990 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
997 $ 4( i4, ',' ), ' result ', i2, ' is', 1p, e10.3 )
998*
999* End of SDRGES3
1000*
1001 END
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine xlaenv(ispec, nvalue)
XLAENV
Definition xlaenv.f:81
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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 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 slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:106
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 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 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 sget53(a, lda, b, ldb, scale, wr, wi, result, info)
SGET53
Definition sget53.f:126
subroutine sget54(n, a, lda, b, ldb, s, lds, t, ldt, u, ldu, v, ldv, work, result)
SGET54
Definition sget54.f:156
subroutine slatm4(itype, n, nz1, nz2, isign, amagn, rcond, triang, idist, iseed, a, lda)
SLATM4
Definition slatm4.f:175