OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cdrgev3.f
Go to the documentation of this file.
1*> \brief \b CDRGEV3
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 CDRGEV3( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12* NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE,
13* ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK, RWORK,
14* 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 RESULT( * ), RWORK( * )
25* COMPLEX A( LDA, * ), ALPHA( * ), ALPHA1( * ),
26* $ B( LDA, * ), BETA( * ), BETA1( * ),
27* $ Q( LDQ, * ), QE( LDQE, * ), S( LDA, * ),
28* $ T( LDA, * ), WORK( * ), Z( LDQ, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> CDRGEV3 checks the nonsymmetric generalized eigenvalue problem driver
38*> routine CGGEV3.
39*>
40*> CGGEV3 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 CDRGEV3 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 CGGEV3:
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*> CDRGEV3 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, CDRGEV3
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 CDRGEV3 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 COMPLEX array, dimension(LDA, max(NN))
257*> Used to hold the original A matrix. Used as input only
258*> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
259*> DOTYPE(MAXTYP+1)=.TRUE.
260*> \endverbatim
261*>
262*> \param[in] LDA
263*> \verbatim
264*> LDA is INTEGER
265*> The leading dimension of A, B, S, and T.
266*> It must be at least 1 and at least max( NN ).
267*> \endverbatim
268*>
269*> \param[in,out] B
270*> \verbatim
271*> B is COMPLEX array, dimension(LDA, max(NN))
272*> Used to hold the original B matrix. Used as input only
273*> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
274*> DOTYPE(MAXTYP+1)=.TRUE.
275*> \endverbatim
276*>
277*> \param[out] S
278*> \verbatim
279*> S is COMPLEX array, dimension (LDA, max(NN))
280*> The Schur form matrix computed from A by CGGEV3. On exit, S
281*> contains the Schur form matrix corresponding to the matrix
282*> in A.
283*> \endverbatim
284*>
285*> \param[out] T
286*> \verbatim
287*> T is COMPLEX array, dimension (LDA, max(NN))
288*> The upper triangular matrix computed from B by CGGEV3.
289*> \endverbatim
290*>
291*> \param[out] Q
292*> \verbatim
293*> Q is COMPLEX array, dimension (LDQ, max(NN))
294*> The (left) eigenvectors matrix computed by CGGEV3.
295*> \endverbatim
296*>
297*> \param[in] LDQ
298*> \verbatim
299*> LDQ is INTEGER
300*> The leading dimension of Q and Z. It must
301*> be at least 1 and at least max( NN ).
302*> \endverbatim
303*>
304*> \param[out] Z
305*> \verbatim
306*> Z is COMPLEX array, dimension( LDQ, max(NN) )
307*> The (right) orthogonal matrix computed by CGGEV3.
308*> \endverbatim
309*>
310*> \param[out] QE
311*> \verbatim
312*> QE is COMPLEX array, dimension( LDQ, max(NN) )
313*> QE holds the computed right or left eigenvectors.
314*> \endverbatim
315*>
316*> \param[in] LDQE
317*> \verbatim
318*> LDQE is INTEGER
319*> The leading dimension of QE. LDQE >= max(1,max(NN)).
320*> \endverbatim
321*>
322*> \param[out] ALPHA
323*> \verbatim
324*> ALPHA is COMPLEX array, dimension (max(NN))
325*> \endverbatim
326*>
327*> \param[out] BETA
328*> \verbatim
329*> BETA is COMPLEX array, dimension (max(NN))
330*>
331*> The generalized eigenvalues of (A,B) computed by CGGEV3.
332*> ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
333*> generalized eigenvalue of A and B.
334*> \endverbatim
335*>
336*> \param[out] ALPHA1
337*> \verbatim
338*> ALPHA1 is COMPLEX array, dimension (max(NN))
339*> \endverbatim
340*>
341*> \param[out] BETA1
342*> \verbatim
343*> BETA1 is COMPLEX array, dimension (max(NN))
344*>
345*> Like ALPHAR, ALPHAI, BETA, these arrays contain the
346*> eigenvalues of A and B, but those computed when CGGEV3 only
347*> computes a partial eigendecomposition, i.e. not the
348*> eigenvalues and left and right eigenvectors.
349*> \endverbatim
350*>
351*> \param[out] WORK
352*> \verbatim
353*> WORK is COMPLEX array, dimension (LWORK)
354*> \endverbatim
355*>
356*> \param[in] LWORK
357*> \verbatim
358*> LWORK is INTEGER
359*> The number of entries in WORK. LWORK >= N*(N+1)
360*> \endverbatim
361*>
362*> \param[out] RWORK
363*> \verbatim
364*> RWORK is REAL array, dimension (8*N)
365*> Real workspace.
366*> \endverbatim
367*>
368*> \param[out] RESULT
369*> \verbatim
370*> RESULT is REAL array, dimension (2)
371*> The values computed by the tests described above.
372*> The values are currently limited to 1/ulp, to avoid overflow.
373*> \endverbatim
374*>
375*> \param[out] INFO
376*> \verbatim
377*> INFO is INTEGER
378*> = 0: successful exit
379*> < 0: if INFO = -i, the i-th argument had an illegal value.
380*> > 0: A routine returned an error code. INFO is the
381*> absolute value of the INFO value returned.
382*> \endverbatim
383*
384* Authors:
385* ========
386*
387*> \author Univ. of Tennessee
388*> \author Univ. of California Berkeley
389*> \author Univ. of Colorado Denver
390*> \author NAG Ltd.
391*
392*> \ingroup complex_eig
393*
394* =====================================================================
395 SUBROUTINE cdrgev3( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
396 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE,
397 $ ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK,
398 $ RWORK, RESULT, INFO )
399*
400* -- LAPACK test routine --
401* -- LAPACK is a software package provided by Univ. of Tennessee, --
402* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
403*
404* .. Scalar Arguments ..
405 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
406 $ NTYPES
407 REAL THRESH
408* ..
409* .. Array Arguments ..
410 LOGICAL DOTYPE( * )
411 INTEGER ISEED( 4 ), NN( * )
412 REAL RESULT( * ), RWORK( * )
413 COMPLEX A( LDA, * ), ALPHA( * ), ALPHA1( * ),
414 $ b( lda, * ), beta( * ), beta1( * ),
415 $ q( ldq, * ), qe( ldqe, * ), s( lda, * ),
416 $ t( lda, * ), work( * ), z( ldq, * )
417* ..
418*
419* =====================================================================
420*
421* .. Parameters ..
422 REAL ZERO, ONE
423 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
424 COMPLEX CZERO, CONE
425 parameter( czero = ( 0.0e+0, 0.0e+0 ),
426 $ cone = ( 1.0e+0, 0.0e+0 ) )
427 INTEGER MAXTYP
428 parameter( maxtyp = 26 )
429* ..
430* .. Local Scalars ..
431 LOGICAL BADNN
432 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
433 $ MAXWRK, MINWRK, MTYPES, N, N1, NB, NERRS,
434 $ nmats, nmax, ntestt
435 REAL SAFMAX, SAFMIN, ULP, ULPINV
436 COMPLEX CTEMP
437* ..
438* .. Local Arrays ..
439 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
440 INTEGER 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 INTEGER ILAENV
449 REAL SLAMCH
450 COMPLEX CLARND
451 EXTERNAL ilaenv, slamch, clarnd
452* ..
453* .. External Subroutines ..
454 EXTERNAL alasvm, cget52, cggev3, clacpy, clarfg, claset,
456* ..
457* .. Intrinsic Functions ..
458 INTRINSIC abs, conjg, 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 lasign / 6*.false., .true., .false., 2*.true.,
479 $ 2*.false., 3*.true., .false., .true.,
480 $ 3*.false., 5*.true., .false. /
481 DATA lbsign / 7*.false., .true., 2*.false.,
482 $ 2*.true., 2*.false., .true., .false., .true.,
483 $ 9*.false. /
484* ..
485* .. Executable Statements ..
486*
487* Check for errors
488*
489 info = 0
490*
491 badnn = .false.
492 nmax = 1
493 DO 10 j = 1, nsizes
494 nmax = max( nmax, nn( j ) )
495 IF( nn( j ).LT.0 )
496 $ badnn = .true.
497 10 CONTINUE
498*
499 IF( nsizes.LT.0 ) THEN
500 info = -1
501 ELSE IF( badnn ) THEN
502 info = -2
503 ELSE IF( ntypes.LT.0 ) THEN
504 info = -3
505 ELSE IF( thresh.LT.zero ) THEN
506 info = -6
507 ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
508 info = -9
509 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
510 info = -14
511 ELSE IF( ldqe.LE.1 .OR. ldqe.LT.nmax ) THEN
512 info = -17
513 END IF
514*
515* Compute workspace
516* (Note: Comments in the code beginning "Workspace:" describe the
517* minimal amount of workspace needed at that point in the code,
518* as well as the preferred amount for good performance.
519* NB refers to the optimal block size for the immediately
520* following subroutine, as returned by ILAENV.
521*
522 minwrk = 1
523 IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
524 minwrk = nmax*( nmax+1 )
525 nb = max( 1, ilaenv( 1, 'CGEQRF', ' ', nmax, nmax, -1, -1 ),
526 $ ilaenv( 1, 'CUNMQR', 'LC', nmax, nmax, nmax, -1 ),
527 $ ilaenv( 1, 'CUNGQR', ' ', nmax, nmax, nmax, -1 ) )
528 maxwrk = max( 2*nmax, nmax*( nb+1 ), nmax*( nmax+1 ) )
529 work( 1 ) = maxwrk
530 END IF
531*
532 IF( lwork.LT.minwrk )
533 $ info = -23
534*
535 IF( info.NE.0 ) THEN
536 CALL xerbla( 'CDRGEV3', -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 ulp = slamch( 'Precision' )
546 safmin = slamch( 'Safe minimum' )
547 safmin = safmin / ulp
548 safmax = one / safmin
549 CALL slabad( safmin, safmax )
550 ulpinv = one / ulp
551*
552* The values RMAGN(2:3) depend on N, see below.
553*
554 rmagn( 0 ) = zero
555 rmagn( 1 ) = one
556*
557* Loop over sizes, types
558*
559 ntestt = 0
560 nerrs = 0
561 nmats = 0
562*
563 DO 220 jsize = 1, nsizes
564 n = nn( jsize )
565 n1 = max( 1, n )
566 rmagn( 2 ) = safmax*ulp / real( n1 )
567 rmagn( 3 ) = safmin*ulpinv*n1
568*
569 IF( nsizes.NE.1 ) THEN
570 mtypes = min( maxtyp, ntypes )
571 ELSE
572 mtypes = min( maxtyp+1, ntypes )
573 END IF
574*
575 DO 210 jtype = 1, mtypes
576 IF( .NOT.dotype( jtype ) )
577 $ GO TO 210
578 nmats = nmats + 1
579*
580* Save ISEED in case of an error.
581*
582 DO 20 j = 1, 4
583 ioldsd( j ) = iseed( j )
584 20 CONTINUE
585*
586* Generate test matrices A and B
587*
588* Description of control parameters:
589*
590* KCLASS: =1 means w/o rotation, =2 means w/ rotation,
591* =3 means random.
592* KATYPE: the "type" to be passed to CLATM4 for computing A.
593* KAZERO: the pattern of zeros on the diagonal for A:
594* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
595* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
596* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
597* non-zero entries.)
598* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
599* =2: large, =3: small.
600* LASIGN: .TRUE. if the diagonal elements of A are to be
601* multiplied by a random magnitude 1 number.
602* KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
603* KTRIAN: =0: don't fill in the upper triangle, =1: do.
604* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
605* RMAGN: used to implement KAMAGN and KBMAGN.
606*
607 IF( mtypes.GT.maxtyp )
608 $ GO TO 100
609 ierr = 0
610 IF( kclass( jtype ).LT.3 ) THEN
611*
612* Generate A (w/o rotation)
613*
614 IF( abs( katype( jtype ) ).EQ.3 ) THEN
615 in = 2*( ( n-1 ) / 2 ) + 1
616 IF( in.NE.n )
617 $ CALL claset( 'Full', n, n, czero, czero, a, lda )
618 ELSE
619 in = n
620 END IF
621 CALL clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
622 $ kz2( kazero( jtype ) ), lasign( jtype ),
623 $ rmagn( kamagn( jtype ) ), ulp,
624 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
625 $ iseed, a, lda )
626 iadd = kadd( kazero( jtype ) )
627 IF( iadd.GT.0 .AND. iadd.LE.n )
628 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
629*
630* Generate B (w/o rotation)
631*
632 IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
633 in = 2*( ( n-1 ) / 2 ) + 1
634 IF( in.NE.n )
635 $ CALL claset( 'Full', n, n, czero, czero, b, lda )
636 ELSE
637 in = n
638 END IF
639 CALL clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
640 $ kz2( kbzero( jtype ) ), lbsign( jtype ),
641 $ rmagn( kbmagn( jtype ) ), one,
642 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
643 $ iseed, b, lda )
644 iadd = kadd( kbzero( jtype ) )
645 IF( iadd.NE.0 .AND. iadd.LE.n )
646 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
647*
648 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
649*
650* Include rotations
651*
652* Generate Q, Z as Householder transformations times
653* a diagonal matrix.
654*
655 DO 40 jc = 1, n - 1
656 DO 30 jr = jc, n
657 q( jr, jc ) = clarnd( 3, iseed )
658 z( jr, jc ) = clarnd( 3, iseed )
659 30 CONTINUE
660 CALL clarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
661 $ work( jc ) )
662 work( 2*n+jc ) = sign( one, real( q( jc, jc ) ) )
663 q( jc, jc ) = cone
664 CALL clarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
665 $ work( n+jc ) )
666 work( 3*n+jc ) = sign( one, real( z( jc, jc ) ) )
667 z( jc, jc ) = cone
668 40 CONTINUE
669 ctemp = clarnd( 3, iseed )
670 q( n, n ) = cone
671 work( n ) = czero
672 work( 3*n ) = ctemp / abs( ctemp )
673 ctemp = clarnd( 3, iseed )
674 z( n, n ) = cone
675 work( 2*n ) = czero
676 work( 4*n ) = ctemp / abs( ctemp )
677*
678* Apply the diagonal matrices
679*
680 DO 60 jc = 1, n
681 DO 50 jr = 1, n
682 a( jr, jc ) = work( 2*n+jr )*
683 $ conjg( work( 3*n+jc ) )*
684 $ a( jr, jc )
685 b( jr, jc ) = work( 2*n+jr )*
686 $ conjg( work( 3*n+jc ) )*
687 $ b( jr, jc )
688 50 CONTINUE
689 60 CONTINUE
690 CALL cunm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
691 $ lda, work( 2*n+1 ), ierr )
692 IF( ierr.NE.0 )
693 $ GO TO 90
694 CALL cunm2r( 'R', 'C', n, n, n-1, z, ldq, work( n+1 ),
695 $ a, lda, work( 2*n+1 ), ierr )
696 IF( ierr.NE.0 )
697 $ GO TO 90
698 CALL cunm2r( 'l', 'n', N, N, N-1, Q, LDQ, WORK, B,
699 $ LDA, WORK( 2*N+1 ), IERR )
700.NE. IF( IERR0 )
701 $ GO TO 90
702 CALL CUNM2R( 'r', 'c', N, N, N-1, Z, LDQ, WORK( N+1 ),
703 $ B, LDA, WORK( 2*N+1 ), IERR )
704.NE. IF( IERR0 )
705 $ GO TO 90
706 END IF
707 ELSE
708*
709* Random matrices
710*
711 DO 80 JC = 1, N
712 DO 70 JR = 1, N
713 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
714 $ CLARND( 4, ISEED )
715 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
716 $ CLARND( 4, ISEED )
717 70 CONTINUE
718 80 CONTINUE
719 END IF
720*
721 90 CONTINUE
722*
723.NE. IF( IERR0 ) THEN
724 WRITE( NOUNIT, FMT = 9999 )'generator', IERR, N, JTYPE,
725 $ IOLDSD
726 INFO = ABS( IERR )
727 RETURN
728 END IF
729*
730 100 CONTINUE
731*
732 DO 110 I = 1, 7
733 RESULT( I ) = -ONE
734 110 CONTINUE
735*
736* Call XLAENV to set the parameters used in CLAQZ0
737*
738 CALL XLAENV( 12, 10 )
739 CALL XLAENV( 13, 12 )
740 CALL XLAENV( 14, 13 )
741 CALL XLAENV( 15, 2 )
742 CALL XLAENV( 17, 10 )
743*
744* Call CGGEV3 to compute eigenvalues and eigenvectors.
745*
746 CALL CLACPY( ' ', N, N, A, LDA, S, LDA )
747 CALL CLACPY( ' ', N, N, B, LDA, T, LDA )
748 CALL CGGEV3( 'v', 'v', N, S, LDA, T, LDA, ALPHA, BETA, Q,
749 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR )
750.NE..AND..NE. IF( IERR0 IERRN+1 ) THEN
751 RESULT( 1 ) = ULPINV
752 WRITE( NOUNIT, FMT = 9999 )'cggev31', IERR, N, JTYPE,
753 $ IOLDSD
754 INFO = ABS( IERR )
755 GO TO 190
756 END IF
757*
758* Do the tests (1) and (2)
759*
760 CALL CGET52( .TRUE., N, A, LDA, B, LDA, Q, LDQ, ALPHA, BETA,
761 $ WORK, RWORK, RESULT( 1 ) )
762.GT. IF( RESULT( 2 )THRESH ) THEN
763 WRITE( NOUNIT, FMT = 9998 )'left', 'cggev31',
764 $ RESULT( 2 ), N, JTYPE, IOLDSD
765 END IF
766*
767* Do the tests (3) and (4)
768*
769 CALL CGET52( .FALSE., N, A, LDA, B, LDA, Z, LDQ, ALPHA,
770 $ BETA, WORK, RWORK, RESULT( 3 ) )
771.GT. IF( RESULT( 4 )THRESH ) THEN
772 WRITE( NOUNIT, FMT = 9998 )'right', 'cggev31',
773 $ RESULT( 4 ), N, JTYPE, IOLDSD
774 END IF
775*
776* Do test (5)
777*
778 CALL CLACPY( ' ', N, N, A, LDA, S, LDA )
779 CALL CLACPY( ' ', N, N, B, LDA, T, LDA )
780 CALL CGGEV3( 'n', 'n', N, S, LDA, T, LDA, ALPHA1, BETA1, Q,
781 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR )
782.NE..AND..NE. IF( IERR0 IERRN+1 ) THEN
783 RESULT( 1 ) = ULPINV
784 WRITE( NOUNIT, FMT = 9999 )'cggev32', IERR, N, JTYPE,
785 $ IOLDSD
786 INFO = ABS( IERR )
787 GO TO 190
788 END IF
789*
790 DO 120 J = 1, N
791.NE..OR..NE. IF( ALPHA( J )ALPHA1( J ) BETA( J )
792 $ BETA1( J ) ) RESULT( 5 ) = ULPINV
793 120 CONTINUE
794*
795* Do the test (6): Compute eigenvalues and left eigenvectors,
796* and test them
797*
798 CALL CLACPY( ' ', N, N, A, LDA, S, LDA )
799 CALL CLACPY( ' ', N, N, B, LDA, T, LDA )
800 CALL CGGEV3( 'v', 'n', N, S, LDA, T, LDA, ALPHA1, BETA1, QE,
801 $ LDQE, Z, LDQ, WORK, LWORK, RWORK, IERR )
802.NE..AND..NE. IF( IERR0 IERRN+1 ) THEN
803 RESULT( 1 ) = ULPINV
804 WRITE( NOUNIT, FMT = 9999 )'cggev33', IERR, N, JTYPE,
805 $ IOLDSD
806 INFO = ABS( IERR )
807 GO TO 190
808 END IF
809
810*
811 DO 130 J = 1, N
812.NE..OR. IF( ALPHA( J )ALPHA1( J )
813.NE. $ BETA( J )BETA1( J ) ) THEN
814 RESULT( 6 ) = ULPINV
815 ENDIF
816 130 CONTINUE
817*
818 DO 150 J = 1, N
819 DO 140 JC = 1, N
820.NE. IF( Q( J, JC )QE( J, JC ) ) THEN
821 RESULT( 6 ) = ULPINV
822 END IF
823 140 CONTINUE
824 150 CONTINUE
825*
826* DO the test (7): Compute eigenvalues and right eigenvectors,
827* and test them
828*
829 CALL CLACPY( ' ', N, N, A, LDA, S, LDA )
830 CALL CLACPY( ' ', N, N, B, LDA, T, LDA )
831 CALL CGGEV3( 'n', 'v', N, S, LDA, T, LDA, ALPHA1, BETA1, Q,
832 $ LDQ, QE, LDQE, WORK, LWORK, RWORK, IERR )
833.NE..AND..NE. IF( IERR0 IERRN+1 ) THEN
834 RESULT( 1 ) = ULPINV
835 WRITE( NOUNIT, FMT = 9999 )'cggev34', IERR, N, JTYPE,
836 $ IOLDSD
837 INFO = ABS( IERR )
838 GO TO 190
839 END IF
840*
841 DO 160 J = 1, N
842.NE..OR..NE. IF( ALPHA( J )ALPHA1( J ) BETA( J )
843 $ BETA1( J ) )RESULT( 7 ) = ULPINV
844 160 CONTINUE
845*
846 DO 180 J = 1, N
847 DO 170 JC = 1, N
848.NE. IF( Z( J, JC )QE( J, JC ) )
849 $ RESULT( 7 ) = ULPINV
850 170 CONTINUE
851 180 CONTINUE
852*
853* End of Loop -- Check for RESULT(j) > THRESH
854*
855 190 CONTINUE
856*
857 NTESTT = NTESTT + 7
858*
859* Print out tests which fail.
860*
861 DO 200 JR = 1, 7
862.GE. IF( RESULT( JR )THRESH ) THEN
863*
864* If this is the first test to fail,
865* print a header to the data file.
866*
867.EQ. IF( NERRS0 ) THEN
868 WRITE( NOUNIT, FMT = 9997 )'cgv'
869*
870* Matrix types
871*
872 WRITE( NOUNIT, FMT = 9996 )
873 WRITE( NOUNIT, FMT = 9995 )
874 WRITE( NOUNIT, FMT = 9994 )'orthogonal'
875*
876* Tests performed
877*
878 WRITE( NOUNIT, FMT = 9993 )
879*
880 END IF
881 NERRS = NERRS + 1
882.LT. IF( RESULT( JR )10000.0 ) THEN
883 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
884 $ RESULT( JR )
885 ELSE
886 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
887 $ RESULT( JR )
888 END IF
889 END IF
890 200 CONTINUE
891*
892 210 CONTINUE
893 220 CONTINUE
894*
895* Summary
896*
897 CALL ALASVM( 'cgv3', NOUNIT, NERRS, NTESTT, 0 )
898*
899 WORK( 1 ) = MAXWRK
900*
901 RETURN
902*
903 9999 FORMAT( ' cdrgev3: ', A, ' returned info=', I6, '.', / 3X, 'n=',
904 $ I6, ', jtype=', I6, ', iseed=(', 3( I5, ',' ), I5, ')' )
905*
906 9998 FORMAT( ' cdrgev3: ', A, ' eigenvectors from ', A,
907 $ ' incorrectly normalized.', / ' bits of error=', 0P, G10.3,
908 $ ',', 3X, 'n=', I4, ', jtype=', I3, ', iseed=(',
909 $ 3( I4, ',' ), I5, ')' )
910*
911 9997 FORMAT( / 1X, A3, ' -- Complex Generalized eigenvalue problem ',
912 $ 'driver' )
913*
914 9996 FORMAT( ' Matrix types (see CDRGEV3 for details): ' )
915*
916 9995 FORMAT( ' Special Matrices:', 23X,
917 $ '(J''=transposed jordan block)',
918 $ / ' 1=(0,0) 2=(i,0) 3=(0,i) 4=(i,i) 5=(j'',j'') ',
919 $ '6=(diag(j'',i), diag(i,j''))', / ' diagonal matrices: ( ',
920 $ 'd=diag(0,1,2,...) )', / ' 7=(d,i) 9=(large*d, small*i',
921 $ ') 11=(large*i, small*d) 13=(large*d, large*i)', /
922 $ ' 8=(i,d) 10=(small*d, large*i) 12=(small*i, large*d) ',
923 $ ' 14=(small*d, small*i)', / ' 15=(d, reversed d)' )
924 9994 FORMAT( ' matrices rotated by random ', A, ' matrices u, v:',
925 $ / ' 16=transposed jordan blocks 19=geometric ',
926 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
927 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
928 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
929 $ / ' large & small matrices:', / ' 22=(large, small) ',
930 $ '23=(small,large) 24=(small,small) 25=(large,large)',
931 $ / ' 26=random o(1) matrices.' )
932*
933 9993 FORMAT( / ' tests performed: ',
934 $ / ' 1 = max | ( b a - a b )''*l | / const.,',
935 $ / ' 2 = | |vr(i)| - 1 | / ulp,',
936 $ / ' 3 = max | ( b a - a b )*r | / const.',
937 $ / ' 4 = | |vl(i)| - 1 | / ulp,',
938 $ / ' 5 = 0 if w same no matter if r or l computed,',
939 $ / ' 6 = 0 if l same no matter if l computed,',
940 $ / ' 7 = 0 if r same no matter if r computed,', / 1X )
941 9992 FORMAT( ' matrix order=', I5, ', type=', I2, ', seed=',
942 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
943 9991 FORMAT( ' matrix order=', I5, ', type=', I2, ', seed=',
944 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 )
945*
946* End of CDRGEV3
947*
948 END
#define alpha
Definition eval.h:35
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine cggev3(jobvl, jobvr, n, a, lda, b, ldb, alpha, beta, vl, ldvl, vr, ldvr, work, lwork, rwork, info)
CGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (...
Definition cggev3.f:216
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 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
subroutine cdrgev3(nsizes, nn, ntypes, dotype, iseed, thresh, nounit, a, lda, b, s, t, q, ldq, z, qe, ldqe, alpha, beta, alpha1, beta1, work, lwork, rwork, result, info)
CDRGEV3
Definition cdrgev3.f:399
subroutine clatm4(itype, n, nz1, nz2, rsign, amagn, rcond, triang, idist, iseed, a, lda)
CLATM4
Definition clatm4.f:171
subroutine cget52(left, n, a, lda, b, ldb, e, lde, alpha, beta, work, rwork, result)
CGET52
Definition cget52.f:161
#define seed()
Definition macros.h:43
#define min(a, b)
Definition macros.h:20
#define max(a, b)
Definition macros.h:21