LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
sdrgev.f
Go to the documentation of this file.
1*> \brief \b SDRGEV
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 SDRGEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12* NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE,
13* ALPHAR, ALPHAI, BETA, ALPHR1, ALPHI1, BETA1,
14* WORK, LWORK, RESULT, INFO )
15*
16* .. Scalar Arguments ..
17* INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
18* \$ NTYPES
19* REAL THRESH
20* ..
21* .. Array Arguments ..
22* LOGICAL DOTYPE( * )
23* INTEGER ISEED( 4 ), NN( * )
24* REAL A( LDA, * ), ALPHAI( * ), ALPHI1( * ),
25* \$ ALPHAR( * ), ALPHR1( * ), B( LDA, * ),
26* \$ BETA( * ), BETA1( * ), Q( LDQ, * ),
27* \$ QE( LDQE, * ), RESULT( * ), S( LDA, * ),
28* \$ T( LDA, * ), WORK( * ), Z( LDQ, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SDRGEV checks the nonsymmetric generalized eigenvalue problem driver
38*> routine SGGEV.
39*>
40*> SGGEV computes for a pair of n-by-n nonsymmetric matrices (A,B) the
41*> generalized eigenvalues and, optionally, the left and right
42*> eigenvectors.
43*>
44*> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
45*> or a ratio alpha/beta = w, such that A - w*B is singular. It is
46*> usually represented as the pair (alpha,beta), as there is reasonable
47*> interpretation for beta=0, and even for both being zero.
48*>
49*> A right generalized eigenvector corresponding to a generalized
50*> eigenvalue w for a pair of matrices (A,B) is a vector r such that
51*> (A - wB) * r = 0. A left generalized eigenvector is a vector l such
52*> that l**H * (A - wB) = 0, where l**H is the conjugate-transpose of l.
53*>
54*> When SDRGEV is called, a number of matrix "sizes" ("n's") and a
55*> number of matrix "types" are specified. For each size ("n")
56*> and each type of matrix, a pair of matrices (A, B) will be generated
57*> and used for testing. For each matrix pair, the following tests
58*> will be performed and compared with the threshold THRESH.
59*>
60*> Results from SGGEV:
61*>
62*> (1) max over all left eigenvalue/-vector pairs (alpha/beta,l) of
63*>
64*> | VL**H * (beta A - alpha B) |/( ulp max(|beta A|, |alpha B|) )
65*>
66*> where VL**H is the conjugate-transpose of VL.
67*>
68*> (2) | |VL(i)| - 1 | / ulp and whether largest component real
69*>
70*> VL(i) denotes the i-th column of VL.
71*>
72*> (3) max over all left eigenvalue/-vector pairs (alpha/beta,r) of
73*>
74*> | (beta A - alpha B) * VR | / ( ulp max(|beta A|, |alpha B|) )
75*>
76*> (4) | |VR(i)| - 1 | / ulp and whether largest component real
77*>
78*> VR(i) denotes the i-th column of VR.
79*>
80*> (5) W(full) = W(partial)
81*> W(full) denotes the eigenvalues computed when both l and r
82*> are also computed, and W(partial) denotes the eigenvalues
83*> computed when only W, only W and r, or only W and l are
84*> computed.
85*>
86*> (6) VL(full) = VL(partial)
87*> VL(full) denotes the left eigenvectors computed when both l
88*> and r are computed, and VL(partial) denotes the result
89*> when only l is computed.
90*>
91*> (7) VR(full) = VR(partial)
92*> VR(full) denotes the right eigenvectors computed when both l
93*> and r are also computed, and VR(partial) denotes the result
94*> when only l is computed.
95*>
96*>
97*> Test Matrices
98*> ---- --------
99*>
100*> The sizes of the test matrices are specified by an array
101*> NN(1:NSIZES); the value of each element NN(j) specifies one size.
102*> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
103*> DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
104*> Currently, the list of possible types is:
105*>
106*> (1) ( 0, 0 ) (a pair of zero matrices)
107*>
108*> (2) ( I, 0 ) (an identity and a zero matrix)
109*>
110*> (3) ( 0, I ) (an identity and a zero matrix)
111*>
112*> (4) ( I, I ) (a pair of identity matrices)
113*>
114*> t t
115*> (5) ( J , J ) (a pair of transposed Jordan blocks)
116*>
117*> t ( I 0 )
118*> (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
119*> ( 0 I ) ( 0 J )
120*> and I is a k x k identity and J a (k+1)x(k+1)
121*> Jordan block; k=(N-1)/2
122*>
123*> (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal
124*> matrix with those diagonal entries.)
125*> (8) ( I, D )
126*>
127*> (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
128*>
129*> (10) ( small*D, big*I )
130*>
131*> (11) ( big*I, small*D )
132*>
133*> (12) ( small*I, big*D )
134*>
135*> (13) ( big*D, big*I )
136*>
137*> (14) ( small*D, small*I )
138*>
139*> (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
140*> D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
141*> t t
142*> (16) Q ( J , J ) Z where Q and Z are random orthogonal matrices.
143*>
144*> (17) Q ( T1, T2 ) Z where T1 and T2 are upper triangular matrices
145*> with random O(1) entries above the diagonal
146*> and diagonal entries diag(T1) =
147*> ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
148*> ( 0, N-3, N-4,..., 1, 0, 0 )
149*>
150*> (18) Q ( T1, T2 ) Z diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
151*> diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
152*> s = machine precision.
153*>
154*> (19) Q ( T1, T2 ) Z diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
155*> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
156*>
157*> N-5
158*> (20) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
159*> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
160*>
161*> (21) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
162*> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
163*> where r1,..., r(N-4) are random.
164*>
165*> (22) Q ( big*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
166*> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
167*>
168*> (23) Q ( small*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
169*> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
170*>
171*> (24) Q ( small*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
172*> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
173*>
174*> (25) Q ( big*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
175*> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
176*>
177*> (26) Q ( T1, T2 ) Z where T1 and T2 are random upper-triangular
178*> matrices.
179*>
180*> \endverbatim
181*
182* Arguments:
183* ==========
184*
185*> \param[in] NSIZES
186*> \verbatim
187*> NSIZES is INTEGER
188*> The number of sizes of matrices to use. If it is zero,
189*> SDRGES does nothing. NSIZES >= 0.
190*> \endverbatim
191*>
192*> \param[in] NN
193*> \verbatim
194*> NN is INTEGER array, dimension (NSIZES)
195*> An array containing the sizes to be used for the matrices.
196*> Zero values will be skipped. NN >= 0.
197*> \endverbatim
198*>
199*> \param[in] NTYPES
200*> \verbatim
201*> NTYPES is INTEGER
202*> The number of elements in DOTYPE. If it is zero, SDRGES
203*> does nothing. It must be at least zero. If it is MAXTYP+1
204*> and NSIZES is 1, then an additional type, MAXTYP+1 is
205*> defined, which is to use whatever matrix is in A. This
206*> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
207*> DOTYPE(MAXTYP+1) is .TRUE. .
208*> \endverbatim
209*>
210*> \param[in] DOTYPE
211*> \verbatim
212*> DOTYPE is LOGICAL array, dimension (NTYPES)
213*> If DOTYPE(j) is .TRUE., then for each size in NN a
214*> matrix of that size and of type j will be generated.
215*> If NTYPES is smaller than the maximum number of types
216*> defined (PARAMETER MAXTYP), then types NTYPES+1 through
217*> MAXTYP will not be generated. If NTYPES is larger
218*> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
219*> will be ignored.
220*> \endverbatim
221*>
222*> \param[in,out] ISEED
223*> \verbatim
224*> ISEED is INTEGER array, dimension (4)
225*> On entry ISEED specifies the seed of the random number
226*> generator. The array elements should be between 0 and 4095;
227*> if not they will be reduced mod 4096. Also, ISEED(4) must
228*> be odd. The random number generator uses a linear
229*> congruential sequence limited to small integers, and so
230*> should produce machine independent random numbers. The
231*> values of ISEED are changed on exit, and can be used in the
232*> next call to SDRGES to continue the same random number
233*> sequence.
234*> \endverbatim
235*>
236*> \param[in] THRESH
237*> \verbatim
238*> THRESH is REAL
239*> A test will count as "failed" if the "error", computed as
240*> described above, exceeds THRESH. Note that the error is
241*> scaled to be O(1), so THRESH should be a reasonably small
242*> multiple of 1, e.g., 10 or 100. In particular, it should
243*> not depend on the precision (single vs. double) or the size
244*> of the matrix. It must be at least zero.
245*> \endverbatim
246*>
247*> \param[in] NOUNIT
248*> \verbatim
249*> NOUNIT is INTEGER
250*> The FORTRAN unit number for printing out error messages
251*> (e.g., if a routine returns IERR not equal to 0.)
252*> \endverbatim
253*>
254*> \param[in,out] A
255*> \verbatim
256*> A is REAL array,
257*> dimension(LDA, max(NN))
258*> Used to hold the original A matrix. Used as input only
259*> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
260*> DOTYPE(MAXTYP+1)=.TRUE.
261*> \endverbatim
262*>
263*> \param[in] LDA
264*> \verbatim
265*> LDA is INTEGER
266*> The leading dimension of A, B, S, and T.
267*> It must be at least 1 and at least max( NN ).
268*> \endverbatim
269*>
270*> \param[in,out] B
271*> \verbatim
272*> B is REAL array,
273*> dimension(LDA, max(NN))
274*> Used to hold the original B matrix. Used as input only
275*> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
276*> DOTYPE(MAXTYP+1)=.TRUE.
277*> \endverbatim
278*>
279*> \param[out] S
280*> \verbatim
281*> S is REAL array,
282*> dimension (LDA, max(NN))
283*> The Schur form matrix computed from A by SGGES. On exit, S
284*> contains the Schur form matrix corresponding to the matrix
285*> in A.
286*> \endverbatim
287*>
288*> \param[out] T
289*> \verbatim
290*> T is REAL array,
291*> dimension (LDA, max(NN))
292*> The upper triangular matrix computed from B by SGGES.
293*> \endverbatim
294*>
295*> \param[out] Q
296*> \verbatim
297*> Q is REAL array,
298*> dimension (LDQ, max(NN))
299*> The (left) eigenvectors matrix computed by SGGEV.
300*> \endverbatim
301*>
302*> \param[in] LDQ
303*> \verbatim
304*> LDQ is INTEGER
305*> The leading dimension of Q and Z. It must
306*> be at least 1 and at least max( NN ).
307*> \endverbatim
308*>
309*> \param[out] Z
310*> \verbatim
311*> Z is REAL array, dimension( LDQ, max(NN) )
312*> The (right) orthogonal matrix computed by SGGES.
313*> \endverbatim
314*>
315*> \param[out] QE
316*> \verbatim
317*> QE is REAL array, dimension( LDQ, max(NN) )
318*> QE holds the computed right or left eigenvectors.
319*> \endverbatim
320*>
321*> \param[in] LDQE
322*> \verbatim
323*> LDQE is INTEGER
324*> The leading dimension of QE. LDQE >= max(1,max(NN)).
325*> \endverbatim
326*>
327*> \param[out] ALPHAR
328*> \verbatim
329*> ALPHAR is REAL array, dimension (max(NN))
330*> \endverbatim
331*>
332*> \param[out] ALPHAI
333*> \verbatim
334*> ALPHAI is REAL array, dimension (max(NN))
335*> \endverbatim
336*>
337*> \param[out] BETA
338*> \verbatim
339*> BETA is REAL array, dimension (max(NN))
340*> \verbatim
341*> The generalized eigenvalues of (A,B) computed by SGGEV.
342*> ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
343*> generalized eigenvalue of A and B.
344*> \endverbatim
345*>
346*> \param[out] ALPHR1
347*> \verbatim
348*> ALPHR1 is REAL array, dimension (max(NN))
349*> \endverbatim
350*>
351*> \param[out] ALPHI1
352*> \verbatim
353*> ALPHI1 is REAL array, dimension (max(NN))
354*> \endverbatim
355*>
356*> \param[out] BETA1
357*> \verbatim
358*> BETA1 is REAL array, dimension (max(NN))
359*>
360*> Like ALPHAR, ALPHAI, BETA, these arrays contain the
361*> eigenvalues of A and B, but those computed when SGGEV only
362*> computes a partial eigendecomposition, i.e. not the
363*> eigenvalues and left and right eigenvectors.
364*> \endverbatim
365*>
366*> \param[out] WORK
367*> \verbatim
368*> WORK is REAL array, dimension (LWORK)
369*> \endverbatim
370*>
371*> \param[in] LWORK
372*> \verbatim
373*> LWORK is INTEGER
374*> The number of entries in WORK. LWORK >= MAX( 8*N, N*(N+1) ).
375*> \endverbatim
376*>
377*> \param[out] RESULT
378*> \verbatim
379*> RESULT is REAL array, dimension (2)
380*> The values computed by the tests described above.
381*> The values are currently limited to 1/ulp, to avoid overflow.
382*> \endverbatim
383*>
384*> \param[out] INFO
385*> \verbatim
386*> INFO is INTEGER
387*> = 0: successful exit
388*> < 0: if INFO = -i, the i-th argument had an illegal value.
389*> > 0: A routine returned an error code. INFO is the
390*> absolute value of the INFO value returned.
391*> \endverbatim
392*
393* Authors:
394* ========
395*
396*> \author Univ. of Tennessee
397*> \author Univ. of California Berkeley
398*> \author Univ. of Colorado Denver
399*> \author NAG Ltd.
400*
401*> \ingroup single_eig
402*
403* =====================================================================
404 SUBROUTINE sdrgev( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
405 \$ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE,
406 \$ ALPHAR, ALPHAI, BETA, ALPHR1, ALPHI1, BETA1,
407 \$ WORK, LWORK, RESULT, INFO )
408*
409* -- LAPACK test routine --
410* -- LAPACK is a software package provided by Univ. of Tennessee, --
411* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
412*
413* .. Scalar Arguments ..
414 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
415 \$ NTYPES
416 REAL THRESH
417* ..
418* .. Array Arguments ..
419 LOGICAL DOTYPE( * )
420 INTEGER ISEED( 4 ), NN( * )
421 REAL A( LDA, * ), ALPHAI( * ), ALPHI1( * ),
422 \$ alphar( * ), alphr1( * ), b( lda, * ),
423 \$ beta( * ), beta1( * ), q( ldq, * ),
424 \$ qe( ldqe, * ), result( * ), s( lda, * ),
425 \$ t( lda, * ), work( * ), z( ldq, * )
426* ..
427*
428* =====================================================================
429*
430* .. Parameters ..
431 REAL ZERO, ONE
432 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
433 INTEGER MAXTYP
434 parameter( maxtyp = 26 )
435* ..
436* .. Local Scalars ..
438 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
439 \$ MAXWRK, MINWRK, MTYPES, N, N1, NERRS, NMATS,
440 \$ nmax, ntestt
441 REAL SAFMAX, SAFMIN, ULP, ULPINV
442* ..
443* .. Local Arrays ..
444 INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
445 \$ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
446 \$ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
447 \$ kbmagn( maxtyp ), kbtype( maxtyp ),
448 \$ kbzero( maxtyp ), kclass( maxtyp ),
449 \$ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
450 REAL RMAGN( 0: 3 )
451* ..
452* .. External Functions ..
453 INTEGER ILAENV
454 REAL SLAMCH, SLARND
455 EXTERNAL ILAENV, SLAMCH, SLARND
456* ..
457* .. External Subroutines ..
458 EXTERNAL alasvm, sget52, sggev, slacpy, slarfg,
460* ..
461* .. Intrinsic Functions ..
462 INTRINSIC abs, max, min, real, sign
463* ..
464* .. Data statements ..
465 DATA kclass / 15*1, 10*2, 1*3 /
466 DATA kz1 / 0, 1, 2, 1, 3, 3 /
467 DATA kz2 / 0, 0, 1, 2, 1, 1 /
468 DATA kadd / 0, 0, 0, 0, 3, 2 /
469 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
470 \$ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
471 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
472 \$ 1, 1, -4, 2, -4, 8*8, 0 /
473 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
474 \$ 4*5, 4*3, 1 /
475 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
476 \$ 4*6, 4*4, 1 /
477 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
478 \$ 2, 1 /
479 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
480 \$ 2, 1 /
481 DATA ktrian / 16*0, 10*1 /
482 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
483 \$ 5*2, 0 /
484 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
485* ..
486* .. Executable Statements ..
487*
488* Check for errors
489*
490 info = 0
491*
493 nmax = 1
494 DO 10 j = 1, nsizes
495 nmax = max( nmax, nn( j ) )
496 IF( nn( j ).LT.0 )
498 10 CONTINUE
499*
500 IF( nsizes.LT.0 ) THEN
501 info = -1
502 ELSE IF( badnn ) THEN
503 info = -2
504 ELSE IF( ntypes.LT.0 ) THEN
505 info = -3
506 ELSE IF( thresh.LT.zero ) THEN
507 info = -6
508 ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
509 info = -9
510 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
511 info = -14
512 ELSE IF( ldqe.LE.1 .OR. ldqe.LT.nmax ) THEN
513 info = -17
514 END IF
515*
516* Compute workspace
517* (Note: Comments in the code beginning "Workspace:" describe the
518* minimal amount of workspace needed at that point in the code,
519* as well as the preferred amount for good performance.
520* NB refers to the optimal block size for the immediately
521* following subroutine, as returned by ILAENV.
522*
523 minwrk = 1
524 IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
525 minwrk = max( 1, 8*nmax, nmax*( nmax+1 ) )
526 maxwrk = 7*nmax + nmax*ilaenv( 1, 'SGEQRF', ' ', nmax, 1, nmax,
527 \$ 0 )
528 maxwrk = max( maxwrk, nmax*( nmax+1 ) )
529 work( 1 ) = maxwrk
530 END IF
531*
532 IF( lwork.LT.minwrk )
533 \$ info = -25
534*
535 IF( info.NE.0 ) THEN
536 CALL xerbla( 'SDRGEV', -info )
537 RETURN
538 END IF
539*
540* Quick return if possible
541*
542 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
543 \$ RETURN
544*
545 safmin = slamch( 'Safe minimum' )
546 ulp = slamch( 'Epsilon' )*slamch( 'Base' )
547 safmin = safmin / ulp
548 safmax = one / safmin
549 ulpinv = one / ulp
550*
551* The values RMAGN(2:3) depend on N, see below.
552*
553 rmagn( 0 ) = zero
554 rmagn( 1 ) = one
555*
556* Loop over sizes, types
557*
558 ntestt = 0
559 nerrs = 0
560 nmats = 0
561*
562 DO 220 jsize = 1, nsizes
563 n = nn( jsize )
564 n1 = max( 1, n )
565 rmagn( 2 ) = safmax*ulp / real( n1 )
566 rmagn( 3 ) = safmin*ulpinv*n1
567*
568 IF( nsizes.NE.1 ) THEN
569 mtypes = min( maxtyp, ntypes )
570 ELSE
571 mtypes = min( maxtyp+1, ntypes )
572 END IF
573*
574 DO 210 jtype = 1, mtypes
575 IF( .NOT.dotype( jtype ) )
576 \$ GO TO 210
577 nmats = nmats + 1
578*
579* Save ISEED in case of an error.
580*
581 DO 20 j = 1, 4
582 ioldsd( j ) = iseed( j )
583 20 CONTINUE
584*
585* Generate test matrices A and B
586*
587* Description of control parameters:
588*
589* KCLASS: =1 means w/o rotation, =2 means w/ rotation,
590* =3 means random.
591* KATYPE: the "type" to be passed to SLATM4 for computing A.
592* KAZERO: the pattern of zeros on the diagonal for A:
593* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
594* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
595* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
596* non-zero entries.)
597* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
598* =2: large, =3: small.
599* IASIGN: 1 if the diagonal elements of A are to be
600* multiplied by a random magnitude 1 number, =2 if
601* randomly chosen diagonal blocks are to be rotated
602* to form 2x2 blocks.
603* KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
604* KTRIAN: =0: don't fill in the upper triangle, =1: do.
605* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
606* RMAGN: used to implement KAMAGN and KBMAGN.
607*
608 IF( mtypes.GT.maxtyp )
609 \$ GO TO 100
610 ierr = 0
611 IF( kclass( jtype ).LT.3 ) THEN
612*
613* Generate A (w/o rotation)
614*
615 IF( abs( katype( jtype ) ).EQ.3 ) THEN
616 in = 2*( ( n-1 ) / 2 ) + 1
617 IF( in.NE.n )
618 \$ CALL slaset( 'Full', n, n, zero, zero, a, lda )
619 ELSE
620 in = n
621 END IF
622 CALL slatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
623 \$ kz2( kazero( jtype ) ), iasign( jtype ),
624 \$ rmagn( kamagn( jtype ) ), ulp,
625 \$ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
626 \$ iseed, a, lda )
630*
631* Generate B (w/o rotation)
632*
633 IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
634 in = 2*( ( n-1 ) / 2 ) + 1
635 IF( in.NE.n )
636 \$ CALL slaset( 'Full', n, n, zero, zero, b, lda )
637 ELSE
638 in = n
639 END IF
640 CALL slatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
641 \$ kz2( kbzero( jtype ) ), ibsign( jtype ),
642 \$ rmagn( kbmagn( jtype ) ), one,
643 \$ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
644 \$ iseed, b, lda )
648*
649 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
650*
651* Include rotations
652*
653* Generate Q, Z as Householder transformations times
654* a diagonal matrix.
655*
656 DO 40 jc = 1, n - 1
657 DO 30 jr = jc, n
658 q( jr, jc ) = slarnd( 3, iseed )
659 z( jr, jc ) = slarnd( 3, iseed )
660 30 CONTINUE
661 CALL slarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
662 \$ work( jc ) )
663 work( 2*n+jc ) = sign( one, q( jc, jc ) )
664 q( jc, jc ) = one
665 CALL slarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
666 \$ work( n+jc ) )
667 work( 3*n+jc ) = sign( one, z( jc, jc ) )
668 z( jc, jc ) = one
669 40 CONTINUE
670 q( n, n ) = one
671 work( n ) = zero
672 work( 3*n ) = sign( one, slarnd( 2, iseed ) )
673 z( n, n ) = one
674 work( 2*n ) = zero
675 work( 4*n ) = sign( one, slarnd( 2, iseed ) )
676*
677* Apply the diagonal matrices
678*
679 DO 60 jc = 1, n
680 DO 50 jr = 1, n
681 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
682 \$ a( jr, jc )
683 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
684 \$ b( jr, jc )
685 50 CONTINUE
686 60 CONTINUE
687 CALL sorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
688 \$ lda, work( 2*n+1 ), ierr )
689 IF( ierr.NE.0 )
690 \$ GO TO 90
691 CALL sorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
692 \$ a, lda, work( 2*n+1 ), ierr )
693 IF( ierr.NE.0 )
694 \$ GO TO 90
695 CALL sorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
696 \$ lda, work( 2*n+1 ), ierr )
697 IF( ierr.NE.0 )
698 \$ GO TO 90
699 CALL sorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
700 \$ b, lda, work( 2*n+1 ), ierr )
701 IF( ierr.NE.0 )
702 \$ GO TO 90
703 END IF
704 ELSE
705*
706* Random matrices
707*
708 DO 80 jc = 1, n
709 DO 70 jr = 1, n
710 a( jr, jc ) = rmagn( kamagn( jtype ) )*
711 \$ slarnd( 2, iseed )
712 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
713 \$ slarnd( 2, iseed )
714 70 CONTINUE
715 80 CONTINUE
716 END IF
717*
718 90 CONTINUE
719*
720 IF( ierr.NE.0 ) THEN
721 WRITE( nounit, fmt = 9999 )'Generator', ierr, n, jtype,
722 \$ ioldsd
723 info = abs( ierr )
724 RETURN
725 END IF
726*
727 100 CONTINUE
728*
729 DO 110 i = 1, 7
730 result( i ) = -one
731 110 CONTINUE
732*
733* Call SGGEV to compute eigenvalues and eigenvectors.
734*
735 CALL slacpy( ' ', n, n, a, lda, s, lda )
736 CALL slacpy( ' ', n, n, b, lda, t, lda )
737 CALL sggev( 'V', 'V', n, s, lda, t, lda, alphar, alphai,
738 \$ beta, q, ldq, z, ldq, work, lwork, ierr )
739 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
740 result( 1 ) = ulpinv
741 WRITE( nounit, fmt = 9999 )'SGGEV1', ierr, n, jtype,
742 \$ ioldsd
743 info = abs( ierr )
744 GO TO 190
745 END IF
746*
747* Do the tests (1) and (2)
748*
749 CALL sget52( .true., n, a, lda, b, lda, q, ldq, alphar,
750 \$ alphai, beta, work, result( 1 ) )
751 IF( result( 2 ).GT.thresh ) THEN
752 WRITE( nounit, fmt = 9998 )'Left', 'SGGEV1',
753 \$ result( 2 ), n, jtype, ioldsd
754 END IF
755*
756* Do the tests (3) and (4)
757*
758 CALL sget52( .false., n, a, lda, b, lda, z, ldq, alphar,
759 \$ alphai, beta, work, result( 3 ) )
760 IF( result( 4 ).GT.thresh ) THEN
761 WRITE( nounit, fmt = 9998 )'Right', 'SGGEV1',
762 \$ result( 4 ), n, jtype, ioldsd
763 END IF
764*
765* Do the test (5)
766*
767 CALL slacpy( ' ', n, n, a, lda, s, lda )
768 CALL slacpy( ' ', n, n, b, lda, t, lda )
769 CALL sggev( 'N', 'N', n, s, lda, t, lda, alphr1, alphi1,
770 \$ beta1, q, ldq, z, ldq, work, lwork, ierr )
771 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
772 result( 1 ) = ulpinv
773 WRITE( nounit, fmt = 9999 )'SGGEV2', ierr, n, jtype,
774 \$ ioldsd
775 info = abs( ierr )
776 GO TO 190
777 END IF
778*
779 DO 120 j = 1, n
780 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
781 \$ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )
782 \$ result( 5 ) = ulpinv
783 120 CONTINUE
784*
785* Do the test (6): Compute eigenvalues and left eigenvectors,
786* and test them
787*
788 CALL slacpy( ' ', n, n, a, lda, s, lda )
789 CALL slacpy( ' ', n, n, b, lda, t, lda )
790 CALL sggev( 'V', 'N', n, s, lda, t, lda, alphr1, alphi1,
791 \$ beta1, qe, ldqe, z, ldq, work, lwork, ierr )
792 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
793 result( 1 ) = ulpinv
794 WRITE( nounit, fmt = 9999 )'SGGEV3', ierr, n, jtype,
795 \$ ioldsd
796 info = abs( ierr )
797 GO TO 190
798 END IF
799*
800 DO 130 j = 1, n
801 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
802 \$ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )
803 \$ result( 6 ) = ulpinv
804 130 CONTINUE
805*
806 DO 150 j = 1, n
807 DO 140 jc = 1, n
808 IF( q( j, jc ).NE.qe( j, jc ) )
809 \$ result( 6 ) = ulpinv
810 140 CONTINUE
811 150 CONTINUE
812*
813* DO the test (7): Compute eigenvalues and right eigenvectors,
814* and test them
815*
816 CALL slacpy( ' ', n, n, a, lda, s, lda )
817 CALL slacpy( ' ', n, n, b, lda, t, lda )
818 CALL sggev( 'N', 'V', n, s, lda, t, lda, alphr1, alphi1,
819 \$ beta1, q, ldq, qe, ldqe, work, lwork, ierr )
820 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
821 result( 1 ) = ulpinv
822 WRITE( nounit, fmt = 9999 )'SGGEV4', ierr, n, jtype,
823 \$ ioldsd
824 info = abs( ierr )
825 GO TO 190
826 END IF
827*
828 DO 160 j = 1, n
829 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
830 \$ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )
831 \$ result( 7 ) = ulpinv
832 160 CONTINUE
833*
834 DO 180 j = 1, n
835 DO 170 jc = 1, n
836 IF( z( j, jc ).NE.qe( j, jc ) )
837 \$ result( 7 ) = ulpinv
838 170 CONTINUE
839 180 CONTINUE
840*
841* End of Loop -- Check for RESULT(j) > THRESH
842*
843 190 CONTINUE
844*
845 ntestt = ntestt + 7
846*
847* Print out tests which fail.
848*
849 DO 200 jr = 1, 7
850 IF( result( jr ).GE.thresh ) THEN
851*
852* If this is the first test to fail,
853* print a header to the data file.
854*
855 IF( nerrs.EQ.0 ) THEN
856 WRITE( nounit, fmt = 9997 )'SGV'
857*
858* Matrix types
859*
860 WRITE( nounit, fmt = 9996 )
861 WRITE( nounit, fmt = 9995 )
862 WRITE( nounit, fmt = 9994 )'Orthogonal'
863*
864* Tests performed
865*
866 WRITE( nounit, fmt = 9993 )
867*
868 END IF
869 nerrs = nerrs + 1
870 IF( result( jr ).LT.10000.0 ) THEN
871 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
872 \$ result( jr )
873 ELSE
874 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
875 \$ result( jr )
876 END IF
877 END IF
878 200 CONTINUE
879*
880 210 CONTINUE
881 220 CONTINUE
882*
883* Summary
884*
885 CALL alasvm( 'SGV', nounit, nerrs, ntestt, 0 )
886*
887 work( 1 ) = maxwrk
888*
889 RETURN
890*
891 9999 FORMAT( ' SDRGEV: ', a, ' returned INFO=', i6, '.', / 3x, 'N=',
892 \$ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
893*
894 9998 FORMAT( ' SDRGEV: ', a, ' Eigenvectors from ', a, ' incorrectly ',
895 \$ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 3x,
896 \$ 'N=', i4, ', JTYPE=', i3, ', ISEED=(', 4( i4, ',' ), i5,
897 \$ ')' )
898*
899 9997 FORMAT( / 1x, a3, ' -- Real Generalized eigenvalue problem driver'
900 \$ )
901*
902 9996 FORMAT( ' Matrix types (see SDRGEV for details): ' )
903*
904 9995 FORMAT( ' Special Matrices:', 23x,
905 \$ '(J''=transposed Jordan block)',
906 \$ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
907 \$ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
908 \$ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
909 \$ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
910 \$ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
911 \$ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
912 9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
913 \$ / ' 16=Transposed Jordan Blocks 19=geometric ',
914 \$ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
915 \$ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
916 \$ 'alpha, beta=0,1 21=random alpha, beta=0,1',
917 \$ / ' Large & Small Matrices:', / ' 22=(large, small) ',
918 \$ '23=(small,large) 24=(small,small) 25=(large,large)',
919 \$ / ' 26=random O(1) matrices.' )
920*
921 9993 FORMAT( / ' Tests performed: ',
922 \$ / ' 1 = max | ( b A - a B )''*l | / const.,',
923 \$ / ' 2 = | |VR(i)| - 1 | / ulp,',
924 \$ / ' 3 = max | ( b A - a B )*r | / const.',
925 \$ / ' 4 = | |VL(i)| - 1 | / ulp,',
926 \$ / ' 5 = 0 if W same no matter if r or l computed,',
927 \$ / ' 6 = 0 if l same no matter if l computed,',
928 \$ / ' 7 = 0 if r same no matter if r computed,', / 1x )
929 9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
930 \$ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
931 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
932 \$ 4( i4, ',' ), ' result ', i2, ' is', 1p, e10.3 )
933*
934* End of SDRGEV
935*
936 END
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sggev(jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
SGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition sggev.f:226
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 sdrgev(nsizes, nn, ntypes, dotype, iseed, thresh, nounit, a, lda, b, s, t, q, ldq, z, qe, ldqe, alphar, alphai, beta, alphr1, alphi1, beta1, work, lwork, result, info)
SDRGEV
Definition sdrgev.f:408
subroutine sget52(left, n, a, lda, b, ldb, e, lde, alphar, alphai, beta, work, result)
SGET52
Definition sget52.f:199
subroutine slatm4(itype, n, nz1, nz2, isign, amagn, rcond, triang, idist, iseed, a, lda)
SLATM4
Definition slatm4.f:175