LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zchkhb2stg.f
Go to the documentation of this file.
1*> \brief \b ZCHKHB2STG
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 ZCHKHB2STG( NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE,
12* ISEED, THRESH, NOUNIT, A, LDA, SD, SE, D1,
13* D2, D3, U, LDU, WORK, LWORK, RWORK RESULT,
14* INFO )
15*
16* .. Scalar Arguments ..
17* INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES,
18* $ NWDTHS
19* DOUBLE PRECISION THRESH
20* ..
21* .. Array Arguments ..
22* LOGICAL DOTYPE( * )
23* INTEGER ISEED( 4 ), KK( * ), NN( * )
24* DOUBLE PRECISION RESULT( * ), RWORK( * ), SD( * ), SE( * ),
25* $ D1( * ), D2( * ), D3( * )
26* COMPLEX*16 A( LDA, * ), U( LDU, * ), WORK( * )
27* ..
28*
29*
30*> \par Purpose:
31* =============
32*>
33*> \verbatim
34*>
35*> ZCHKHB2STG tests the reduction of a Hermitian band matrix to tridiagonal
36*> from, used with the Hermitian eigenvalue problem.
37*>
38*> ZHBTRD factors a Hermitian band matrix A as U S U* , where * means
39*> conjugate transpose, S is symmetric tridiagonal, and U is unitary.
40*> ZHBTRD can use either just the lower or just the upper triangle
41*> of A; ZCHKHB2STG checks both cases.
42*>
43*> ZHETRD_HB2ST factors a Hermitian band matrix A as U S U* ,
44*> where * means conjugate transpose, S is symmetric tridiagonal, and U is
45*> unitary. ZHETRD_HB2ST can use either just the lower or just
46*> the upper triangle of A; ZCHKHB2STG checks both cases.
47*>
48*> DSTEQR factors S as Z D1 Z'.
49*> D1 is the matrix of eigenvalues computed when Z is not computed
50*> and from the S resulting of DSBTRD "U" (used as reference for DSYTRD_SB2ST)
51*> D2 is the matrix of eigenvalues computed when Z is not computed
52*> and from the S resulting of DSYTRD_SB2ST "U".
53*> D3 is the matrix of eigenvalues computed when Z is not computed
54*> and from the S resulting of DSYTRD_SB2ST "L".
55*>
56*> When ZCHKHB2STG is called, a number of matrix "sizes" ("n's"), a number
57*> of bandwidths ("k's"), and a number of matrix "types" are
58*> specified. For each size ("n"), each bandwidth ("k") less than or
59*> equal to "n", and each type of matrix, one matrix will be generated
60*> and used to test the hermitian banded reduction routine. For each
61*> matrix, a number of tests will be performed:
62*>
63*> (1) | A - V S V* | / ( |A| n ulp ) computed by ZHBTRD with
64*> UPLO='U'
65*>
66*> (2) | I - UU* | / ( n ulp )
67*>
68*> (3) | A - V S V* | / ( |A| n ulp ) computed by ZHBTRD with
69*> UPLO='L'
70*>
71*> (4) | I - UU* | / ( n ulp )
72*>
73*> (5) | D1 - D2 | / ( |D1| ulp ) where D1 is computed by
74*> DSBTRD with UPLO='U' and
75*> D2 is computed by
76*> ZHETRD_HB2ST with UPLO='U'
77*>
78*> (6) | D1 - D3 | / ( |D1| ulp ) where D1 is computed by
79*> DSBTRD with UPLO='U' and
80*> D3 is computed by
81*> ZHETRD_HB2ST with UPLO='L'
82*>
83*> The "sizes" are specified by an array NN(1:NSIZES); the value of
84*> each element NN(j) specifies one size.
85*> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
86*> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
87*> Currently, the list of possible types is:
88*>
89*> (1) The zero matrix.
90*> (2) The identity matrix.
91*>
92*> (3) A diagonal matrix with evenly spaced entries
93*> 1, ..., ULP and random signs.
94*> (ULP = (first number larger than 1) - 1 )
95*> (4) A diagonal matrix with geometrically spaced entries
96*> 1, ..., ULP and random signs.
97*> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
98*> and random signs.
99*>
100*> (6) Same as (4), but multiplied by SQRT( overflow threshold )
101*> (7) Same as (4), but multiplied by SQRT( underflow threshold )
102*>
103*> (8) A matrix of the form U* D U, where U is unitary and
104*> D has evenly spaced entries 1, ..., ULP with random signs
105*> on the diagonal.
106*>
107*> (9) A matrix of the form U* D U, where U is unitary and
108*> D has geometrically spaced entries 1, ..., ULP with random
109*> signs on the diagonal.
110*>
111*> (10) A matrix of the form U* D U, where U is unitary and
112*> D has "clustered" entries 1, ULP,..., ULP with random
113*> signs on the diagonal.
114*>
115*> (11) Same as (8), but multiplied by SQRT( overflow threshold )
116*> (12) Same as (8), but multiplied by SQRT( underflow threshold )
117*>
118*> (13) Hermitian matrix with random entries chosen from (-1,1).
119*> (14) Same as (13), but multiplied by SQRT( overflow threshold )
120*> (15) Same as (13), but multiplied by SQRT( underflow threshold )
121*> \endverbatim
122*
123* Arguments:
124* ==========
125*
126*> \param[in] NSIZES
127*> \verbatim
128*> NSIZES is INTEGER
129*> The number of sizes of matrices to use. If it is zero,
130*> ZCHKHB2STG does nothing. It must be at least zero.
131*> \endverbatim
132*>
133*> \param[in] NN
134*> \verbatim
135*> NN is INTEGER array, dimension (NSIZES)
136*> An array containing the sizes to be used for the matrices.
137*> Zero values will be skipped. The values must be at least
138*> zero.
139*> \endverbatim
140*>
141*> \param[in] NWDTHS
142*> \verbatim
143*> NWDTHS is INTEGER
144*> The number of bandwidths to use. If it is zero,
145*> ZCHKHB2STG does nothing. It must be at least zero.
146*> \endverbatim
147*>
148*> \param[in] KK
149*> \verbatim
150*> KK is INTEGER array, dimension (NWDTHS)
151*> An array containing the bandwidths to be used for the band
152*> matrices. The values must be at least zero.
153*> \endverbatim
154*>
155*> \param[in] NTYPES
156*> \verbatim
157*> NTYPES is INTEGER
158*> The number of elements in DOTYPE. If it is zero, ZCHKHB2STG
159*> does nothing. It must be at least zero. If it is MAXTYP+1
160*> and NSIZES is 1, then an additional type, MAXTYP+1 is
161*> defined, which is to use whatever matrix is in A. This
162*> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
163*> DOTYPE(MAXTYP+1) is .TRUE. .
164*> \endverbatim
165*>
166*> \param[in] DOTYPE
167*> \verbatim
168*> DOTYPE is LOGICAL array, dimension (NTYPES)
169*> If DOTYPE(j) is .TRUE., then for each size in NN a
170*> matrix of that size and of type j will be generated.
171*> If NTYPES is smaller than the maximum number of types
172*> defined (PARAMETER MAXTYP), then types NTYPES+1 through
173*> MAXTYP will not be generated. If NTYPES is larger
174*> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
175*> will be ignored.
176*> \endverbatim
177*>
178*> \param[in,out] ISEED
179*> \verbatim
180*> ISEED is INTEGER array, dimension (4)
181*> On entry ISEED specifies the seed of the random number
182*> generator. The array elements should be between 0 and 4095;
183*> if not they will be reduced mod 4096. Also, ISEED(4) must
184*> be odd. The random number generator uses a linear
185*> congruential sequence limited to small integers, and so
186*> should produce machine independent random numbers. The
187*> values of ISEED are changed on exit, and can be used in the
188*> next call to ZCHKHB2STG to continue the same random number
189*> sequence.
190*> \endverbatim
191*>
192*> \param[in] THRESH
193*> \verbatim
194*> THRESH is DOUBLE PRECISION
195*> A test will count as "failed" if the "error", computed as
196*> described above, exceeds THRESH. Note that the error
197*> is scaled to be O(1), so THRESH should be a reasonably
198*> small multiple of 1, e.g., 10 or 100. In particular,
199*> it should not depend on the precision (single vs. double)
200*> or the size of the matrix. It must be at least zero.
201*> \endverbatim
202*>
203*> \param[in] NOUNIT
204*> \verbatim
205*> NOUNIT is INTEGER
206*> The FORTRAN unit number for printing out error messages
207*> (e.g., if a routine returns IINFO not equal to 0.)
208*> \endverbatim
209*>
210*> \param[in,out] A
211*> \verbatim
212*> A is COMPLEX*16 array, dimension
213*> (LDA, max(NN))
214*> Used to hold the matrix whose eigenvalues are to be
215*> computed.
216*> \endverbatim
217*>
218*> \param[in] LDA
219*> \verbatim
220*> LDA is INTEGER
221*> The leading dimension of A. It must be at least 2 (not 1!)
222*> and at least max( KK )+1.
223*> \endverbatim
224*>
225*> \param[out] SD
226*> \verbatim
227*> SD is DOUBLE PRECISION array, dimension (max(NN))
228*> Used to hold the diagonal of the tridiagonal matrix computed
229*> by ZHBTRD.
230*> \endverbatim
231*>
232*> \param[out] SE
233*> \verbatim
234*> SE is DOUBLE PRECISION array, dimension (max(NN))
235*> Used to hold the off-diagonal of the tridiagonal matrix
236*> computed by ZHBTRD.
237*> \endverbatim
238*>
239*> \param[out] D1
240*> \verbatim
241*> D1 is DOUBLE PRECISION array, dimension (max(NN))
242*> \endverbatim
243*>
244*> \param[out] D2
245*> \verbatim
246*> D2 is DOUBLE PRECISION array, dimension (max(NN))
247*> \endverbatim
248*>*> \param[out] D3
249*> \verbatim
250*> D3 is DOUBLE PRECISION array, dimension (max(NN))
251*> \endverbatim
252*>
253*> \param[out] U
254*> \verbatim
255*> U is COMPLEX*16 array, dimension (LDU, max(NN))
256*> Used to hold the unitary matrix computed by ZHBTRD.
257*> \endverbatim
258*>
259*> \param[in] LDU
260*> \verbatim
261*> LDU is INTEGER
262*> The leading dimension of U. It must be at least 1
263*> and at least max( NN ).
264*> \endverbatim
265*>
266*> \param[out] WORK
267*> \verbatim
268*> WORK is COMPLEX*16 array, dimension (LWORK)
269*> \endverbatim
270*>
271*> \param[in] LWORK
272*> \verbatim
273*> LWORK is INTEGER
274*> The number of entries in WORK. This must be at least
275*> max( LDA+1, max(NN)+1 )*max(NN).
276*> \endverbatim
277*>
278*> \param[out] RWORK
279*> \verbatim
280*> RWORK is DOUBLE PRECISION array
281*> \endverbatim
282*>
283*> \param[out] RESULT
284*> \verbatim
285*> RESULT is DOUBLE PRECISION array, dimension (4)
286*> The values computed by the tests described above.
287*> The values are currently limited to 1/ulp, to avoid
288*> overflow.
289*> \endverbatim
290*>
291*> \param[out] INFO
292*> \verbatim
293*> INFO is INTEGER
294*> If 0, then everything ran OK.
295*>
296*>-----------------------------------------------------------------------
297*>
298*> Some Local Variables and Parameters:
299*> ---- ----- --------- --- ----------
300*> ZERO, ONE Real 0 and 1.
301*> MAXTYP The number of types defined.
302*> NTEST The number of tests performed, or which can
303*> be performed so far, for the current matrix.
304*> NTESTT The total number of tests performed so far.
305*> NMAX Largest value in NN.
306*> NMATS The number of matrices generated so far.
307*> NERRS The number of tests which have exceeded THRESH
308*> so far.
309*> COND, IMODE Values to be passed to the matrix generators.
310*> ANORM Norm of A; passed to matrix generators.
311*>
312*> OVFL, UNFL Overflow and underflow thresholds.
313*> ULP, ULPINV Finest relative precision and its inverse.
314*> RTOVFL, RTUNFL Square roots of the previous 2 values.
315*> The following four arrays decode JTYPE:
316*> KTYPE(j) The general type (1-10) for type "j".
317*> KMODE(j) The MODE value to be passed to the matrix
318*> generator for type "j".
319*> KMAGN(j) The order of magnitude ( O(1),
320*> O(overflow^(1/2) ), O(underflow^(1/2) )
321*> \endverbatim
322*
323* Authors:
324* ========
325*
326*> \author Univ. of Tennessee
327*> \author Univ. of California Berkeley
328*> \author Univ. of Colorado Denver
329*> \author NAG Ltd.
330*
331*> \ingroup complex16_eig
332*
333* =====================================================================
334 SUBROUTINE zchkhb2stg( NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE,
335 $ ISEED, THRESH, NOUNIT, A, LDA, SD, SE, D1,
336 $ D2, D3, U, LDU, WORK, LWORK, RWORK, RESULT,
337 $ INFO )
338*
339* -- LAPACK test routine --
340* -- LAPACK is a software package provided by Univ. of Tennessee, --
341* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
342*
343* .. Scalar Arguments ..
344 INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES,
345 $ NWDTHS
346 DOUBLE PRECISION THRESH
347* ..
348* .. Array Arguments ..
349 LOGICAL DOTYPE( * )
350 INTEGER ISEED( 4 ), KK( * ), NN( * )
351 DOUBLE PRECISION RESULT( * ), RWORK( * ), SD( * ), SE( * ),
352 $ d1( * ), d2( * ), d3( * )
353 COMPLEX*16 A( LDA, * ), U( LDU, * ), WORK( * )
354* ..
355*
356* =====================================================================
357*
358* .. Parameters ..
359 COMPLEX*16 CZERO, CONE
360 PARAMETER ( CZERO = ( 0.0d+0, 0.0d+0 ),
361 $ cone = ( 1.0d+0, 0.0d+0 ) )
362 DOUBLE PRECISION ZERO, ONE, TWO, TEN
363 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
364 $ ten = 10.0d+0 )
365 DOUBLE PRECISION HALF
366 parameter( half = one / two )
367 INTEGER MAXTYP
368 parameter( maxtyp = 15 )
369* ..
370* .. Local Scalars ..
371 LOGICAL BADNN, BADNNB
372 INTEGER I, IINFO, IMODE, ITYPE, J, JC, JCOL, JR, JSIZE,
373 $ JTYPE, JWIDTH, K, KMAX, LH, LW, MTYPES, N,
374 $ nerrs, nmats, nmax, ntest, ntestt
375 DOUBLE PRECISION ANINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
376 $ TEMP1, TEMP2, TEMP3, TEMP4, ULP, ULPINV, UNFL
377* ..
378* .. Local Arrays ..
379 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KMAGN( MAXTYP ),
380 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
381* ..
382* .. External Functions ..
383 DOUBLE PRECISION DLAMCH
384 EXTERNAL DLAMCH
385* ..
386* .. External Subroutines ..
387 EXTERNAL dlasum, xerbla, zhbt21, zhbtrd, zlacpy, zlaset,
389* ..
390* .. Intrinsic Functions ..
391 INTRINSIC abs, dble, dconjg, max, min, sqrt
392* ..
393* .. Data statements ..
394 DATA ktype / 1, 2, 5*4, 5*5, 3*8 /
395 DATA kmagn / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
396 $ 2, 3 /
397 DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
398 $ 0, 0 /
399* ..
400* .. Executable Statements ..
401*
402* Check for errors
403*
404 ntestt = 0
405 info = 0
406*
407* Important constants
408*
409 badnn = .false.
410 nmax = 1
411 DO 10 j = 1, nsizes
412 nmax = max( nmax, nn( j ) )
413 IF( nn( j ).LT.0 )
414 $ badnn = .true.
415 10 CONTINUE
416*
417 badnnb = .false.
418 kmax = 0
419 DO 20 j = 1, nsizes
420 kmax = max( kmax, kk( j ) )
421 IF( kk( j ).LT.0 )
422 $ badnnb = .true.
423 20 CONTINUE
424 kmax = min( nmax-1, kmax )
425*
426* Check for errors
427*
428 IF( nsizes.LT.0 ) THEN
429 info = -1
430 ELSE IF( badnn ) THEN
431 info = -2
432 ELSE IF( nwdths.LT.0 ) THEN
433 info = -3
434 ELSE IF( badnnb ) THEN
435 info = -4
436 ELSE IF( ntypes.LT.0 ) THEN
437 info = -5
438 ELSE IF( lda.LT.kmax+1 ) THEN
439 info = -11
440 ELSE IF( ldu.LT.nmax ) THEN
441 info = -15
442 ELSE IF( ( max( lda, nmax )+1 )*nmax.GT.lwork ) THEN
443 info = -17
444 END IF
445*
446 IF( info.NE.0 ) THEN
447 CALL xerbla( 'ZCHKHB2STG', -info )
448 RETURN
449 END IF
450*
451* Quick return if possible
452*
453 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 .OR. nwdths.EQ.0 )
454 $ RETURN
455*
456* More Important constants
457*
458 unfl = dlamch( 'Safe minimum' )
459 ovfl = one / unfl
460 ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
461 ulpinv = one / ulp
462 rtunfl = sqrt( unfl )
463 rtovfl = sqrt( ovfl )
464*
465* Loop over sizes, types
466*
467 nerrs = 0
468 nmats = 0
469*
470 DO 190 jsize = 1, nsizes
471 n = nn( jsize )
472 aninv = one / dble( max( 1, n ) )
473*
474 DO 180 jwidth = 1, nwdths
475 k = kk( jwidth )
476 IF( k.GT.n )
477 $ GO TO 180
478 k = max( 0, min( n-1, k ) )
479*
480 IF( nsizes.NE.1 ) THEN
481 mtypes = min( maxtyp, ntypes )
482 ELSE
483 mtypes = min( maxtyp+1, ntypes )
484 END IF
485*
486 DO 170 jtype = 1, mtypes
487 IF( .NOT.dotype( jtype ) )
488 $ GO TO 170
489 nmats = nmats + 1
490 ntest = 0
491*
492 DO 30 j = 1, 4
493 ioldsd( j ) = iseed( j )
494 30 CONTINUE
495*
496* Compute "A".
497* Store as "Upper"; later, we will copy to other format.
498*
499* Control parameters:
500*
501* KMAGN KMODE KTYPE
502* =1 O(1) clustered 1 zero
503* =2 large clustered 2 identity
504* =3 small exponential (none)
505* =4 arithmetic diagonal, (w/ eigenvalues)
506* =5 random log hermitian, w/ eigenvalues
507* =6 random (none)
508* =7 random diagonal
509* =8 random hermitian
510* =9 positive definite
511* =10 diagonally dominant tridiagonal
512*
513 IF( mtypes.GT.maxtyp )
514 $ GO TO 100
515*
516 itype = ktype( jtype )
517 imode = kmode( jtype )
518*
519* Compute norm
520*
521 GO TO ( 40, 50, 60 )kmagn( jtype )
522*
523 40 CONTINUE
524 anorm = one
525 GO TO 70
526*
527 50 CONTINUE
528 anorm = ( rtovfl*ulp )*aninv
529 GO TO 70
530*
531 60 CONTINUE
532 anorm = rtunfl*n*ulpinv
533 GO TO 70
534*
535 70 CONTINUE
536*
537 CALL zlaset( 'Full', lda, n, czero, czero, a, lda )
538 iinfo = 0
539 IF( jtype.LE.15 ) THEN
540 cond = ulpinv
541 ELSE
542 cond = ulpinv*aninv / ten
543 END IF
544*
545* Special Matrices -- Identity & Jordan block
546*
547* Zero
548*
549 IF( itype.EQ.1 ) THEN
550 iinfo = 0
551*
552 ELSE IF( itype.EQ.2 ) THEN
553*
554* Identity
555*
556 DO 80 jcol = 1, n
557 a( k+1, jcol ) = anorm
558 80 CONTINUE
559*
560 ELSE IF( itype.EQ.4 ) THEN
561*
562* Diagonal Matrix, [Eigen]values Specified
563*
564 CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode,
565 $ cond, anorm, 0, 0, 'Q', a( k+1, 1 ), lda,
566 $ work, iinfo )
567*
568 ELSE IF( itype.EQ.5 ) THEN
569*
570* Hermitian, eigenvalues specified
571*
572 CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode,
573 $ cond, anorm, k, k, 'Q', a, lda, work,
574 $ iinfo )
575*
576 ELSE IF( itype.EQ.7 ) THEN
577*
578* Diagonal, random eigenvalues
579*
580 CALL zlatmr( n, n, 'S', iseed, 'H', work, 6, one,
581 $ cone, 'T', 'N', work( n+1 ), 1, one,
582 $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
583 $ zero, anorm, 'Q', a( k+1, 1 ), lda,
584 $ idumma, iinfo )
585*
586 ELSE IF( itype.EQ.8 ) THEN
587*
588* Hermitian, random eigenvalues
589*
590 CALL zlatmr( n, n, 'S', iseed, 'H', work, 6, one,
591 $ cone, 'T', 'N', work( n+1 ), 1, one,
592 $ work( 2*n+1 ), 1, one, 'N', idumma, k, k,
593 $ zero, anorm, 'Q', a, lda, idumma, iinfo )
594*
595 ELSE IF( itype.EQ.9 ) THEN
596*
597* Positive definite, eigenvalues specified.
598*
599 CALL zlatms( n, n, 'S', iseed, 'P', rwork, imode,
600 $ cond, anorm, k, k, 'Q', a, lda,
601 $ work( n+1 ), iinfo )
602*
603 ELSE IF( itype.EQ.10 ) THEN
604*
605* Positive definite tridiagonal, eigenvalues specified.
606*
607 IF( n.GT.1 )
608 $ k = max( 1, k )
609 CALL zlatms( n, n, 'S', iseed, 'P', rwork, imode,
610 $ cond, anorm, 1, 1, 'Q', a( k, 1 ), lda,
611 $ work, iinfo )
612 DO 90 i = 2, n
613 temp1 = abs( a( k, i ) ) /
614 $ sqrt( abs( a( k+1, i-1 )*a( k+1, i ) ) )
615 IF( temp1.GT.half ) THEN
616 a( k, i ) = half*sqrt( abs( a( k+1,
617 $ i-1 )*a( k+1, i ) ) )
618 END IF
619 90 CONTINUE
620*
621 ELSE
622*
623 iinfo = 1
624 END IF
625*
626 IF( iinfo.NE.0 ) THEN
627 WRITE( nounit, fmt = 9999 )'Generator', iinfo, n,
628 $ jtype, ioldsd
629 info = abs( iinfo )
630 RETURN
631 END IF
632*
633 100 CONTINUE
634*
635* Call ZHBTRD to compute S and U from upper triangle.
636*
637 CALL zlacpy( ' ', k+1, n, a, lda, work, lda )
638*
639 ntest = 1
640 CALL zhbtrd( 'V', 'U', n, k, work, lda, sd, se, u, ldu,
641 $ work( lda*n+1 ), iinfo )
642*
643 IF( iinfo.NE.0 ) THEN
644 WRITE( nounit, fmt = 9999 )'ZHBTRD(U)', iinfo, n,
645 $ jtype, ioldsd
646 info = abs( iinfo )
647 IF( iinfo.LT.0 ) THEN
648 RETURN
649 ELSE
650 result( 1 ) = ulpinv
651 GO TO 150
652 END IF
653 END IF
654*
655* Do tests 1 and 2
656*
657 CALL zhbt21( 'Upper', n, k, 1, a, lda, sd, se, u, ldu,
658 $ work, rwork, result( 1 ) )
659*
660* Before converting A into lower for DSBTRD, run DSYTRD_SB2ST
661* otherwise matrix A will be converted to lower and then need
662* to be converted back to upper in order to run the upper case
663* ofDSYTRD_SB2ST
664*
665* Compute D1 the eigenvalues resulting from the tridiagonal
666* form using the DSBTRD and used as reference to compare
667* with the DSYTRD_SB2ST routine
668*
669* Compute D1 from the DSBTRD and used as reference for the
670* DSYTRD_SB2ST
671*
672 CALL dcopy( n, sd, 1, d1, 1 )
673 IF( n.GT.0 )
674 $ CALL dcopy( n-1, se, 1, rwork, 1 )
675*
676 CALL zsteqr( 'N', n, d1, rwork, work, ldu,
677 $ rwork( n+1 ), iinfo )
678 IF( iinfo.NE.0 ) THEN
679 WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n,
680 $ jtype, ioldsd
681 info = abs( iinfo )
682 IF( iinfo.LT.0 ) THEN
683 RETURN
684 ELSE
685 result( 5 ) = ulpinv
686 GO TO 150
687 END IF
688 END IF
689*
690* DSYTRD_SB2ST Upper case is used to compute D2.
691* Note to set SD and SE to zero to be sure not reusing
692* the one from above. Compare it with D1 computed
693* using the DSBTRD.
694*
695 CALL dlaset( 'Full', n, 1, zero, zero, sd, n )
696 CALL dlaset( 'Full', n, 1, zero, zero, se, n )
697 CALL zlacpy( ' ', k+1, n, a, lda, u, ldu )
698 lh = max(1, 4*n)
699 lw = lwork - lh
700 CALL zhetrd_hb2st( 'N', 'N', "U", n, k, u, ldu, sd, se,
701 $ work, lh, work( lh+1 ), lw, iinfo )
702*
703* Compute D2 from the DSYTRD_SB2ST Upper case
704*
705 CALL dcopy( n, sd, 1, d2, 1 )
706 IF( n.GT.0 )
707 $ CALL dcopy( n-1, se, 1, rwork, 1 )
708*
709 CALL zsteqr( 'N', n, d2, rwork, work, ldu,
710 $ rwork( n+1 ), iinfo )
711 IF( iinfo.NE.0 ) THEN
712 WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n,
713 $ jtype, ioldsd
714 info = abs( iinfo )
715 IF( iinfo.LT.0 ) THEN
716 RETURN
717 ELSE
718 result( 5 ) = ulpinv
719 GO TO 150
720 END IF
721 END IF
722*
723* Convert A from Upper-Triangle-Only storage to
724* Lower-Triangle-Only storage.
725*
726 DO 120 jc = 1, n
727 DO 110 jr = 0, min( k, n-jc )
728 a( jr+1, jc ) = dconjg( a( k+1-jr, jc+jr ) )
729 110 CONTINUE
730 120 CONTINUE
731 DO 140 jc = n + 1 - k, n
732 DO 130 jr = min( k, n-jc ) + 1, k
733 a( jr+1, jc ) = zero
734 130 CONTINUE
735 140 CONTINUE
736*
737* Call ZHBTRD to compute S and U from lower triangle
738*
739 CALL zlacpy( ' ', k+1, n, a, lda, work, lda )
740*
741 ntest = 3
742 CALL zhbtrd( 'V', 'L', n, k, work, lda, sd, se, u, ldu,
743 $ work( lda*n+1 ), iinfo )
744*
745 IF( iinfo.NE.0 ) THEN
746 WRITE( nounit, fmt = 9999 )'ZHBTRD(L)', iinfo, n,
747 $ jtype, ioldsd
748 info = abs( iinfo )
749 IF( iinfo.LT.0 ) THEN
750 RETURN
751 ELSE
752 result( 3 ) = ulpinv
753 GO TO 150
754 END IF
755 END IF
756 ntest = 4
757*
758* Do tests 3 and 4
759*
760 CALL zhbt21( 'Lower', n, k, 1, a, lda, sd, se, u, ldu,
761 $ work, rwork, result( 3 ) )
762*
763* DSYTRD_SB2ST Lower case is used to compute D3.
764* Note to set SD and SE to zero to be sure not reusing
765* the one from above. Compare it with D1 computed
766* using the DSBTRD.
767*
768 CALL dlaset( 'Full', n, 1, zero, zero, sd, n )
769 CALL dlaset( 'Full', n, 1, zero, zero, se, n )
770 CALL zlacpy( ' ', k+1, n, a, lda, u, ldu )
771 lh = max(1, 4*n)
772 lw = lwork - lh
773 CALL zhetrd_hb2st( 'N', 'N', "L", n, k, u, ldu, sd, se,
774 $ work, lh, work( lh+1 ), lw, iinfo )
775*
776* Compute D3 from the 2-stage Upper case
777*
778 CALL dcopy( n, sd, 1, d3, 1 )
779 IF( n.GT.0 )
780 $ CALL dcopy( n-1, se, 1, rwork, 1 )
781*
782 CALL zsteqr( 'N', n, d3, rwork, work, ldu,
783 $ rwork( n+1 ), iinfo )
784 IF( iinfo.NE.0 ) THEN
785 WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n,
786 $ jtype, ioldsd
787 info = abs( iinfo )
788 IF( iinfo.LT.0 ) THEN
789 RETURN
790 ELSE
791 result( 6 ) = ulpinv
792 GO TO 150
793 END IF
794 END IF
795*
796*
797* Do Tests 3 and 4 which are similar to 11 and 12 but with the
798* D1 computed using the standard 1-stage reduction as reference
799*
800 ntest = 6
801 temp1 = zero
802 temp2 = zero
803 temp3 = zero
804 temp4 = zero
805*
806 DO 151 j = 1, n
807 temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
808 temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
809 temp3 = max( temp3, abs( d1( j ) ), abs( d3( j ) ) )
810 temp4 = max( temp4, abs( d1( j )-d3( j ) ) )
811 151 CONTINUE
812*
813 result(5) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
814 result(6) = temp4 / max( unfl, ulp*max( temp3, temp4 ) )
815*
816* End of Loop -- Check for RESULT(j) > THRESH
817*
818 150 CONTINUE
819 ntestt = ntestt + ntest
820*
821* Print out tests which fail.
822*
823 DO 160 jr = 1, ntest
824 IF( result( jr ).GE.thresh ) THEN
825*
826* If this is the first test to fail,
827* print a header to the data file.
828*
829 IF( nerrs.EQ.0 ) THEN
830 WRITE( nounit, fmt = 9998 )'ZHB'
831 WRITE( nounit, fmt = 9997 )
832 WRITE( nounit, fmt = 9996 )
833 WRITE( nounit, fmt = 9995 )'Hermitian'
834 WRITE( nounit, fmt = 9994 )'unitary', '*',
835 $ 'conjugate transpose', ( '*', j = 1, 6 )
836 END IF
837 nerrs = nerrs + 1
838 WRITE( nounit, fmt = 9993 )n, k, ioldsd, jtype,
839 $ jr, result( jr )
840 END IF
841 160 CONTINUE
842*
843 170 CONTINUE
844 180 CONTINUE
845 190 CONTINUE
846*
847* Summary
848*
849 CALL dlasum( 'ZHB', nounit, nerrs, ntestt )
850 RETURN
851*
852 9999 FORMAT( ' ZCHKHB2STG: ', a, ' returned INFO=', i6, '.', / 9x,
853 $ 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5,
854 $ ')' )
855 9998 FORMAT( / 1x, a3,
856 $ ' -- Complex Hermitian Banded Tridiagonal Reduction Routines'
857 $ )
858 9997 FORMAT( ' Matrix types (see DCHK23 for details): ' )
859*
860 9996 FORMAT( / ' Special Matrices:',
861 $ / ' 1=Zero matrix. ',
862 $ ' 5=Diagonal: clustered entries.',
863 $ / ' 2=Identity matrix. ',
864 $ ' 6=Diagonal: large, evenly spaced.',
865 $ / ' 3=Diagonal: evenly spaced entries. ',
866 $ ' 7=Diagonal: small, evenly spaced.',
867 $ / ' 4=Diagonal: geometr. spaced entries.' )
868 9995 FORMAT( ' Dense ', a, ' Banded Matrices:',
869 $ / ' 8=Evenly spaced eigenvals. ',
870 $ ' 12=Small, evenly spaced eigenvals.',
871 $ / ' 9=Geometrically spaced eigenvals. ',
872 $ ' 13=Matrix with random O(1) entries.',
873 $ / ' 10=Clustered eigenvalues. ',
874 $ ' 14=Matrix with large random entries.',
875 $ / ' 11=Large, evenly spaced eigenvals. ',
876 $ ' 15=Matrix with small random entries.' )
877*
878 9994 FORMAT( / ' Tests performed: (S is Tridiag, U is ', a, ',',
879 $ / 20x, a, ' means ', a, '.', / ' UPLO=''U'':',
880 $ / ' 1= | A - U S U', a1, ' | / ( |A| n ulp ) ',
881 $ ' 2= | I - U U', a1, ' | / ( n ulp )', / ' UPLO=''L'':',
882 $ / ' 3= | A - U S U', a1, ' | / ( |A| n ulp ) ',
883 $ ' 4= | I - U U', a1, ' | / ( n ulp )' / ' Eig check:',
884 $ /' 5= | D1 - D2', '', ' | / ( |D1| ulp ) ',
885 $ ' 6= | D1 - D3', '', ' | / ( |D1| ulp ) ' )
886 9993 FORMAT( ' N=', i5, ', K=', i4, ', seed=', 4( i4, ',' ), ' type ',
887 $ i2, ', test(', i2, ')=', g10.3 )
888*
889* End of ZCHKHB2STG
890*
891 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dlasum(type, iounit, ie, nrun)
DLASUM
Definition dlasum.f:43
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine zhbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
ZHBTRD
Definition zhbtrd.f:163
subroutine zhetrd_hb2st(stage1, vect, uplo, n, kd, ab, ldab, d, e, hous, lhous, work, lwork, info)
ZHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric tridiagonal form T
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:132
subroutine zchkhb2stg(nsizes, nn, nwdths, kk, ntypes, dotype, iseed, thresh, nounit, a, lda, sd, se, d1, d2, d3, u, ldu, work, lwork, rwork, result, info)
ZCHKHB2STG
Definition zchkhb2stg.f:338
subroutine zhbt21(uplo, n, ka, ks, a, lda, d, e, u, ldu, work, rwork, result)
ZHBT21
Definition zhbt21.f:152
subroutine zlatmr(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)
ZLATMR
Definition zlatmr.f:490
subroutine zlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
ZLATMS
Definition zlatms.f:332