LAPACK 3.12.0
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 "U" (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 "U".
 D3 is the matrix of eigenvalues computed when Z is not computed
 and from the S resulting of DSYTRD_SB2ST "L".

 When ZCHKHB2STG is called, a number of matrix "sizes" ("n's"), a number
 of bandwidths ("k's"), and a number of matrix "types" are
 specified.  For each size ("n"), each bandwidth ("k") less than or
 equal to "n", 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 "sizes" are specified by an array NN(1:NSIZES); the value of
 each element NN(j) specifies one size.
 The "types" are specified by a logical array DOTYPE( 1:NTYPES );
 if DOTYPE(j) is .TRUE., then matrix type "j" 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 "clustered" 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 "clustered" 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 "failed" if the "error", 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 "j".
       KMODE(j)        The MODE value to be passed to the matrix
                       generator for type "j".
       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: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
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: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 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: