LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
cchkhs.f
Go to the documentation of this file.
1*> \brief \b CCHKHS
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 CCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12* NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, W1,
13* W3, EVECTL, EVECTR, EVECTY, EVECTX, UU, TAU,
14* WORK, NWORK, RWORK, IWORK, SELECT, RESULT,
15* INFO )
16*
17* .. Scalar Arguments ..
18* INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK
19* REAL THRESH
20* ..
21* .. Array Arguments ..
22* LOGICAL DOTYPE( * ), SELECT( * )
23* INTEGER ISEED( 4 ), IWORK( * ), NN( * )
24* REAL RESULT( 14 ), RWORK( * )
25* COMPLEX A( LDA, * ), EVECTL( LDU, * ),
26* \$ EVECTR( LDU, * ), EVECTX( LDU, * ),
27* \$ EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ),
28* \$ T2( LDA, * ), TAU( * ), U( LDU, * ),
29* \$ UU( LDU, * ), UZ( LDU, * ), W1( * ), W3( * ),
30* \$ WORK( * ), Z( LDU, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> CCHKHS checks the nonsymmetric eigenvalue problem routines.
40*>
41*> CGEHRD factors A as U H U' , where ' means conjugate
42*> transpose, H is hessenberg, and U is unitary.
43*>
44*> CUNGHR generates the unitary matrix U.
45*>
46*> CUNMHR multiplies a matrix by the unitary matrix U.
47*>
48*> CHSEQR factors H as Z T Z' , where Z is unitary and T
49*> is upper triangular. It also computes the eigenvalues,
50*> w(1), ..., w(n); we define a diagonal matrix W whose
51*> (diagonal) entries are the eigenvalues.
52*>
53*> CTREVC computes the left eigenvector matrix L and the
54*> right eigenvector matrix R for the matrix T. The
55*> columns of L are the complex conjugates of the left
56*> eigenvectors of T. The columns of R are the right
57*> eigenvectors of T. L is lower triangular, and R is
58*> upper triangular.
59*>
60*> CHSEIN computes the left eigenvector matrix Y and the
61*> right eigenvector matrix X for the matrix H. The
62*> columns of Y are the complex conjugates of the left
63*> eigenvectors of H. The columns of X are the right
64*> eigenvectors of H. Y is lower triangular, and X is
65*> upper triangular.
66*>
67*> When CCHKHS is called, a number of matrix "sizes" ("n's") and a
68*> number of matrix "types" are specified. For each size ("n")
69*> and each type of matrix, one matrix will be generated and used
70*> to test the nonsymmetric eigenroutines. For each matrix, 14
71*> tests will be performed:
72*>
73*> (1) | A - U H U**H | / ( |A| n ulp )
74*>
75*> (2) | I - UU**H | / ( n ulp )
76*>
77*> (3) | H - Z T Z**H | / ( |H| n ulp )
78*>
79*> (4) | I - ZZ**H | / ( n ulp )
80*>
81*> (5) | A - UZ H (UZ)**H | / ( |A| n ulp )
82*>
83*> (6) | I - UZ (UZ)**H | / ( n ulp )
84*>
85*> (7) | T(Z computed) - T(Z not computed) | / ( |T| ulp )
86*>
87*> (8) | W(Z computed) - W(Z not computed) | / ( |W| ulp )
88*>
89*> (9) | TR - RW | / ( |T| |R| ulp )
90*>
91*> (10) | L**H T - W**H L | / ( |T| |L| ulp )
92*>
93*> (11) | HX - XW | / ( |H| |X| ulp )
94*>
95*> (12) | Y**H H - W**H Y | / ( |H| |Y| ulp )
96*>
97*> (13) | AX - XW | / ( |A| |X| ulp )
98*>
99*> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp )
100*>
101*> The "sizes" are specified by an array NN(1:NSIZES); the value of
102*> each element NN(j) specifies one size.
103*> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
104*> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
105*> Currently, the list of possible types is:
106*>
107*> (1) The zero matrix.
108*> (2) The identity matrix.
109*> (3) A (transposed) Jordan block, with 1's on the diagonal.
110*>
111*> (4) A diagonal matrix with evenly spaced entries
112*> 1, ..., ULP and random complex angles.
113*> (ULP = (first number larger than 1) - 1 )
114*> (5) A diagonal matrix with geometrically spaced entries
115*> 1, ..., ULP and random complex angles.
116*> (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
117*> and random complex angles.
118*>
119*> (7) Same as (4), but multiplied by SQRT( overflow threshold )
120*> (8) Same as (4), but multiplied by SQRT( underflow threshold )
121*>
122*> (9) A matrix of the form U' T U, where U is unitary and
123*> T has evenly spaced entries 1, ..., ULP with random complex
124*> angles on the diagonal and random O(1) entries in the upper
125*> triangle.
126*>
127*> (10) A matrix of the form U' T U, where U is unitary and
128*> T has geometrically spaced entries 1, ..., ULP with random
129*> complex angles on the diagonal and random O(1) entries in
130*> the upper triangle.
131*>
132*> (11) A matrix of the form U' T U, where U is unitary and
133*> T has "clustered" entries 1, ULP,..., ULP with random
134*> complex angles on the diagonal and random O(1) entries in
135*> the upper triangle.
136*>
137*> (12) A matrix of the form U' T U, where U is unitary and
138*> T has complex eigenvalues randomly chosen from
139*> ULP < |z| < 1 and random O(1) entries in the upper
140*> triangle.
141*>
142*> (13) A matrix of the form X' T X, where X has condition
143*> SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
144*> with random complex angles on the diagonal and random O(1)
145*> entries in the upper triangle.
146*>
147*> (14) A matrix of the form X' T X, where X has condition
148*> SQRT( ULP ) and T has geometrically spaced entries
149*> 1, ..., ULP with random complex angles on the diagonal
150*> and random O(1) entries in the upper triangle.
151*>
152*> (15) A matrix of the form X' T X, where X has condition
153*> SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
154*> with random complex angles on the diagonal and random O(1)
155*> entries in the upper triangle.
156*>
157*> (16) A matrix of the form X' T X, where X has condition
158*> SQRT( ULP ) and T has complex eigenvalues randomly chosen
159*> from ULP < |z| < 1 and random O(1) entries in the upper
160*> triangle.
161*>
162*> (17) Same as (16), but multiplied by SQRT( overflow threshold )
163*> (18) Same as (16), but multiplied by SQRT( underflow threshold )
164*>
165*> (19) Nonsymmetric matrix with random entries chosen from |z| < 1
166*> (20) Same as (19), but multiplied by SQRT( overflow threshold )
167*> (21) Same as (19), but multiplied by SQRT( underflow threshold )
168*> \endverbatim
169*
170* Arguments:
171* ==========
172*
173*> \verbatim
174*> NSIZES - INTEGER
175*> The number of sizes of matrices to use. If it is zero,
176*> CCHKHS does nothing. It must be at least zero.
177*> Not modified.
178*>
179*> NN - INTEGER array, dimension (NSIZES)
180*> An array containing the sizes to be used for the matrices.
181*> Zero values will be skipped. The values must be at least
182*> zero.
183*> Not modified.
184*>
185*> NTYPES - INTEGER
186*> The number of elements in DOTYPE. If it is zero, CCHKHS
187*> does nothing. It must be at least zero. If it is MAXTYP+1
188*> and NSIZES is 1, then an additional type, MAXTYP+1 is
189*> defined, which is to use whatever matrix is in A. This
190*> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
191*> DOTYPE(MAXTYP+1) is .TRUE. .
192*> Not modified.
193*>
194*> DOTYPE - LOGICAL array, dimension (NTYPES)
195*> If DOTYPE(j) is .TRUE., then for each size in NN a
196*> matrix of that size and of type j will be generated.
197*> If NTYPES is smaller than the maximum number of types
198*> defined (PARAMETER MAXTYP), then types NTYPES+1 through
199*> MAXTYP will not be generated. If NTYPES is larger
200*> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
201*> will be ignored.
202*> Not modified.
203*>
204*> ISEED - INTEGER array, dimension (4)
205*> On entry ISEED specifies the seed of the random number
206*> generator. The array elements should be between 0 and 4095;
207*> if not they will be reduced mod 4096. Also, ISEED(4) must
208*> be odd. The random number generator uses a linear
209*> congruential sequence limited to small integers, and so
210*> should produce machine independent random numbers. The
211*> values of ISEED are changed on exit, and can be used in the
212*> next call to CCHKHS to continue the same random number
213*> sequence.
214*> Modified.
215*>
216*> THRESH - REAL
217*> A test will count as "failed" if the "error", computed as
218*> described above, exceeds THRESH. Note that the error
219*> is scaled to be O(1), so THRESH should be a reasonably
220*> small multiple of 1, e.g., 10 or 100. In particular,
221*> it should not depend on the precision (single vs. double)
222*> or the size of the matrix. It must be at least zero.
223*> Not modified.
224*>
225*> NOUNIT - INTEGER
226*> The FORTRAN unit number for printing out error messages
227*> (e.g., if a routine returns IINFO not equal to 0.)
228*> Not modified.
229*>
230*> A - COMPLEX array, dimension (LDA,max(NN))
231*> Used to hold the matrix whose eigenvalues are to be
232*> computed. On exit, A contains the last matrix actually
233*> used.
234*> Modified.
235*>
236*> LDA - INTEGER
237*> The leading dimension of A, H, T1 and T2. It must be at
238*> least 1 and at least max( NN ).
239*> Not modified.
240*>
241*> H - COMPLEX array, dimension (LDA,max(NN))
242*> The upper hessenberg matrix computed by CGEHRD. On exit,
243*> H contains the Hessenberg form of the matrix in A.
244*> Modified.
245*>
246*> T1 - COMPLEX array, dimension (LDA,max(NN))
247*> The Schur (="quasi-triangular") matrix computed by CHSEQR
248*> if Z is computed. On exit, T1 contains the Schur form of
249*> the matrix in A.
250*> Modified.
251*>
252*> T2 - COMPLEX array, dimension (LDA,max(NN))
253*> The Schur matrix computed by CHSEQR when Z is not computed.
254*> This should be identical to T1.
255*> Modified.
256*>
257*> LDU - INTEGER
258*> The leading dimension of U, Z, UZ and UU. It must be at
259*> least 1 and at least max( NN ).
260*> Not modified.
261*>
262*> U - COMPLEX array, dimension (LDU,max(NN))
263*> The unitary matrix computed by CGEHRD.
264*> Modified.
265*>
266*> Z - COMPLEX array, dimension (LDU,max(NN))
267*> The unitary matrix computed by CHSEQR.
268*> Modified.
269*>
270*> UZ - COMPLEX array, dimension (LDU,max(NN))
271*> The product of U times Z.
272*> Modified.
273*>
274*> W1 - COMPLEX array, dimension (max(NN))
275*> The eigenvalues of A, as computed by a full Schur
276*> decomposition H = Z T Z'. On exit, W1 contains the
277*> eigenvalues of the matrix in A.
278*> Modified.
279*>
280*> W3 - COMPLEX array, dimension (max(NN))
281*> The eigenvalues of A, as computed by a partial Schur
282*> decomposition (Z not computed, T only computed as much
283*> as is necessary for determining eigenvalues). On exit,
284*> W3 contains the eigenvalues of the matrix in A, possibly
285*> perturbed by CHSEIN.
286*> Modified.
287*>
288*> EVECTL - COMPLEX array, dimension (LDU,max(NN))
289*> The conjugate transpose of the (upper triangular) left
290*> eigenvector matrix for the matrix in T1.
291*> Modified.
292*>
293*> EVECTR - COMPLEX array, dimension (LDU,max(NN))
294*> The (upper triangular) right eigenvector matrix for the
295*> matrix in T1.
296*> Modified.
297*>
298*> EVECTY - COMPLEX array, dimension (LDU,max(NN))
299*> The conjugate transpose of the left eigenvector matrix
300*> for the matrix in H.
301*> Modified.
302*>
303*> EVECTX - COMPLEX array, dimension (LDU,max(NN))
304*> The right eigenvector matrix for the matrix in H.
305*> Modified.
306*>
307*> UU - COMPLEX array, dimension (LDU,max(NN))
308*> Details of the unitary matrix computed by CGEHRD.
309*> Modified.
310*>
311*> TAU - COMPLEX array, dimension (max(NN))
312*> Further details of the unitary matrix computed by CGEHRD.
313*> Modified.
314*>
315*> WORK - COMPLEX array, dimension (NWORK)
316*> Workspace.
317*> Modified.
318*>
319*> NWORK - INTEGER
320*> The number of entries in WORK. NWORK >= 4*NN(j)*NN(j) + 2.
321*>
322*> RWORK - REAL array, dimension (max(NN))
323*> Workspace. Could be equivalenced to IWORK, but not SELECT.
324*> Modified.
325*>
326*> IWORK - INTEGER array, dimension (max(NN))
327*> Workspace.
328*> Modified.
329*>
330*> SELECT - LOGICAL array, dimension (max(NN))
331*> Workspace. Could be equivalenced to IWORK, but not RWORK.
332*> Modified.
333*>
334*> RESULT - REAL array, dimension (14)
335*> The values computed by the fourteen tests described above.
336*> The values are currently limited to 1/ulp, to avoid
337*> overflow.
338*> Modified.
339*>
340*> INFO - INTEGER
341*> If 0, then everything ran OK.
342*> -1: NSIZES < 0
343*> -2: Some NN(j) < 0
344*> -3: NTYPES < 0
345*> -6: THRESH < 0
346*> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
347*> -14: LDU < 1 or LDU < NMAX.
348*> -26: NWORK too small.
349*> If CLATMR, CLATMS, or CLATME returns an error code, the
350*> absolute value of it is returned.
351*> If 1, then CHSEQR could not find all the shifts.
352*> If 2, then the EISPACK code (for small blocks) failed.
353*> If >2, then 30*N iterations were not enough to find an
354*> eigenvalue or to decompose the problem.
355*> Modified.
356*>
357*>-----------------------------------------------------------------------
358*>
359*> Some Local Variables and Parameters:
360*> ---- ----- --------- --- ----------
361*>
362*> ZERO, ONE Real 0 and 1.
363*> MAXTYP The number of types defined.
364*> MTEST The number of tests defined: care must be taken
365*> that (1) the size of RESULT, (2) the number of
366*> tests actually performed, and (3) MTEST agree.
367*> NTEST The number of tests performed on this matrix
368*> so far. This should be less than MTEST, and
369*> equal to it by the last test. It will be less
370*> if any of the routines being tested indicates
371*> that it could not compute the matrices that
372*> would be tested.
373*> NMAX Largest value in NN.
374*> NMATS The number of matrices generated so far.
375*> NERRS The number of tests which have exceeded THRESH
376*> so far (computed by SLAFTS).
377*> COND, CONDS,
378*> IMODE Values to be passed to the matrix generators.
379*> ANORM Norm of A; passed to matrix generators.
380*>
381*> OVFL, UNFL Overflow and underflow thresholds.
382*> ULP, ULPINV Finest relative precision and its inverse.
383*> RTOVFL, RTUNFL,
384*> RTULP, RTULPI Square roots of the previous 4 values.
385*>
386*> The following four arrays decode JTYPE:
387*> KTYPE(j) The general type (1-10) for type "j".
388*> KMODE(j) The MODE value to be passed to the matrix
389*> generator for type "j".
390*> KMAGN(j) The order of magnitude ( O(1),
391*> O(overflow^(1/2) ), O(underflow^(1/2) )
392*> KCONDS(j) Selects whether CONDS is to be 1 or
393*> 1/sqrt(ulp). (0 means irrelevant.)
394*> \endverbatim
395*
396* Authors:
397* ========
398*
399*> \author Univ. of Tennessee
400*> \author Univ. of California Berkeley
401*> \author Univ. of Colorado Denver
402*> \author NAG Ltd.
403*
404*> \ingroup complex_eig
405*
406* =====================================================================
407 SUBROUTINE cchkhs( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
408 \$ NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, W1,
409 \$ W3, EVECTL, EVECTR, EVECTY, EVECTX, UU, TAU,
410 \$ WORK, NWORK, RWORK, IWORK, SELECT, RESULT,
411 \$ INFO )
412*
413* -- LAPACK test routine --
414* -- LAPACK is a software package provided by Univ. of Tennessee, --
415* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
416*
417* .. Scalar Arguments ..
418 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK
419 REAL THRESH
420* ..
421* .. Array Arguments ..
422 LOGICAL DOTYPE( * ), SELECT( * )
423 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
424 REAL RESULT( 14 ), RWORK( * )
425 COMPLEX A( LDA, * ), EVECTL( LDU, * ),
426 \$ evectr( ldu, * ), evectx( ldu, * ),
427 \$ evecty( ldu, * ), h( lda, * ), t1( lda, * ),
428 \$ t2( lda, * ), tau( * ), u( ldu, * ),
429 \$ uu( ldu, * ), uz( ldu, * ), w1( * ), w3( * ),
430 \$ work( * ), z( ldu, * )
431* ..
432*
433* =====================================================================
434*
435* .. Parameters ..
436 REAL ZERO, ONE
437 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
438 COMPLEX CZERO, CONE
439 PARAMETER ( CZERO = ( 0.0e+0, 0.0e+0 ),
440 \$ cone = ( 1.0e+0, 0.0e+0 ) )
441 INTEGER MAXTYP
442 parameter( maxtyp = 21 )
443* ..
444* .. Local Scalars ..
446 INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL,
447 \$ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS,
448 \$ NMATS, NMAX, NTEST, NTESTT
449 REAL ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP,
450 \$ rtulpi, rtunfl, temp1, temp2, ulp, ulpinv, unfl
451* ..
452* .. Local Arrays ..
453 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
454 \$ KMAGN( MAXTYP ), KMODE( MAXTYP ),
455 \$ KTYPE( MAXTYP )
456 REAL DUMMA( 4 )
457 COMPLEX CDUMMA( 4 )
458* ..
459* .. External Functions ..
460 REAL SLAMCH
461 EXTERNAL slamch
462* ..
463* .. External Subroutines ..
464 EXTERNAL ccopy, cgehrd, cgemm, cget10, cget22, chsein,
467 \$ slasum, xerbla
468* ..
469* .. Intrinsic Functions ..
470 INTRINSIC abs, max, min, real, sqrt
471* ..
472* .. Data statements ..
473 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
474 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
475 \$ 3, 1, 2, 3 /
476 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
477 \$ 1, 5, 5, 5, 4, 3, 1 /
478 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
479* ..
480* .. Executable Statements ..
481*
482* Check for errors
483*
484 ntestt = 0
485 info = 0
486*
488 nmax = 0
489 DO 10 j = 1, nsizes
490 nmax = max( nmax, nn( j ) )
491 IF( nn( j ).LT.0 )
493 10 CONTINUE
494*
495* Check for errors
496*
497 IF( nsizes.LT.0 ) THEN
498 info = -1
499 ELSE IF( badnn ) THEN
500 info = -2
501 ELSE IF( ntypes.LT.0 ) THEN
502 info = -3
503 ELSE IF( thresh.LT.zero ) THEN
504 info = -6
505 ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
506 info = -9
507 ELSE IF( ldu.LE.1 .OR. ldu.LT.nmax ) THEN
508 info = -14
509 ELSE IF( 4*nmax*nmax+2.GT.nwork ) THEN
510 info = -26
511 END IF
512*
513 IF( info.NE.0 ) THEN
514 CALL xerbla( 'CCHKHS', -info )
515 RETURN
516 END IF
517*
518* Quick return if possible
519*
520 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
521 \$ RETURN
522*
523* More important constants
524*
525 unfl = slamch( 'Safe minimum' )
526 ovfl = slamch( 'Overflow' )
527 CALL slabad( unfl, ovfl )
528 ulp = slamch( 'Epsilon' )*slamch( 'Base' )
529 ulpinv = one / ulp
530 rtunfl = sqrt( unfl )
531 rtovfl = sqrt( ovfl )
532 rtulp = sqrt( ulp )
533 rtulpi = one / rtulp
534*
535* Loop over sizes, types
536*
537 nerrs = 0
538 nmats = 0
539*
540 DO 260 jsize = 1, nsizes
541 n = nn( jsize )
542 IF( n.EQ.0 )
543 \$ GO TO 260
544 n1 = max( 1, n )
545 aninv = one / real( n1 )
546*
547 IF( nsizes.NE.1 ) THEN
548 mtypes = min( maxtyp, ntypes )
549 ELSE
550 mtypes = min( maxtyp+1, ntypes )
551 END IF
552*
553 DO 250 jtype = 1, mtypes
554 IF( .NOT.dotype( jtype ) )
555 \$ GO TO 250
556 nmats = nmats + 1
557 ntest = 0
558*
559* Save ISEED in case of an error.
560*
561 DO 20 j = 1, 4
562 ioldsd( j ) = iseed( j )
563 20 CONTINUE
564*
565* Initialize RESULT
566*
567 DO 30 j = 1, 14
568 result( j ) = zero
569 30 CONTINUE
570*
571* Compute "A"
572*
573* Control parameters:
574*
575* KMAGN KCONDS KMODE KTYPE
576* =1 O(1) 1 clustered 1 zero
577* =2 large large clustered 2 identity
578* =3 small exponential Jordan
579* =4 arithmetic diagonal, (w/ eigenvalues)
580* =5 random log hermitian, w/ eigenvalues
581* =6 random general, w/ eigenvalues
582* =7 random diagonal
583* =8 random hermitian
584* =9 random general
585* =10 random triangular
586*
587 IF( mtypes.GT.maxtyp )
588 \$ GO TO 100
589*
590 itype = ktype( jtype )
591 imode = kmode( jtype )
592*
593* Compute norm
594*
595 GO TO ( 40, 50, 60 )kmagn( jtype )
596*
597 40 CONTINUE
598 anorm = one
599 GO TO 70
600*
601 50 CONTINUE
602 anorm = ( rtovfl*ulp )*aninv
603 GO TO 70
604*
605 60 CONTINUE
606 anorm = rtunfl*n*ulpinv
607 GO TO 70
608*
609 70 CONTINUE
610*
611 CALL claset( 'Full', lda, n, czero, czero, a, lda )
612 iinfo = 0
613 cond = ulpinv
614*
615* Special Matrices
616*
617 IF( itype.EQ.1 ) THEN
618*
619* Zero
620*
621 iinfo = 0
622 ELSE IF( itype.EQ.2 ) THEN
623*
624* Identity
625*
626 DO 80 jcol = 1, n
627 a( jcol, jcol ) = anorm
628 80 CONTINUE
629*
630 ELSE IF( itype.EQ.3 ) THEN
631*
632* Jordan Block
633*
634 DO 90 jcol = 1, n
635 a( jcol, jcol ) = anorm
636 IF( jcol.GT.1 )
637 \$ a( jcol, jcol-1 ) = one
638 90 CONTINUE
639*
640 ELSE IF( itype.EQ.4 ) THEN
641*
642* Diagonal Matrix, [Eigen]values Specified
643*
644 CALL clatmr( n, n, 'D', iseed, 'N', work, imode, cond,
645 \$ cone, 'T', 'N', work( n+1 ), 1, one,
646 \$ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
647 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
648*
649 ELSE IF( itype.EQ.5 ) THEN
650*
651* Hermitian, eigenvalues specified
652*
653 CALL clatms( n, n, 'D', iseed, 'H', rwork, imode, cond,
654 \$ anorm, n, n, 'N', a, lda, work, iinfo )
655*
656 ELSE IF( itype.EQ.6 ) THEN
657*
658* General, eigenvalues specified
659*
660 IF( kconds( jtype ).EQ.1 ) THEN
661 conds = one
662 ELSE IF( kconds( jtype ).EQ.2 ) THEN
663 conds = rtulpi
664 ELSE
665 conds = zero
666 END IF
667*
668 CALL clatme( n, 'D', iseed, work, imode, cond, cone,
669 \$ 'T', 'T', 'T', rwork, 4, conds, n, n, anorm,
670 \$ a, lda, work( n+1 ), iinfo )
671*
672 ELSE IF( itype.EQ.7 ) THEN
673*
674* Diagonal, random eigenvalues
675*
676 CALL clatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
677 \$ 'T', 'N', work( n+1 ), 1, one,
678 \$ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
679 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
680*
681 ELSE IF( itype.EQ.8 ) THEN
682*
683* Hermitian, random eigenvalues
684*
685 CALL clatmr( n, n, 'D', iseed, 'H', work, 6, one, cone,
686 \$ 'T', 'N', work( n+1 ), 1, one,
687 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
688 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
689*
690 ELSE IF( itype.EQ.9 ) THEN
691*
692* General, random eigenvalues
693*
694 CALL clatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
695 \$ 'T', 'N', work( n+1 ), 1, one,
696 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
697 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
698*
699 ELSE IF( itype.EQ.10 ) THEN
700*
701* Triangular, random eigenvalues
702*
703 CALL clatmr( n, n, 'D', iseed, 'N', work, 6, one, cone,
704 \$ 'T', 'N', work( n+1 ), 1, one,
705 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
706 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
707*
708 ELSE
709*
710 iinfo = 1
711 END IF
712*
713 IF( iinfo.NE.0 ) THEN
714 WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
715 \$ ioldsd
716 info = abs( iinfo )
717 RETURN
718 END IF
719*
720 100 CONTINUE
721*
722* Call CGEHRD to compute H and U, do tests.
723*
724 CALL clacpy( ' ', n, n, a, lda, h, lda )
725 ntest = 1
726*
727 ilo = 1
728 ihi = n
729*
730 CALL cgehrd( n, ilo, ihi, h, lda, work, work( n+1 ),
731 \$ nwork-n, iinfo )
732*
733 IF( iinfo.NE.0 ) THEN
734 result( 1 ) = ulpinv
735 WRITE( nounit, fmt = 9999 )'CGEHRD', iinfo, n, jtype,
736 \$ ioldsd
737 info = abs( iinfo )
738 GO TO 240
739 END IF
740*
741 DO 120 j = 1, n - 1
742 uu( j+1, j ) = czero
743 DO 110 i = j + 2, n
744 u( i, j ) = h( i, j )
745 uu( i, j ) = h( i, j )
746 h( i, j ) = czero
747 110 CONTINUE
748 120 CONTINUE
749 CALL ccopy( n-1, work, 1, tau, 1 )
750 CALL cunghr( n, ilo, ihi, u, ldu, work, work( n+1 ),
751 \$ nwork-n, iinfo )
752 ntest = 2
753*
754 CALL chst01( n, ilo, ihi, a, lda, h, lda, u, ldu, work,
755 \$ nwork, rwork, result( 1 ) )
756*
757* Call CHSEQR to compute T1, T2 and Z, do tests.
758*
759* Eigenvalues only (W3)
760*
761 CALL clacpy( ' ', n, n, h, lda, t2, lda )
762 ntest = 3
763 result( 3 ) = ulpinv
764*
765 CALL chseqr( 'E', 'N', n, ilo, ihi, t2, lda, w3, uz, ldu,
766 \$ work, nwork, iinfo )
767 IF( iinfo.NE.0 ) THEN
768 WRITE( nounit, fmt = 9999 )'CHSEQR(E)', iinfo, n, jtype,
769 \$ ioldsd
770 IF( iinfo.LE.n+2 ) THEN
771 info = abs( iinfo )
772 GO TO 240
773 END IF
774 END IF
775*
776* Eigenvalues (W1) and Full Schur Form (T2)
777*
778 CALL clacpy( ' ', n, n, h, lda, t2, lda )
779*
780 CALL chseqr( 'S', 'N', n, ilo, ihi, t2, lda, w1, uz, ldu,
781 \$ work, nwork, iinfo )
782 IF( iinfo.NE.0 .AND. iinfo.LE.n+2 ) THEN
783 WRITE( nounit, fmt = 9999 )'CHSEQR(S)', iinfo, n, jtype,
784 \$ ioldsd
785 info = abs( iinfo )
786 GO TO 240
787 END IF
788*
789* Eigenvalues (W1), Schur Form (T1), and Schur Vectors (UZ)
790*
791 CALL clacpy( ' ', n, n, h, lda, t1, lda )
792 CALL clacpy( ' ', n, n, u, ldu, uz, ldu )
793*
794 CALL chseqr( 'S', 'V', n, ilo, ihi, t1, lda, w1, uz, ldu,
795 \$ work, nwork, iinfo )
796 IF( iinfo.NE.0 .AND. iinfo.LE.n+2 ) THEN
797 WRITE( nounit, fmt = 9999 )'CHSEQR(V)', iinfo, n, jtype,
798 \$ ioldsd
799 info = abs( iinfo )
800 GO TO 240
801 END IF
802*
803* Compute Z = U' UZ
804*
805 CALL cgemm( 'C', 'N', n, n, n, cone, u, ldu, uz, ldu, czero,
806 \$ z, ldu )
807 ntest = 8
808*
809* Do Tests 3: | H - Z T Z' | / ( |H| n ulp )
810* and 4: | I - Z Z' | / ( n ulp )
811*
812 CALL chst01( n, ilo, ihi, h, lda, t1, lda, z, ldu, work,
813 \$ nwork, rwork, result( 3 ) )
814*
815* Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp )
816* and 6: | I - UZ (UZ)' | / ( n ulp )
817*
818 CALL chst01( n, ilo, ihi, a, lda, t1, lda, uz, ldu, work,
819 \$ nwork, rwork, result( 5 ) )
820*
821* Do Test 7: | T2 - T1 | / ( |T| n ulp )
822*
823 CALL cget10( n, n, t2, lda, t1, lda, work, rwork,
824 \$ result( 7 ) )
825*
826* Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp )
827*
828 temp1 = zero
829 temp2 = zero
830 DO 130 j = 1, n
831 temp1 = max( temp1, abs( w1( j ) ), abs( w3( j ) ) )
832 temp2 = max( temp2, abs( w1( j )-w3( j ) ) )
833 130 CONTINUE
834*
835 result( 8 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
836*
837* Compute the Left and Right Eigenvectors of T
838*
839* Compute the Right eigenvector Matrix:
840*
841 ntest = 9
842 result( 9 ) = ulpinv
843*
844* Select every other eigenvector
845*
846 DO 140 j = 1, n
847 SELECT( j ) = .false.
848 140 CONTINUE
849 DO 150 j = 1, n, 2
850 SELECT( j ) = .true.
851 150 CONTINUE
852 CALL ctrevc( 'Right', 'All', SELECT, n, t1, lda, cdumma,
853 \$ ldu, evectr, ldu, n, in, work, rwork, iinfo )
854 IF( iinfo.NE.0 ) THEN
855 WRITE( nounit, fmt = 9999 )'CTREVC(R,A)', iinfo, n,
856 \$ jtype, ioldsd
857 info = abs( iinfo )
858 GO TO 240
859 END IF
860*
861* Test 9: | TR - RW | / ( |T| |R| ulp )
862*
863 CALL cget22( 'N', 'N', 'N', n, t1, lda, evectr, ldu, w1,
864 \$ work, rwork, dumma( 1 ) )
865 result( 9 ) = dumma( 1 )
866 IF( dumma( 2 ).GT.thresh ) THEN
867 WRITE( nounit, fmt = 9998 )'Right', 'CTREVC',
868 \$ dumma( 2 ), n, jtype, ioldsd
869 END IF
870*
871* Compute selected right eigenvectors and confirm that
872* they agree with previous right eigenvectors
873*
874 CALL ctrevc( 'Right', 'Some', SELECT, n, t1, lda, cdumma,
875 \$ ldu, evectl, ldu, n, in, work, rwork, iinfo )
876 IF( iinfo.NE.0 ) THEN
877 WRITE( nounit, fmt = 9999 )'CTREVC(R,S)', iinfo, n,
878 \$ jtype, ioldsd
879 info = abs( iinfo )
880 GO TO 240
881 END IF
882*
883 k = 1
884 match = .true.
885 DO 170 j = 1, n
886 IF( SELECT( j ) ) THEN
887 DO 160 jj = 1, n
888 IF( evectr( jj, j ).NE.evectl( jj, k ) ) THEN
889 match = .false.
890 GO TO 180
891 END IF
892 160 CONTINUE
893 k = k + 1
894 END IF
895 170 CONTINUE
896 180 CONTINUE
897 IF( .NOT.match )
898 \$ WRITE( nounit, fmt = 9997 )'Right', 'CTREVC', n, jtype,
899 \$ ioldsd
900*
901* Compute the Left eigenvector Matrix:
902*
903 ntest = 10
904 result( 10 ) = ulpinv
905 CALL ctrevc( 'Left', 'All', SELECT, n, t1, lda, evectl, ldu,
906 \$ cdumma, ldu, n, in, work, rwork, iinfo )
907 IF( iinfo.NE.0 ) THEN
908 WRITE( nounit, fmt = 9999 )'CTREVC(L,A)', iinfo, n,
909 \$ jtype, ioldsd
910 info = abs( iinfo )
911 GO TO 240
912 END IF
913*
914* Test 10: | LT - WL | / ( |T| |L| ulp )
915*
916 CALL cget22( 'C', 'N', 'C', n, t1, lda, evectl, ldu, w1,
917 \$ work, rwork, dumma( 3 ) )
918 result( 10 ) = dumma( 3 )
919 IF( dumma( 4 ).GT.thresh ) THEN
920 WRITE( nounit, fmt = 9998 )'Left', 'CTREVC', dumma( 4 ),
921 \$ n, jtype, ioldsd
922 END IF
923*
924* Compute selected left eigenvectors and confirm that
925* they agree with previous left eigenvectors
926*
927 CALL ctrevc( 'Left', 'Some', SELECT, n, t1, lda, evectr,
928 \$ ldu, cdumma, ldu, n, in, work, rwork, iinfo )
929 IF( iinfo.NE.0 ) THEN
930 WRITE( nounit, fmt = 9999 )'CTREVC(L,S)', iinfo, n,
931 \$ jtype, ioldsd
932 info = abs( iinfo )
933 GO TO 240
934 END IF
935*
936 k = 1
937 match = .true.
938 DO 200 j = 1, n
939 IF( SELECT( j ) ) THEN
940 DO 190 jj = 1, n
941 IF( evectl( jj, j ).NE.evectr( jj, k ) ) THEN
942 match = .false.
943 GO TO 210
944 END IF
945 190 CONTINUE
946 k = k + 1
947 END IF
948 200 CONTINUE
949 210 CONTINUE
950 IF( .NOT.match )
951 \$ WRITE( nounit, fmt = 9997 )'Left', 'CTREVC', n, jtype,
952 \$ ioldsd
953*
954* Call CHSEIN for Right eigenvectors of H, do test 11
955*
956 ntest = 11
957 result( 11 ) = ulpinv
958 DO 220 j = 1, n
959 SELECT( j ) = .true.
960 220 CONTINUE
961*
962 CALL chsein( 'Right', 'Qr', 'Ninitv', SELECT, n, h, lda, w3,
963 \$ cdumma, ldu, evectx, ldu, n1, in, work, rwork,
964 \$ iwork, iwork, iinfo )
965 IF( iinfo.NE.0 ) THEN
966 WRITE( nounit, fmt = 9999 )'CHSEIN(R)', iinfo, n, jtype,
967 \$ ioldsd
968 info = abs( iinfo )
969 IF( iinfo.LT.0 )
970 \$ GO TO 240
971 ELSE
972*
973* Test 11: | HX - XW | / ( |H| |X| ulp )
974*
975* (from inverse iteration)
976*
977 CALL cget22( 'N', 'N', 'N', n, h, lda, evectx, ldu, w3,
978 \$ work, rwork, dumma( 1 ) )
979 IF( dumma( 1 ).LT.ulpinv )
980 \$ result( 11 ) = dumma( 1 )*aninv
981 IF( dumma( 2 ).GT.thresh ) THEN
982 WRITE( nounit, fmt = 9998 )'Right', 'CHSEIN',
983 \$ dumma( 2 ), n, jtype, ioldsd
984 END IF
985 END IF
986*
987* Call CHSEIN for Left eigenvectors of H, do test 12
988*
989 ntest = 12
990 result( 12 ) = ulpinv
991 DO 230 j = 1, n
992 SELECT( j ) = .true.
993 230 CONTINUE
994*
995 CALL chsein( 'Left', 'Qr', 'Ninitv', SELECT, n, h, lda, w3,
996 \$ evecty, ldu, cdumma, ldu, n1, in, work, rwork,
997 \$ iwork, iwork, iinfo )
998 IF( iinfo.NE.0 ) THEN
999 WRITE( nounit, fmt = 9999 )'CHSEIN(L)', iinfo, n, jtype,
1000 \$ ioldsd
1001 info = abs( iinfo )
1002 IF( iinfo.LT.0 )
1003 \$ GO TO 240
1004 ELSE
1005*
1006* Test 12: | YH - WY | / ( |H| |Y| ulp )
1007*
1008* (from inverse iteration)
1009*
1010 CALL cget22( 'C', 'N', 'C', n, h, lda, evecty, ldu, w3,
1011 \$ work, rwork, dumma( 3 ) )
1012 IF( dumma( 3 ).LT.ulpinv )
1013 \$ result( 12 ) = dumma( 3 )*aninv
1014 IF( dumma( 4 ).GT.thresh ) THEN
1015 WRITE( nounit, fmt = 9998 )'Left', 'CHSEIN',
1016 \$ dumma( 4 ), n, jtype, ioldsd
1017 END IF
1018 END IF
1019*
1020* Call CUNMHR for Right eigenvectors of A, do test 13
1021*
1022 ntest = 13
1023 result( 13 ) = ulpinv
1024*
1025 CALL cunmhr( 'Left', 'No transpose', n, n, ilo, ihi, uu,
1026 \$ ldu, tau, evectx, ldu, work, nwork, iinfo )
1027 IF( iinfo.NE.0 ) THEN
1028 WRITE( nounit, fmt = 9999 )'CUNMHR(L)', iinfo, n, jtype,
1029 \$ ioldsd
1030 info = abs( iinfo )
1031 IF( iinfo.LT.0 )
1032 \$ GO TO 240
1033 ELSE
1034*
1035* Test 13: | AX - XW | / ( |A| |X| ulp )
1036*
1037* (from inverse iteration)
1038*
1039 CALL cget22( 'N', 'N', 'N', n, a, lda, evectx, ldu, w3,
1040 \$ work, rwork, dumma( 1 ) )
1041 IF( dumma( 1 ).LT.ulpinv )
1042 \$ result( 13 ) = dumma( 1 )*aninv
1043 END IF
1044*
1045* Call CUNMHR for Left eigenvectors of A, do test 14
1046*
1047 ntest = 14
1048 result( 14 ) = ulpinv
1049*
1050 CALL cunmhr( 'Left', 'No transpose', n, n, ilo, ihi, uu,
1051 \$ ldu, tau, evecty, ldu, work, nwork, iinfo )
1052 IF( iinfo.NE.0 ) THEN
1053 WRITE( nounit, fmt = 9999 )'CUNMHR(L)', iinfo, n, jtype,
1054 \$ ioldsd
1055 info = abs( iinfo )
1056 IF( iinfo.LT.0 )
1057 \$ GO TO 240
1058 ELSE
1059*
1060* Test 14: | YA - WY | / ( |A| |Y| ulp )
1061*
1062* (from inverse iteration)
1063*
1064 CALL cget22( 'C', 'N', 'C', n, a, lda, evecty, ldu, w3,
1065 \$ work, rwork, dumma( 3 ) )
1066 IF( dumma( 3 ).LT.ulpinv )
1067 \$ result( 14 ) = dumma( 3 )*aninv
1068 END IF
1069*
1070* End of Loop -- Check for RESULT(j) > THRESH
1071*
1072 240 CONTINUE
1073*
1074 ntestt = ntestt + ntest
1075 CALL slafts( 'CHS', n, n, jtype, ntest, result, ioldsd,
1076 \$ thresh, nounit, nerrs )
1077*
1078 250 CONTINUE
1079 260 CONTINUE
1080*
1081* Summary
1082*
1083 CALL slasum( 'CHS', nounit, nerrs, ntestt )
1084*
1085 RETURN
1086*
1087 9999 FORMAT( ' CCHKHS: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
1088 \$ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1089 9998 FORMAT( ' CCHKHS: ', a, ' Eigenvectors from ', a, ' incorrectly ',
1090 \$ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
1091 \$ 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5,
1092 \$ ')' )
1093 9997 FORMAT( ' CCHKHS: Selected ', a, ' Eigenvectors from ', a,
1094 \$ ' do not match other eigenvectors ', 9x, 'N=', i6,
1095 \$ ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1096*
1097* End of CCHKHS
1098*
1099 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine chst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RWORK, RESULT)
CHST01
Definition: chst01.f:140
subroutine cget10(M, N, A, LDA, B, LDB, WORK, RWORK, RESULT)
CGET10
Definition: cget10.f:99
subroutine cget22(TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, W, WORK, RWORK, RESULT)
CGET22
Definition: cget22.f:144
subroutine cchkhs(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, W1, W3, EVECTL, EVECTR, EVECTY, EVECTX, UU, TAU, WORK, NWORK, RWORK, IWORK, SELECT, RESULT, INFO)
CCHKHS
Definition: cchkhs.f:412
subroutine clatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
CLATMS
Definition: clatms.f:332
subroutine clatme(N, DIST, ISEED, D, MODE, COND, DMAX, RSIGN, UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, A, LDA, WORK, INFO)
CLATME
Definition: clatme.f:301
subroutine clatmr(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, RSIGN, GRADE, DL, MODEL, CONDL, DR, MODER, CONDR, PIVTNG, IPIVOT, KL, KU, SPARSE, ANORM, PACK, A, LDA, IWORK, INFO)
CLATMR
Definition: clatmr.f:490
subroutine cgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CGEHRD
Definition: cgehrd.f:167
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 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 cunmhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMHR
Definition: cunmhr.f:179
subroutine cunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CUNGHR
Definition: cunghr.f:126
subroutine chseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
CHSEQR
Definition: chseqr.f:299
subroutine chsein(SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL, IFAILR, INFO)
CHSEIN
Definition: chsein.f:245
subroutine ctrevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTREVC
Definition: ctrevc.f:218
subroutine slafts(TYPE, M, N, IMAT, NTESTS, RESULT, ISEED, THRESH, IOUNIT, IE)
SLAFTS
Definition: slafts.f:99
subroutine slasum(TYPE, IOUNIT, IE, NRUN)
SLASUM
Definition: slasum.f:41