LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zchkhb2stg()

subroutine zchkhb2stg ( integer nsizes,
integer, dimension( * ) nn,
integer nwdths,
integer, dimension( * ) kk,
integer ntypes,
logical, dimension( * ) dotype,
integer, dimension( 4 ) iseed,
double precision thresh,
integer nounit,
complex*16, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) sd,
double precision, dimension( * ) se,
double precision, dimension( * ) d1,
double precision, dimension( * ) d2,
double precision, dimension( * ) d3,
complex*16, dimension( ldu, * ) u,
integer ldu,
complex*16, dimension( * ) work,
integer lwork,
double precision, dimension( * ) rwork,
double precision, dimension( * ) result,
integer info )

ZCHKHB2STG

Purpose:
!>
!> ZCHKHB2STG tests the reduction of a Hermitian band matrix to tridiagonal
!> from, used with the Hermitian eigenvalue problem.
!>
!> ZHBTRD factors a Hermitian band matrix A as  U S U* , where * means
!> conjugate transpose, S is symmetric tridiagonal, and U is unitary.
!> ZHBTRD can use either just the lower or just the upper triangle
!> of A; ZCHKHB2STG checks both cases.
!>
!> ZHETRD_HB2ST factors a Hermitian band matrix A as  U S U* ,
!> where * means conjugate transpose, S is symmetric tridiagonal, and U is
!> unitary. ZHETRD_HB2ST can use either just the lower or just
!> the upper triangle of A; ZCHKHB2STG checks both cases.
!>
!> DSTEQR factors S as  Z D1 Z'.
!> D1 is the matrix of eigenvalues computed when Z is not computed
!> and from the S resulting of DSBTRD  (used as reference for DSYTRD_SB2ST)
!> D2 is the matrix of eigenvalues computed when Z is not computed
!> and from the S resulting of DSYTRD_SB2ST .
!> D3 is the matrix of eigenvalues computed when Z is not computed
!> and from the S resulting of DSYTRD_SB2ST .
!>
!> When ZCHKHB2STG is called, a number of matrix  (), a number
!> of bandwidths (), and a number of matrix  are
!> specified.  For each size (), each bandwidth () less than or
!> equal to , and each type of matrix, one matrix will be generated
!> and used to test the hermitian banded reduction routine.  For each
!> matrix, a number of tests will be performed:
!>
!> (1)     | A - V S V* | / ( |A| n ulp )  computed by ZHBTRD with
!>                                         UPLO='U'
!>
!> (2)     | I - UU* | / ( n ulp )
!>
!> (3)     | A - V S V* | / ( |A| n ulp )  computed by ZHBTRD with
!>                                         UPLO='L'
!>
!> (4)     | I - UU* | / ( n ulp )
!>
!> (5)     | D1 - D2 | / ( |D1| ulp )      where D1 is computed by
!>                                         DSBTRD with UPLO='U' and
!>                                         D2 is computed by
!>                                         ZHETRD_HB2ST with UPLO='U'
!>
!> (6)     | D1 - D3 | / ( |D1| ulp )      where D1 is computed by
!>                                         DSBTRD with UPLO='U' and
!>                                         D3 is computed by
!>                                         ZHETRD_HB2ST with UPLO='L'
!>
!> The  are specified by an array NN(1:NSIZES); the value of
!> each element NN(j) specifies one size.
!> The  are specified by a logical array DOTYPE( 1:NTYPES );
!> if DOTYPE(j) is .TRUE., then matrix type  will be generated.
!> Currently, the list of possible types is:
!>
!> (1)  The zero matrix.
!> (2)  The identity matrix.
!>
!> (3)  A diagonal matrix with evenly spaced entries
!>      1, ..., ULP  and random signs.
!>      (ULP = (first number larger than 1) - 1 )
!> (4)  A diagonal matrix with geometrically spaced entries
!>      1, ..., ULP  and random signs.
!> (5)  A diagonal matrix with  entries 1, ULP, ..., ULP
!>      and random signs.
!>
!> (6)  Same as (4), but multiplied by SQRT( overflow threshold )
!> (7)  Same as (4), but multiplied by SQRT( underflow threshold )
!>
!> (8)  A matrix of the form  U* D U, where U is unitary and
!>      D has evenly spaced entries 1, ..., ULP with random signs
!>      on the diagonal.
!>
!> (9)  A matrix of the form  U* D U, where U is unitary and
!>      D has geometrically spaced entries 1, ..., ULP with random
!>      signs on the diagonal.
!>
!> (10) A matrix of the form  U* D U, where U is unitary and
!>      D has  entries 1, ULP,..., ULP with random
!>      signs on the diagonal.
!>
!> (11) Same as (8), but multiplied by SQRT( overflow threshold )
!> (12) Same as (8), but multiplied by SQRT( underflow threshold )
!>
!> (13) Hermitian matrix with random entries chosen from (-1,1).
!> (14) Same as (13), but multiplied by SQRT( overflow threshold )
!> (15) Same as (13), but multiplied by SQRT( underflow threshold )
!> 
Parameters
[in]NSIZES
!>          NSIZES is INTEGER
!>          The number of sizes of matrices to use.  If it is zero,
!>          ZCHKHB2STG does nothing.  It must be at least zero.
!> 
[in]NN
!>          NN is INTEGER array, dimension (NSIZES)
!>          An array containing the sizes to be used for the matrices.
!>          Zero values will be skipped.  The values must be at least
!>          zero.
!> 
[in]NWDTHS
!>          NWDTHS is INTEGER
!>          The number of bandwidths to use.  If it is zero,
!>          ZCHKHB2STG does nothing.  It must be at least zero.
!> 
[in]KK
!>          KK is INTEGER array, dimension (NWDTHS)
!>          An array containing the bandwidths to be used for the band
!>          matrices.  The values must be at least zero.
!> 
[in]NTYPES
!>          NTYPES is INTEGER
!>          The number of elements in DOTYPE.   If it is zero, ZCHKHB2STG
!>          does nothing.  It must be at least zero.  If it is MAXTYP+1
!>          and NSIZES is 1, then an additional type, MAXTYP+1 is
!>          defined, which is to use whatever matrix is in A.  This
!>          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
!>          DOTYPE(MAXTYP+1) is .TRUE. .
!> 
[in]DOTYPE
!>          DOTYPE is LOGICAL array, dimension (NTYPES)
!>          If DOTYPE(j) is .TRUE., then for each size in NN a
!>          matrix of that size and of type j will be generated.
!>          If NTYPES is smaller than the maximum number of types
!>          defined (PARAMETER MAXTYP), then types NTYPES+1 through
!>          MAXTYP will not be generated.  If NTYPES is larger
!>          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
!>          will be ignored.
!> 
[in,out]ISEED
!>          ISEED is INTEGER array, dimension (4)
!>          On entry ISEED specifies the seed of the random number
!>          generator. The array elements should be between 0 and 4095;
!>          if not they will be reduced mod 4096.  Also, ISEED(4) must
!>          be odd.  The random number generator uses a linear
!>          congruential sequence limited to small integers, and so
!>          should produce machine independent random numbers. The
!>          values of ISEED are changed on exit, and can be used in the
!>          next call to ZCHKHB2STG to continue the same random number
!>          sequence.
!> 
[in]THRESH
!>          THRESH is DOUBLE PRECISION
!>          A test will count as  if the , computed as
!>          described above, exceeds THRESH.  Note that the error
!>          is scaled to be O(1), so THRESH should be a reasonably
!>          small multiple of 1, e.g., 10 or 100.  In particular,
!>          it should not depend on the precision (single vs. double)
!>          or the size of the matrix.  It must be at least zero.
!> 
[in]NOUNIT
!>          NOUNIT is INTEGER
!>          The FORTRAN unit number for printing out error messages
!>          (e.g., if a routine returns IINFO not equal to 0.)
!> 
[in,out]A
!>          A is COMPLEX*16 array, dimension
!>                            (LDA, max(NN))
!>          Used to hold the matrix whose eigenvalues are to be
!>          computed.
!> 
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of A.  It must be at least 2 (not 1!)
!>          and at least max( KK )+1.
!> 
[out]SD
!>          SD is DOUBLE PRECISION array, dimension (max(NN))
!>          Used to hold the diagonal of the tridiagonal matrix computed
!>          by ZHBTRD.
!> 
[out]SE
!>          SE is DOUBLE PRECISION array, dimension (max(NN))
!>          Used to hold the off-diagonal of the tridiagonal matrix
!>          computed by ZHBTRD.
!> 
[out]D1
!>          D1 is DOUBLE PRECISION array, dimension (max(NN))
!> 
[out]D2
!>          D2 is DOUBLE PRECISION array, dimension (max(NN))
!> 
*>
[out]D3
!>          D3 is DOUBLE PRECISION array, dimension (max(NN))
!> 
[out]U
!>          U is COMPLEX*16 array, dimension (LDU, max(NN))
!>          Used to hold the unitary matrix computed by ZHBTRD.
!> 
[in]LDU
!>          LDU is INTEGER
!>          The leading dimension of U.  It must be at least 1
!>          and at least max( NN ).
!> 
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (LWORK)
!> 
[in]LWORK
!>          LWORK is INTEGER
!>          The number of entries in WORK.  This must be at least
!>          max( LDA+1, max(NN)+1 )*max(NN).
!> 
[out]RWORK
!>          RWORK is DOUBLE PRECISION array
!> 
[out]RESULT
!>          RESULT is DOUBLE PRECISION array, dimension (4)
!>          The values computed by the tests described above.
!>          The values are currently limited to 1/ulp, to avoid
!>          overflow.
!> 
[out]INFO
!>          INFO is INTEGER
!>          If 0, then everything ran OK.
!>
!>-----------------------------------------------------------------------
!>
!>       Some Local Variables and Parameters:
!>       ---- ----- --------- --- ----------
!>       ZERO, ONE       Real 0 and 1.
!>       MAXTYP          The number of types defined.
!>       NTEST           The number of tests performed, or which can
!>                       be performed so far, for the current matrix.
!>       NTESTT          The total number of tests performed so far.
!>       NMAX            Largest value in NN.
!>       NMATS           The number of matrices generated so far.
!>       NERRS           The number of tests which have exceeded THRESH
!>                       so far.
!>       COND, IMODE     Values to be passed to the matrix generators.
!>       ANORM           Norm of A; passed to matrix generators.
!>
!>       OVFL, UNFL      Overflow and underflow thresholds.
!>       ULP, ULPINV     Finest relative precision and its inverse.
!>       RTOVFL, RTUNFL  Square roots of the previous 2 values.
!>               The following four arrays decode JTYPE:
!>       KTYPE(j)        The general type (1-10) for type .
!>       KMODE(j)        The MODE value to be passed to the matrix
!>                       generator for type .
!>       KMAGN(j)        The order of magnitude ( O(1),
!>                       O(overflow^(1/2) ), O(underflow^(1/2) )
!> 
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 334 of file zchkhb2stg.f.

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*
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:161
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:101
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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:108
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:104
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:130
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
Here is the call graph for this function:
Here is the caller graph for this function: