LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cdrges.f
Go to the documentation of this file.
1*> \brief \b CDRGES
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 CDRGES( 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*> CDRGES checks the nonsymmetric generalized eigenvalue (Schur form)
35*> problem driver CGGES.
36*>
37*> CGGES 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 CDRGES 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*> SDRGES 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, SDRGES
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 SDRGES 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 CGGES. 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 CGGES.
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 CGGES.
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 CGGES.
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 CGGES.
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 cdrges( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
379 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHA,
380 $ BETA, WORK, LWORK, RWORK, RESULT, BWORK, INFO )
381*
382* -- LAPACK test routine --
383* -- LAPACK is a software package provided by Univ. of Tennessee, --
384* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
385*
386* .. Scalar Arguments ..
387 INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
388 REAL THRESH
389* ..
390* .. Array Arguments ..
391 LOGICAL BWORK( * ), DOTYPE( * )
392 INTEGER ISEED( 4 ), NN( * )
393 REAL RESULT( 13 ), RWORK( * )
394 COMPLEX A( LDA, * ), ALPHA( * ), B( LDA, * ),
395 $ beta( * ), q( ldq, * ), s( lda, * ),
396 $ t( lda, * ), work( * ), z( ldq, * )
397* ..
398*
399* =====================================================================
400*
401* .. Parameters ..
402 REAL ZERO, ONE
403 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
404 COMPLEX CZERO, CONE
405 parameter( czero = ( 0.0e+0, 0.0e+0 ),
406 $ cone = ( 1.0e+0, 0.0e+0 ) )
407 INTEGER MAXTYP
408 PARAMETER ( MAXTYP = 26 )
409* ..
410* .. Local Scalars ..
411 LOGICAL BADNN, ILABAD
412 CHARACTER SORT
413 INTEGER I, IADD, IINFO, IN, ISORT, J, JC, JR, JSIZE,
414 $ jtype, knteig, maxwrk, minwrk, mtypes, n, n1,
415 $ nb, nerrs, nmats, nmax, ntest, ntestt, rsub,
416 $ sdim
417 REAL SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV
418 COMPLEX CTEMP, X
419* ..
420* .. Local Arrays ..
421 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
422 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
423 $ katype( maxtyp ), kazero( maxtyp ),
424 $ kbmagn( maxtyp ), kbtype( maxtyp ),
425 $ kbzero( maxtyp ), kclass( maxtyp ),
426 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
427 REAL RMAGN( 0: 3 )
428* ..
429* .. External Functions ..
430 LOGICAL CLCTES
431 INTEGER ILAENV
432 REAL SLAMCH
433 COMPLEX CLARND
434 EXTERNAL clctes, ilaenv, slamch, clarnd
435* ..
436* .. External Subroutines ..
437 EXTERNAL alasvm, cget51, cget54, cgges, clacpy, clarfg,
439* ..
440* .. Intrinsic Functions ..
441 INTRINSIC abs, aimag, conjg, max, min, real, sign
442* ..
443* .. Statement Functions ..
444 REAL ABS1
445* ..
446* .. Statement Function definitions ..
447 abs1( x ) = abs( real( x ) ) + abs( aimag( x ) )
448* ..
449* .. Data statements ..
450 DATA kclass / 15*1, 10*2, 1*3 /
451 DATA kz1 / 0, 1, 2, 1, 3, 3 /
452 DATA kz2 / 0, 0, 1, 2, 1, 1 /
453 DATA kadd / 0, 0, 0, 0, 3, 2 /
454 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
455 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
456 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
457 $ 1, 1, -4, 2, -4, 8*8, 0 /
458 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
459 $ 4*5, 4*3, 1 /
460 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
461 $ 4*6, 4*4, 1 /
462 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
463 $ 2, 1 /
464 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
465 $ 2, 1 /
466 DATA ktrian / 16*0, 10*1 /
467 DATA lasign / 6*.false., .true., .false., 2*.true.,
468 $ 2*.false., 3*.true., .false., .true.,
469 $ 3*.false., 5*.true., .false. /
470 DATA lbsign / 7*.false., .true., 2*.false.,
471 $ 2*.true., 2*.false., .true., .false., .true.,
472 $ 9*.false. /
473* ..
474* .. Executable Statements ..
475*
476* Check for errors
477*
478 info = 0
479*
480 badnn = .false.
481 nmax = 1
482 DO 10 j = 1, nsizes
483 nmax = max( nmax, nn( j ) )
484 IF( nn( j ).LT.0 )
485 $ badnn = .true.
486 10 CONTINUE
487*
488 IF( nsizes.LT.0 ) THEN
489 info = -1
490 ELSE IF( badnn ) THEN
491 info = -2
492 ELSE IF( ntypes.LT.0 ) THEN
493 info = -3
494 ELSE IF( thresh.LT.zero ) THEN
495 info = -6
496 ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
497 info = -9
498 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
499 info = -14
500 END IF
501*
502* Compute workspace
503* (Note: Comments in the code beginning "Workspace:" describe the
504* minimal amount of workspace needed at that point in the code,
505* as well as the preferred amount for good performance.
506* NB refers to the optimal block size for the immediately
507* following subroutine, as returned by ILAENV.
508*
509 minwrk = 1
510 IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
511 minwrk = 3*nmax*nmax
512 nb = max( 1, ilaenv( 1, 'CGEQRF', ' ', nmax, nmax, -1, -1 ),
513 $ ilaenv( 1, 'CUNMQR', 'LC', nmax, nmax, nmax, -1 ),
514 $ ilaenv( 1, 'CUNGQR', ' ', nmax, nmax, nmax, -1 ) )
515 maxwrk = max( nmax+nmax*nb, 3*nmax*nmax )
516 work( 1 ) = maxwrk
517 END IF
518*
519 IF( lwork.LT.minwrk )
520 $ info = -19
521*
522 IF( info.NE.0 ) THEN
523 CALL xerbla( 'CDRGES', -info )
524 RETURN
525 END IF
526*
527* Quick return if possible
528*
529 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
530 $ RETURN
531*
532 ulp = slamch( 'Precision' )
533 safmin = slamch( 'Safe minimum' )
534 safmin = safmin / ulp
535 safmax = one / safmin
536 ulpinv = one / ulp
537*
538* The values RMAGN(2:3) depend on N, see below.
539*
540 rmagn( 0 ) = zero
541 rmagn( 1 ) = one
542*
543* Loop over matrix sizes
544*
545 ntestt = 0
546 nerrs = 0
547 nmats = 0
548*
549 DO 190 jsize = 1, nsizes
550 n = nn( jsize )
551 n1 = max( 1, n )
552 rmagn( 2 ) = safmax*ulp / real( n1 )
553 rmagn( 3 ) = safmin*ulpinv*real( n1 )
554*
555 IF( nsizes.NE.1 ) THEN
556 mtypes = min( maxtyp, ntypes )
557 ELSE
558 mtypes = min( maxtyp+1, ntypes )
559 END IF
560*
561* Loop over matrix types
562*
563 DO 180 jtype = 1, mtypes
564 IF( .NOT.dotype( jtype ) )
565 $ GO TO 180
566 nmats = nmats + 1
567 ntest = 0
568*
569* Save ISEED in case of an error.
570*
571 DO 20 j = 1, 4
572 ioldsd( j ) = iseed( j )
573 20 CONTINUE
574*
575* Initialize RESULT
576*
577 DO 30 j = 1, 13
578 result( j ) = zero
579 30 CONTINUE
580*
581* Generate test matrices A and B
582*
583* Description of control parameters:
584*
585* KCLASS: =1 means w/o rotation, =2 means w/ rotation,
586* =3 means random.
587* KATYPE: the "type" to be passed to CLATM4 for computing A.
588* KAZERO: the pattern of zeros on the diagonal for A:
589* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
590* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
591* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
592* non-zero entries.)
593* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
594* =2: large, =3: small.
595* LASIGN: .TRUE. if the diagonal elements of A are to be
596* multiplied by a random magnitude 1 number.
597* KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
598* KTRIAN: =0: don't fill in the upper triangle, =1: do.
599* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
600* RMAGN: used to implement KAMAGN and KBMAGN.
601*
602 IF( mtypes.GT.maxtyp )
603 $ GO TO 110
604 iinfo = 0
605 IF( kclass( jtype ).LT.3 ) THEN
606*
607* Generate A (w/o rotation)
608*
609 IF( abs( katype( jtype ) ).EQ.3 ) THEN
610 in = 2*( ( n-1 ) / 2 ) + 1
611 IF( in.NE.n )
612 $ CALL claset( 'Full', n, n, czero, czero, a, lda )
613 ELSE
614 in = n
615 END IF
616 CALL clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
617 $ kz2( kazero( jtype ) ), lasign( jtype ),
618 $ rmagn( kamagn( jtype ) ), ulp,
619 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
620 $ iseed, a, lda )
621 iadd = kadd( kazero( jtype ) )
622 IF( iadd.GT.0 .AND. iadd.LE.n )
623 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
624*
625* Generate B (w/o rotation)
626*
627 IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
628 in = 2*( ( n-1 ) / 2 ) + 1
629 IF( in.NE.n )
630 $ CALL claset( 'Full', n, n, czero, czero, b, lda )
631 ELSE
632 in = n
633 END IF
634 CALL clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
635 $ kz2( kbzero( jtype ) ), lbsign( jtype ),
636 $ rmagn( kbmagn( jtype ) ), one,
637 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
638 $ iseed, b, lda )
639 iadd = kadd( kbzero( jtype ) )
640 IF( iadd.NE.0 .AND. iadd.LE.n )
641 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
642*
643 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
644*
645* Include rotations
646*
647* Generate Q, Z as Householder transformations times
648* a diagonal matrix.
649*
650 DO 50 jc = 1, n - 1
651 DO 40 jr = jc, n
652 q( jr, jc ) = clarnd( 3, iseed )
653 z( jr, jc ) = clarnd( 3, iseed )
654 40 CONTINUE
655 CALL clarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
656 $ work( jc ) )
657 work( 2*n+jc ) = sign( one, real( q( jc, jc ) ) )
658 q( jc, jc ) = cone
659 CALL clarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
660 $ work( n+jc ) )
661 work( 3*n+jc ) = sign( one, real( z( jc, jc ) ) )
662 z( jc, jc ) = cone
663 50 CONTINUE
664 ctemp = clarnd( 3, iseed )
665 q( n, n ) = cone
666 work( n ) = czero
667 work( 3*n ) = ctemp / abs( ctemp )
668 ctemp = clarnd( 3, iseed )
669 z( n, n ) = cone
670 work( 2*n ) = czero
671 work( 4*n ) = ctemp / abs( ctemp )
672*
673* Apply the diagonal matrices
674*
675 DO 70 jc = 1, n
676 DO 60 jr = 1, n
677 a( jr, jc ) = work( 2*n+jr )*
678 $ conjg( work( 3*n+jc ) )*
679 $ a( jr, jc )
680 b( jr, jc ) = work( 2*n+jr )*
681 $ conjg( work( 3*n+jc ) )*
682 $ b( jr, jc )
683 60 CONTINUE
684 70 CONTINUE
685 CALL cunm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
686 $ lda, work( 2*n+1 ), iinfo )
687 IF( iinfo.NE.0 )
688 $ GO TO 100
689 CALL cunm2r( 'R', 'C', n, n, n-1, z, ldq, work( n+1 ),
690 $ a, lda, work( 2*n+1 ), iinfo )
691 IF( iinfo.NE.0 )
692 $ GO TO 100
693 CALL cunm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
694 $ lda, work( 2*n+1 ), iinfo )
695 IF( iinfo.NE.0 )
696 $ GO TO 100
697 CALL cunm2r( 'R', 'C', n, n, n-1, z, ldq, work( n+1 ),
698 $ b, lda, work( 2*n+1 ), iinfo )
699 IF( iinfo.NE.0 )
700 $ GO TO 100
701 END IF
702 ELSE
703*
704* Random matrices
705*
706 DO 90 jc = 1, n
707 DO 80 jr = 1, n
708 a( jr, jc ) = rmagn( kamagn( jtype ) )*
709 $ clarnd( 4, iseed )
710 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
711 $ clarnd( 4, iseed )
712 80 CONTINUE
713 90 CONTINUE
714 END IF
715*
716 100 CONTINUE
717*
718 IF( iinfo.NE.0 ) THEN
719 WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
720 $ ioldsd
721 info = abs( iinfo )
722 RETURN
723 END IF
724*
725 110 CONTINUE
726*
727 DO 120 i = 1, 13
728 result( i ) = -one
729 120 CONTINUE
730*
731* Test with and without sorting of eigenvalues
732*
733 DO 150 isort = 0, 1
734 IF( isort.EQ.0 ) THEN
735 sort = 'N'
736 rsub = 0
737 ELSE
738 sort = 'S'
739 rsub = 5
740 END IF
741*
742* Call CGGES to compute H, T, Q, Z, alpha, and beta.
743*
744 CALL clacpy( 'Full', n, n, a, lda, s, lda )
745 CALL clacpy( 'Full', n, n, b, lda, t, lda )
746 ntest = 1 + rsub + isort
747 result( 1+rsub+isort ) = ulpinv
748 CALL cgges( 'V', 'V', sort, clctes, n, s, lda, t, lda,
749 $ sdim, alpha, beta, q, ldq, z, ldq, work,
750 $ lwork, rwork, bwork, iinfo )
751 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
752 result( 1+rsub+isort ) = ulpinv
753 WRITE( nounit, fmt = 9999 )'CGGES', iinfo, n, jtype,
754 $ ioldsd
755 info = abs( iinfo )
756 GO TO 160
757 END IF
758*
759 ntest = 4 + rsub
760*
761* Do tests 1--4 (or tests 7--9 when reordering )
762*
763 IF( isort.EQ.0 ) THEN
764 CALL cget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
765 $ work, rwork, result( 1 ) )
766 CALL cget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
767 $ work, rwork, result( 2 ) )
768 ELSE
769 CALL cget54( n, a, lda, b, lda, s, lda, t, lda, q,
770 $ ldq, z, ldq, work, result( 2+rsub ) )
771 END IF
772*
773 CALL cget51( 3, n, b, lda, t, lda, q, ldq, q, ldq, work,
774 $ rwork, result( 3+rsub ) )
775 CALL cget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
776 $ rwork, result( 4+rsub ) )
777*
778* Do test 5 and 6 (or Tests 10 and 11 when reordering):
779* check Schur form of A and compare eigenvalues with
780* diagonals.
781*
782 ntest = 6 + rsub
783 temp1 = zero
784*
785 DO 130 j = 1, n
786 ilabad = .false.
787 temp2 = ( abs1( alpha( j )-s( j, j ) ) /
788 $ max( safmin, abs1( alpha( j ) ), abs1( s( j,
789 $ j ) ) )+abs1( beta( j )-t( j, j ) ) /
790 $ max( safmin, abs1( beta( j ) ), abs1( t( j,
791 $ j ) ) ) ) / ulp
792*
793 IF( j.LT.n ) THEN
794 IF( s( j+1, j ).NE.zero ) THEN
795 ilabad = .true.
796 result( 5+rsub ) = ulpinv
797 END IF
798 END IF
799 IF( j.GT.1 ) THEN
800 IF( s( j, j-1 ).NE.zero ) THEN
801 ilabad = .true.
802 result( 5+rsub ) = ulpinv
803 END IF
804 END IF
805 temp1 = max( temp1, temp2 )
806 IF( ilabad ) THEN
807 WRITE( nounit, fmt = 9998 )j, n, jtype, ioldsd
808 END IF
809 130 CONTINUE
810 result( 6+rsub ) = temp1
811*
812 IF( isort.GE.1 ) THEN
813*
814* Do test 12
815*
816 ntest = 12
817 result( 12 ) = zero
818 knteig = 0
819 DO 140 i = 1, n
820 IF( clctes( alpha( i ), beta( i ) ) )
821 $ knteig = knteig + 1
822 140 CONTINUE
823 IF( sdim.NE.knteig )
824 $ result( 13 ) = ulpinv
825 END IF
826*
827 150 CONTINUE
828*
829* End of Loop -- Check for RESULT(j) > THRESH
830*
831 160 CONTINUE
832*
833 ntestt = ntestt + ntest
834*
835* Print out tests which fail.
836*
837 DO 170 jr = 1, ntest
838 IF( result( jr ).GE.thresh ) THEN
839*
840* If this is the first test to fail,
841* print a header to the data file.
842*
843 IF( nerrs.EQ.0 ) THEN
844 WRITE( nounit, fmt = 9997 )'CGS'
845*
846* Matrix types
847*
848 WRITE( nounit, fmt = 9996 )
849 WRITE( nounit, fmt = 9995 )
850 WRITE( nounit, fmt = 9994 )'Unitary'
851*
852* Tests performed
853*
854 WRITE( nounit, fmt = 9993 )'unitary', '''',
855 $ 'transpose', ( '''', j = 1, 8 )
856*
857 END IF
858 nerrs = nerrs + 1
859 IF( result( jr ).LT.10000.0 ) THEN
860 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
861 $ result( jr )
862 ELSE
863 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
864 $ result( jr )
865 END IF
866 END IF
867 170 CONTINUE
868*
869 180 CONTINUE
870 190 CONTINUE
871*
872* Summary
873*
874 CALL alasvm( 'CGS', nounit, nerrs, ntestt, 0 )
875*
876 work( 1 ) = maxwrk
877*
878 RETURN
879*
880 9999 FORMAT( ' CDRGES: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
881 $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
882*
883 9998 FORMAT( ' CDRGES: S not in Schur form at eigenvalue ', i6, '.',
884 $ / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
885 $ i5, ')' )
886*
887 9997 FORMAT( / 1x, a3, ' -- Complex Generalized Schur from problem ',
888 $ 'driver' )
889*
890 9996 FORMAT( ' Matrix types (see CDRGES for details): ' )
891*
892 9995 FORMAT( ' Special Matrices:', 23x,
893 $ '(J''=transposed Jordan block)',
894 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
895 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
896 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
897 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
898 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
899 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
900 9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
901 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
902 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
903 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
904 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
905 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
906 $ '23=(small,large) 24=(small,small) 25=(large,large)',
907 $ / ' 26=random O(1) matrices.' )
908*
909 9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
910 $ 'Q and Z are ', a, ',', / 19x,
911 $ 'l and r are the appropriate left and right', / 19x,
912 $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
913 $ ' means ', a, '.)', / ' Without ordering: ',
914 $ / ' 1 = | A - Q S Z', a,
915 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
916 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
917 $ ' | / ( n ulp ) 4 = | I - ZZ', a,
918 $ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
919 $ / ' 6 = difference between (alpha,beta)',
920 $ ' and diagonals of (S,T)', / ' With ordering: ',
921 $ / ' 7 = | (A,B) - Q (S,T) Z', a, ' | / ( |(A,B)| n ulp )',
922 $ / ' 8 = | I - QQ', a,
923 $ ' | / ( n ulp ) 9 = | I - ZZ', a,
924 $ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
925 $ / ' 11 = difference between (alpha,beta) and diagonals',
926 $ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
927 $ 'selected eigenvalues', / )
928 9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
929 $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
930 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
931 $ 4( i4, ',' ), ' result ', i2, ' is', 1p, e10.3 )
932*
933* End of CDRGES
934*
935 END
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cdrges(nsizes, nn, ntypes, dotype, iseed, thresh, nounit, a, lda, b, s, t, q, ldq, z, alpha, beta, work, lwork, rwork, result, bwork, info)
CDRGES
Definition cdrges.f:381
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 clatm4(itype, n, nz1, nz2, rsign, amagn, rcond, triang, idist, iseed, a, lda)
CLATM4
Definition clatm4.f:171
subroutine cgges(jobvsl, jobvsr, sort, selctg, n, a, lda, b, ldb, sdim, alpha, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, rwork, bwork, info)
CGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE m...
Definition cgges.f:270
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 clarfg(n, alpha, x, incx, tau)
CLARFG generates an elementary reflector (Householder matrix).
Definition clarfg.f:106
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 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