LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ 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 
)

ZCHKHBSTG

Purpose:
 ZCHKHBSTG 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; ZCHKHBSTG 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; ZCHKHBSTG 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 ZCHKHBSTG 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,
          ZCHKHBSTG 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,
          ZCHKHBSTG 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, ZCHKHBSTG
          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 ZCHKHBSTG 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]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.
Date
June 2017

Definition at line 325 of file zchkhb2stg.f.

325 *
326 * -- LAPACK test routine (version 3.7.1) --
327 * -- LAPACK is a software package provided by Univ. of Tennessee, --
328 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
329 * June 2017
330 *
331 * .. Scalar Arguments ..
332  INTEGER info, lda, ldu, lwork, nounit, nsizes, ntypes,
333  $ nwdths
334  DOUBLE PRECISION thresh
335 * ..
336 * .. Array Arguments ..
337  LOGICAL dotype( * )
338  INTEGER iseed( 4 ), kk( * ), nn( * )
339  DOUBLE PRECISION result( * ), rwork( * ), sd( * ), se( * ),
340  $ d1( * ), d2( * ), d3( * )
341  COMPLEX*16 a( lda, * ), u( ldu, * ), work( * )
342 * ..
343 *
344 * =====================================================================
345 *
346 * .. Parameters ..
347  COMPLEX*16 czero, cone
348  parameter( czero = ( 0.0d+0, 0.0d+0 ),
349  $ cone = ( 1.0d+0, 0.0d+0 ) )
350  DOUBLE PRECISION zero, one, two, ten
351  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
352  $ ten = 10.0d+0 )
353  DOUBLE PRECISION half
354  parameter( half = one / two )
355  INTEGER maxtyp
356  parameter( maxtyp = 15 )
357 * ..
358 * .. Local Scalars ..
359  LOGICAL badnn, badnnb
360  INTEGER i, iinfo, imode, itype, j, jc, jcol, jr, jsize,
361  $ jtype, jwidth, k, kmax, lh, lw, mtypes, n,
362  $ nerrs, nmats, nmax, ntest, ntestt
363  DOUBLE PRECISION aninv, anorm, cond, ovfl, rtovfl, rtunfl,
364  $ temp1, temp2, temp3, temp4, ulp, ulpinv, unfl
365 * ..
366 * .. Local Arrays ..
367  INTEGER idumma( 1 ), ioldsd( 4 ), kmagn( maxtyp ),
368  $ kmode( maxtyp ), ktype( maxtyp )
369 * ..
370 * .. External Functions ..
371  DOUBLE PRECISION dlamch
372  EXTERNAL dlamch
373 * ..
374 * .. External Subroutines ..
375  EXTERNAL dlasum, xerbla, zhbt21, zhbtrd, zlacpy, zlaset,
376  $ zlatmr, zlatms, zhetrd_hb2st, zsteqr
377 * ..
378 * .. Intrinsic Functions ..
379  INTRINSIC abs, dble, dconjg, max, min, sqrt
380 * ..
381 * .. Data statements ..
382  DATA ktype / 1, 2, 5*4, 5*5, 3*8 /
383  DATA kmagn / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
384  $ 2, 3 /
385  DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
386  $ 0, 0 /
387 * ..
388 * .. Executable Statements ..
389 *
390 * Check for errors
391 *
392  ntestt = 0
393  info = 0
394 *
395 * Important constants
396 *
397  badnn = .false.
398  nmax = 1
399  DO 10 j = 1, nsizes
400  nmax = max( nmax, nn( j ) )
401  IF( nn( j ).LT.0 )
402  $ badnn = .true.
403  10 CONTINUE
404 *
405  badnnb = .false.
406  kmax = 0
407  DO 20 j = 1, nsizes
408  kmax = max( kmax, kk( j ) )
409  IF( kk( j ).LT.0 )
410  $ badnnb = .true.
411  20 CONTINUE
412  kmax = min( nmax-1, kmax )
413 *
414 * Check for errors
415 *
416  IF( nsizes.LT.0 ) THEN
417  info = -1
418  ELSE IF( badnn ) THEN
419  info = -2
420  ELSE IF( nwdths.LT.0 ) THEN
421  info = -3
422  ELSE IF( badnnb ) THEN
423  info = -4
424  ELSE IF( ntypes.LT.0 ) THEN
425  info = -5
426  ELSE IF( lda.LT.kmax+1 ) THEN
427  info = -11
428  ELSE IF( ldu.LT.nmax ) THEN
429  info = -15
430  ELSE IF( ( max( lda, nmax )+1 )*nmax.GT.lwork ) THEN
431  info = -17
432  END IF
433 *
434  IF( info.NE.0 ) THEN
435  CALL xerbla( 'ZCHKHBSTG', -info )
436  RETURN
437  END IF
438 *
439 * Quick return if possible
440 *
441  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 .OR. nwdths.EQ.0 )
442  $ RETURN
443 *
444 * More Important constants
445 *
446  unfl = dlamch( 'Safe minimum' )
447  ovfl = one / unfl
448  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
449  ulpinv = one / ulp
450  rtunfl = sqrt( unfl )
451  rtovfl = sqrt( ovfl )
452 *
453 * Loop over sizes, types
454 *
455  nerrs = 0
456  nmats = 0
457 *
458  DO 190 jsize = 1, nsizes
459  n = nn( jsize )
460  aninv = one / dble( max( 1, n ) )
461 *
462  DO 180 jwidth = 1, nwdths
463  k = kk( jwidth )
464  IF( k.GT.n )
465  $ GO TO 180
466  k = max( 0, min( n-1, k ) )
467 *
468  IF( nsizes.NE.1 ) THEN
469  mtypes = min( maxtyp, ntypes )
470  ELSE
471  mtypes = min( maxtyp+1, ntypes )
472  END IF
473 *
474  DO 170 jtype = 1, mtypes
475  IF( .NOT.dotype( jtype ) )
476  $ GO TO 170
477  nmats = nmats + 1
478  ntest = 0
479 *
480  DO 30 j = 1, 4
481  ioldsd( j ) = iseed( j )
482  30 CONTINUE
483 *
484 * Compute "A".
485 * Store as "Upper"; later, we will copy to other format.
486 *
487 * Control parameters:
488 *
489 * KMAGN KMODE KTYPE
490 * =1 O(1) clustered 1 zero
491 * =2 large clustered 2 identity
492 * =3 small exponential (none)
493 * =4 arithmetic diagonal, (w/ eigenvalues)
494 * =5 random log hermitian, w/ eigenvalues
495 * =6 random (none)
496 * =7 random diagonal
497 * =8 random hermitian
498 * =9 positive definite
499 * =10 diagonally dominant tridiagonal
500 *
501  IF( mtypes.GT.maxtyp )
502  $ GO TO 100
503 *
504  itype = ktype( jtype )
505  imode = kmode( jtype )
506 *
507 * Compute norm
508 *
509  GO TO ( 40, 50, 60 )kmagn( jtype )
510 *
511  40 CONTINUE
512  anorm = one
513  GO TO 70
514 *
515  50 CONTINUE
516  anorm = ( rtovfl*ulp )*aninv
517  GO TO 70
518 *
519  60 CONTINUE
520  anorm = rtunfl*n*ulpinv
521  GO TO 70
522 *
523  70 CONTINUE
524 *
525  CALL zlaset( 'Full', lda, n, czero, czero, a, lda )
526  iinfo = 0
527  IF( jtype.LE.15 ) THEN
528  cond = ulpinv
529  ELSE
530  cond = ulpinv*aninv / ten
531  END IF
532 *
533 * Special Matrices -- Identity & Jordan block
534 *
535 * Zero
536 *
537  IF( itype.EQ.1 ) THEN
538  iinfo = 0
539 *
540  ELSE IF( itype.EQ.2 ) THEN
541 *
542 * Identity
543 *
544  DO 80 jcol = 1, n
545  a( k+1, jcol ) = anorm
546  80 CONTINUE
547 *
548  ELSE IF( itype.EQ.4 ) THEN
549 *
550 * Diagonal Matrix, [Eigen]values Specified
551 *
552  CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode,
553  $ cond, anorm, 0, 0, 'Q', a( k+1, 1 ), lda,
554  $ work, iinfo )
555 *
556  ELSE IF( itype.EQ.5 ) THEN
557 *
558 * Hermitian, eigenvalues specified
559 *
560  CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode,
561  $ cond, anorm, k, k, 'Q', a, lda, work,
562  $ iinfo )
563 *
564  ELSE IF( itype.EQ.7 ) THEN
565 *
566 * Diagonal, random eigenvalues
567 *
568  CALL zlatmr( n, n, 'S', iseed, 'H', work, 6, one,
569  $ cone, 'T', 'N', work( n+1 ), 1, one,
570  $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
571  $ zero, anorm, 'Q', a( k+1, 1 ), lda,
572  $ idumma, iinfo )
573 *
574  ELSE IF( itype.EQ.8 ) THEN
575 *
576 * Hermitian, random eigenvalues
577 *
578  CALL zlatmr( n, n, 'S', iseed, 'H', work, 6, one,
579  $ cone, 'T', 'N', work( n+1 ), 1, one,
580  $ work( 2*n+1 ), 1, one, 'N', idumma, k, k,
581  $ zero, anorm, 'Q', a, lda, idumma, iinfo )
582 *
583  ELSE IF( itype.EQ.9 ) THEN
584 *
585 * Positive definite, eigenvalues specified.
586 *
587  CALL zlatms( n, n, 'S', iseed, 'P', rwork, imode,
588  $ cond, anorm, k, k, 'Q', a, lda,
589  $ work( n+1 ), iinfo )
590 *
591  ELSE IF( itype.EQ.10 ) THEN
592 *
593 * Positive definite tridiagonal, eigenvalues specified.
594 *
595  IF( n.GT.1 )
596  $ k = max( 1, k )
597  CALL zlatms( n, n, 'S', iseed, 'P', rwork, imode,
598  $ cond, anorm, 1, 1, 'Q', a( k, 1 ), lda,
599  $ work, iinfo )
600  DO 90 i = 2, n
601  temp1 = abs( a( k, i ) ) /
602  $ sqrt( abs( a( k+1, i-1 )*a( k+1, i ) ) )
603  IF( temp1.GT.half ) THEN
604  a( k, i ) = half*sqrt( abs( a( k+1,
605  $ i-1 )*a( k+1, i ) ) )
606  END IF
607  90 CONTINUE
608 *
609  ELSE
610 *
611  iinfo = 1
612  END IF
613 *
614  IF( iinfo.NE.0 ) THEN
615  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n,
616  $ jtype, ioldsd
617  info = abs( iinfo )
618  RETURN
619  END IF
620 *
621  100 CONTINUE
622 *
623 * Call ZHBTRD to compute S and U from upper triangle.
624 *
625  CALL zlacpy( ' ', k+1, n, a, lda, work, lda )
626 *
627  ntest = 1
628  CALL zhbtrd( 'V', 'U', n, k, work, lda, sd, se, u, ldu,
629  $ work( lda*n+1 ), iinfo )
630 *
631  IF( iinfo.NE.0 ) THEN
632  WRITE( nounit, fmt = 9999 )'ZHBTRD(U)', iinfo, n,
633  $ jtype, ioldsd
634  info = abs( iinfo )
635  IF( iinfo.LT.0 ) THEN
636  RETURN
637  ELSE
638  result( 1 ) = ulpinv
639  GO TO 150
640  END IF
641  END IF
642 *
643 * Do tests 1 and 2
644 *
645  CALL zhbt21( 'Upper', n, k, 1, a, lda, sd, se, u, ldu,
646  $ work, rwork, result( 1 ) )
647 *
648 * Before converting A into lower for DSBTRD, run DSYTRD_SB2ST
649 * otherwise matrix A will be converted to lower and then need
650 * to be converted back to upper in order to run the upper case
651 * ofDSYTRD_SB2ST
652 *
653 * Compute D1 the eigenvalues resulting from the tridiagonal
654 * form using the DSBTRD and used as reference to compare
655 * with the DSYTRD_SB2ST routine
656 *
657 * Compute D1 from the DSBTRD and used as reference for the
658 * DSYTRD_SB2ST
659 *
660  CALL dcopy( n, sd, 1, d1, 1 )
661  IF( n.GT.0 )
662  $ CALL dcopy( n-1, se, 1, rwork, 1 )
663 *
664  CALL zsteqr( 'N', n, d1, rwork, work, ldu,
665  $ rwork( n+1 ), iinfo )
666  IF( iinfo.NE.0 ) THEN
667  WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n,
668  $ jtype, ioldsd
669  info = abs( iinfo )
670  IF( iinfo.LT.0 ) THEN
671  RETURN
672  ELSE
673  result( 5 ) = ulpinv
674  GO TO 150
675  END IF
676  END IF
677 *
678 * DSYTRD_SB2ST Upper case is used to compute D2.
679 * Note to set SD and SE to zero to be sure not reusing
680 * the one from above. Compare it with D1 computed
681 * using the DSBTRD.
682 *
683  CALL dlaset( 'Full', n, 1, zero, zero, sd, 1 )
684  CALL dlaset( 'Full', n, 1, zero, zero, se, 1 )
685  CALL zlacpy( ' ', k+1, n, a, lda, u, ldu )
686  lh = max(1, 4*n)
687  lw = lwork - lh
688  CALL zhetrd_hb2st( 'N', 'N', "U", n, k, u, ldu, sd, se,
689  $ work, lh, work( lh+1 ), lw, iinfo )
690 *
691 * Compute D2 from the DSYTRD_SB2ST Upper case
692 *
693  CALL dcopy( n, sd, 1, d2, 1 )
694  IF( n.GT.0 )
695  $ CALL dcopy( n-1, se, 1, rwork, 1 )
696 *
697  CALL zsteqr( 'N', n, d2, rwork, work, ldu,
698  $ rwork( n+1 ), iinfo )
699  IF( iinfo.NE.0 ) THEN
700  WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n,
701  $ jtype, ioldsd
702  info = abs( iinfo )
703  IF( iinfo.LT.0 ) THEN
704  RETURN
705  ELSE
706  result( 5 ) = ulpinv
707  GO TO 150
708  END IF
709  END IF
710 *
711 * Convert A from Upper-Triangle-Only storage to
712 * Lower-Triangle-Only storage.
713 *
714  DO 120 jc = 1, n
715  DO 110 jr = 0, min( k, n-jc )
716  a( jr+1, jc ) = dconjg( a( k+1-jr, jc+jr ) )
717  110 CONTINUE
718  120 CONTINUE
719  DO 140 jc = n + 1 - k, n
720  DO 130 jr = min( k, n-jc ) + 1, k
721  a( jr+1, jc ) = zero
722  130 CONTINUE
723  140 CONTINUE
724 *
725 * Call ZHBTRD to compute S and U from lower triangle
726 *
727  CALL zlacpy( ' ', k+1, n, a, lda, work, lda )
728 *
729  ntest = 3
730  CALL zhbtrd( 'V', 'L', n, k, work, lda, sd, se, u, ldu,
731  $ work( lda*n+1 ), iinfo )
732 *
733  IF( iinfo.NE.0 ) THEN
734  WRITE( nounit, fmt = 9999 )'ZHBTRD(L)', iinfo, n,
735  $ jtype, ioldsd
736  info = abs( iinfo )
737  IF( iinfo.LT.0 ) THEN
738  RETURN
739  ELSE
740  result( 3 ) = ulpinv
741  GO TO 150
742  END IF
743  END IF
744  ntest = 4
745 *
746 * Do tests 3 and 4
747 *
748  CALL zhbt21( 'Lower', n, k, 1, a, lda, sd, se, u, ldu,
749  $ work, rwork, result( 3 ) )
750 *
751 * DSYTRD_SB2ST Lower case is used to compute D3.
752 * Note to set SD and SE to zero to be sure not reusing
753 * the one from above. Compare it with D1 computed
754 * using the DSBTRD.
755 *
756  CALL dlaset( 'Full', n, 1, zero, zero, sd, 1 )
757  CALL dlaset( 'Full', n, 1, zero, zero, se, 1 )
758  CALL zlacpy( ' ', k+1, n, a, lda, u, ldu )
759  lh = max(1, 4*n)
760  lw = lwork - lh
761  CALL zhetrd_hb2st( 'N', 'N', "L", n, k, u, ldu, sd, se,
762  $ work, lh, work( lh+1 ), lw, iinfo )
763 *
764 * Compute D3 from the 2-stage Upper case
765 *
766  CALL dcopy( n, sd, 1, d3, 1 )
767  IF( n.GT.0 )
768  $ CALL dcopy( n-1, se, 1, rwork, 1 )
769 *
770  CALL zsteqr( 'N', n, d3, rwork, work, ldu,
771  $ rwork( n+1 ), iinfo )
772  IF( iinfo.NE.0 ) THEN
773  WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n,
774  $ jtype, ioldsd
775  info = abs( iinfo )
776  IF( iinfo.LT.0 ) THEN
777  RETURN
778  ELSE
779  result( 6 ) = ulpinv
780  GO TO 150
781  END IF
782  END IF
783 *
784 *
785 * Do Tests 3 and 4 which are similar to 11 and 12 but with the
786 * D1 computed using the standard 1-stage reduction as reference
787 *
788  ntest = 6
789  temp1 = zero
790  temp2 = zero
791  temp3 = zero
792  temp4 = zero
793 *
794  DO 151 j = 1, n
795  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
796  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
797  temp3 = max( temp3, abs( d1( j ) ), abs( d3( j ) ) )
798  temp4 = max( temp4, abs( d1( j )-d3( j ) ) )
799  151 CONTINUE
800 *
801  result(5) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
802  result(6) = temp4 / max( unfl, ulp*max( temp3, temp4 ) )
803 *
804 * End of Loop -- Check for RESULT(j) > THRESH
805 *
806  150 CONTINUE
807  ntestt = ntestt + ntest
808 *
809 * Print out tests which fail.
810 *
811  DO 160 jr = 1, ntest
812  IF( result( jr ).GE.thresh ) THEN
813 *
814 * If this is the first test to fail,
815 * print a header to the data file.
816 *
817  IF( nerrs.EQ.0 ) THEN
818  WRITE( nounit, fmt = 9998 )'ZHB'
819  WRITE( nounit, fmt = 9997 )
820  WRITE( nounit, fmt = 9996 )
821  WRITE( nounit, fmt = 9995 )'Hermitian'
822  WRITE( nounit, fmt = 9994 )'unitary', '*',
823  $ 'conjugate transpose', ( '*', j = 1, 6 )
824  END IF
825  nerrs = nerrs + 1
826  WRITE( nounit, fmt = 9993 )n, k, ioldsd, jtype,
827  $ jr, result( jr )
828  END IF
829  160 CONTINUE
830 *
831  170 CONTINUE
832  180 CONTINUE
833  190 CONTINUE
834 *
835 * Summary
836 *
837  CALL dlasum( 'ZHB', nounit, nerrs, ntestt )
838  RETURN
839 *
840  9999 FORMAT( ' ZCHKHBSTG: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
841  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
842  9998 FORMAT( / 1x, a3,
843  $ ' -- Complex Hermitian Banded Tridiagonal Reduction Routines'
844  $ )
845  9997 FORMAT( ' Matrix types (see DCHK23 for details): ' )
846 *
847  9996 FORMAT( / ' Special Matrices:',
848  $ / ' 1=Zero matrix. ',
849  $ ' 5=Diagonal: clustered entries.',
850  $ / ' 2=Identity matrix. ',
851  $ ' 6=Diagonal: large, evenly spaced.',
852  $ / ' 3=Diagonal: evenly spaced entries. ',
853  $ ' 7=Diagonal: small, evenly spaced.',
854  $ / ' 4=Diagonal: geometr. spaced entries.' )
855  9995 FORMAT( ' Dense ', a, ' Banded Matrices:',
856  $ / ' 8=Evenly spaced eigenvals. ',
857  $ ' 12=Small, evenly spaced eigenvals.',
858  $ / ' 9=Geometrically spaced eigenvals. ',
859  $ ' 13=Matrix with random O(1) entries.',
860  $ / ' 10=Clustered eigenvalues. ',
861  $ ' 14=Matrix with large random entries.',
862  $ / ' 11=Large, evenly spaced eigenvals. ',
863  $ ' 15=Matrix with small random entries.' )
864 *
865  9994 FORMAT( / ' Tests performed: (S is Tridiag, U is ', a, ',',
866  $ / 20x, a, ' means ', a, '.', / ' UPLO=''U'':',
867  $ / ' 1= | A - U S U', a1, ' | / ( |A| n ulp ) ',
868  $ ' 2= | I - U U', a1, ' | / ( n ulp )', / ' UPLO=''L'':',
869  $ / ' 3= | A - U S U', a1, ' | / ( |A| n ulp ) ',
870  $ ' 4= | I - U U', a1, ' | / ( n ulp )' / ' Eig check:',
871  $ /' 5= | D1 - D2', '', ' | / ( |D1| ulp ) ',
872  $ ' 6= | D1 - D3', '', ' | / ( |D1| ulp ) ' )
873  9993 FORMAT( ' N=', i5, ', K=', i4, ', seed=', 4( i4, ',' ), ' type ',
874  $ i2, ', test(', i2, ')=', g10.3 )
875 *
876 * End of ZCHKHBSTG
877 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
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:492
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:134
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:112
subroutine zhbt21(UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK, RWORK, RESULT)
ZHBT21
Definition: zhbt21.f:152
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
ZLATMS
Definition: zlatms.f:334
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:108
subroutine dlasum(TYPE, IOUNIT, IE, NRUN)
DLASUM
Definition: dlasum.f:45
subroutine zhbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
ZHBTRD
Definition: zhbtrd.f:165
Here is the call graph for this function:
Here is the caller graph for this function: