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