409 SUBROUTINE schkhs( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
410 $ NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, WR1,
411 $ WI1, WR2, WI2, WR3, WI3, EVECTL, EVECTR,
412 $ EVECTY, EVECTX, UU, TAU, WORK, NWORK, IWORK,
413 $ SELECT, RESULT, INFO )
421 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK
425 LOGICAL DOTYPE( * ), SELECT( * )
426 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
427 REAL A( LDA, * ), EVECTL( LDU, * ),
428 $ EVECTR( LDU, * ), EVECTX( LDU, * ),
429 $ evecty( ldu, * ), h( lda, * ), result( 14 ),
430 $ t1( lda, * ), t2( lda, * ), tau( * ),
431 $ u( ldu, * ), uu( ldu, * ), uz( ldu, * ),
432 $ wi1( * ), wi2( * ), wi3( * ), work( * ),
433 $ wr1( * ), wr2( * ), wr3( * ), z( ldu, * )
440 PARAMETER ( ZERO = 0.0, one = 1.0 )
442 PARAMETER ( MAXTYP = 21 )
446 INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL,
447 $ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS,
448 $ NMATS, NMAX, NSELC, NSELR, NTEST, NTESTT
449 REAL ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP,
450 $ rtulpi, rtunfl, temp1, temp2, ulp, ulpinv, unfl
453 CHARACTER ADUMMA( 1 )
454 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
455 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
470 INTRINSIC abs, max, min, real, sqrt
473 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
474 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
476 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
477 $ 1, 5, 5, 5, 4, 3, 1 /
478 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
490 nmax = max( nmax, nn( j ) )
497 IF( nsizes.LT.0 )
THEN
499 ELSE IF( badnn )
THEN
501 ELSE IF( ntypes.LT.0 )
THEN
503 ELSE IF( thresh.LT.zero )
THEN
505 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
507 ELSE IF( ldu.LE.1 .OR. ldu.LT.nmax )
THEN
509 ELSE IF( 4*nmax*nmax+2.GT.nwork )
THEN
514 CALL xerbla(
'SCHKHS', -info )
520 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
525 unfl = slamch(
'Safe minimum' )
526 ovfl = slamch(
'Overflow' )
528 ulp = slamch(
'Epsilon' )*slamch(
'Base' )
530 rtunfl = sqrt( unfl )
531 rtovfl = sqrt( ovfl )
540 DO 270 jsize = 1, nsizes
545 aninv = one / real( n1 )
547 IF( nsizes.NE.1 )
THEN
548 mtypes = min( maxtyp, ntypes )
550 mtypes = min( maxtyp+1, ntypes )
553 DO 260 jtype = 1, mtypes
554 IF( .NOT.dotype( jtype ) )
562 ioldsd( j ) = iseed( j )
587 IF( mtypes.GT.maxtyp )
590 itype = ktype( jtype )
591 imode = kmode( jtype )
595 GO TO ( 40, 50, 60 )kmagn( jtype )
602 anorm = ( rtovfl*ulp )*aninv
606 anorm = rtunfl*n*ulpinv
611 CALL slaset(
'Full', lda, n, zero, zero, a, lda )
617 IF( itype.EQ.1 )
THEN
623 ELSE IF( itype.EQ.2 )
THEN
628 a( jcol, jcol ) = anorm
631 ELSE IF( itype.EQ.3 )
THEN
636 a( jcol, jcol ) = anorm
638 $ a( jcol, jcol-1 ) = one
641 ELSE IF( itype.EQ.4 )
THEN
645 CALL slatms( n, n,
'S', iseed,
'S', work, imode, cond,
646 $ anorm, 0, 0,
'N', a, lda, work( n+1 ),
649 ELSE IF( itype.EQ.5 )
THEN
653 CALL slatms( n, n,
'S', iseed,
'S', work, imode, cond,
654 $ anorm, n, n,
'N', a, lda, work( n+1 ),
657 ELSE IF( itype.EQ.6 )
THEN
661 IF( kconds( jtype ).EQ.1 )
THEN
663 ELSE IF( kconds( jtype ).EQ.2 )
THEN
670 CALL slatme( n,
'S', iseed, work, imode, cond, one,
671 $ adumma,
'T',
'T',
'T', work( n+1 ), 4,
672 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
675 ELSE IF( itype.EQ.7 )
THEN
679 CALL slatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
680 $
'T',
'N', work( n+1 ), 1, one,
681 $ work( 2*n+1 ), 1, one,
'N', idumma, 0, 0,
682 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
684 ELSE IF( itype.EQ.8 )
THEN
688 CALL slatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
689 $
'T',
'N', work( n+1 ), 1, one,
690 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
691 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
693 ELSE IF( itype.EQ.9 )
THEN
697 CALL slatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
698 $
'T',
'N', work( n+1 ), 1, one,
699 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
700 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
702 ELSE IF( itype.EQ.10 )
THEN
706 CALL slatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
707 $
'T',
'N', work( n+1 ), 1, one,
708 $ work( 2*n+1 ), 1, one,
'N', idumma, n, 0,
709 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
716 IF( iinfo.NE.0 )
THEN
717 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
727 CALL slacpy(
' ', n, n, a, lda, h, lda )
734 CALL sgehrd( n, ilo, ihi, h, lda, work, work( n+1 ),
737 IF( iinfo.NE.0 )
THEN
739 WRITE( nounit, fmt = 9999 )
'SGEHRD', iinfo, n, jtype,
748 u( i, j ) = h( i, j )
749 uu( i, j ) = h( i, j )
753 CALL scopy( n-1, work, 1, tau, 1 )
754 CALL sorghr( n, ilo, ihi, u, ldu, work, work( n+1 ),
758 CALL shst01( n, ilo, ihi, a, lda, h, lda, u, ldu, work,
759 $ nwork, result( 1 ) )
765 CALL slacpy(
' ', n, n, h, lda, t2, lda )
769 CALL shseqr(
'E',
'N', n, ilo, ihi, t2, lda, wr3, wi3, uz,
770 $ ldu, work, nwork, iinfo )
771 IF( iinfo.NE.0 )
THEN
772 WRITE( nounit, fmt = 9999 )
'SHSEQR(E)', iinfo, n, jtype,
774 IF( iinfo.LE.n+2 )
THEN
782 CALL slacpy(
' ', n, n, h, lda, t2, lda )
784 CALL shseqr(
'S',
'N', n, ilo, ihi, t2, lda, wr2, wi2, uz,
785 $ ldu, work, nwork, iinfo )
786 IF( iinfo.NE.0 .AND. iinfo.LE.n+2 )
THEN
787 WRITE( nounit, fmt = 9999 )
'SHSEQR(S)', iinfo, n, jtype,
796 CALL slacpy(
' ', n, n, h, lda, t1, lda )
797 CALL slacpy(
' ', n, n, u, ldu, uz, ldu )
799 CALL shseqr(
'S',
'V', n, ilo, ihi, t1, lda, wr1, wi1, uz,
800 $ ldu, work, nwork, iinfo )
801 IF( iinfo.NE.0 .AND. iinfo.LE.n+2 )
THEN
802 WRITE( nounit, fmt = 9999 )
'SHSEQR(V)', iinfo, n, jtype,
810 CALL sgemm(
'T',
'N', n, n, n, one, u, ldu, uz, ldu, zero,
817 CALL shst01( n, ilo, ihi, h, lda, t1, lda, z, ldu, work,
818 $ nwork, result( 3 ) )
823 CALL shst01( n, ilo, ihi, a, lda, t1, lda, uz, ldu, work,
824 $ nwork, result( 5 ) )
828 CALL sget10( n, n, t2, lda, t1, lda, work, result( 7 ) )
835 temp1 = max( temp1, abs( wr1( j ) )+abs( wi1( j ) ),
836 $ abs( wr2( j ) )+abs( wi2( j ) ) )
837 temp2 = max( temp2, abs( wr1( j )-wr2( j ) )+
838 $ abs( wi1( j )-wi2( j ) ) )
841 result( 8 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
856 IF( wi1( j ).EQ.zero )
THEN
857 IF( nselr.LT.max( n / 4, 1 ) )
THEN
861 SELECT( j ) = .false.
865 IF( nselc.LT.max( n / 4, 1 ) )
THEN
868 SELECT( j-1 ) = .false.
870 SELECT( j ) = .false.
871 SELECT( j-1 ) = .false.
878 CALL strevc(
'Right',
'All',
SELECT, n, t1, lda, dumma, ldu,
879 $ evectr, ldu, n, in, work, iinfo )
880 IF( iinfo.NE.0 )
THEN
881 WRITE( nounit, fmt = 9999 )
'STREVC(R,A)', iinfo, n,
889 CALL sget22(
'N',
'N',
'N', n, t1, lda, evectr, ldu, wr1,
890 $ wi1, work, dumma( 1 ) )
891 result( 9 ) = dumma( 1 )
892 IF( dumma( 2 ).GT.thresh )
THEN
893 WRITE( nounit, fmt = 9998 )
'Right',
'STREVC',
894 $ dumma( 2 ), n, jtype, ioldsd
900 CALL strevc(
'Right',
'Some',
SELECT, n, t1, lda, dumma,
901 $ ldu, evectl, ldu, n, in, work, iinfo )
902 IF( iinfo.NE.0 )
THEN
903 WRITE( nounit, fmt = 9999 )
'STREVC(R,S)', iinfo, n,
912 IF(
SELECT( j ) .AND. wi1( j ).EQ.zero )
THEN
914 IF( evectr( jj, j ).NE.evectl( jj, k ) )
THEN
920 ELSE IF(
SELECT( j ) .AND. wi1( j ).NE.zero )
THEN
922 IF( evectr( jj, j ).NE.evectl( jj, k ) .OR.
923 $ evectr( jj, j+1 ).NE.evectl( jj, k+1 ) )
THEN
933 $
WRITE( nounit, fmt = 9997 )
'Right',
'STREVC', n, jtype,
939 result( 10 ) = ulpinv
940 CALL strevc(
'Left',
'All',
SELECT, n, t1, lda, evectl, ldu,
941 $ dumma, ldu, n, in, work, iinfo )
942 IF( iinfo.NE.0 )
THEN
943 WRITE( nounit, fmt = 9999 )
'STREVC(L,A)', iinfo, n,
951 CALL sget22(
'Trans',
'N',
'Conj', n, t1, lda, evectl, ldu,
952 $ wr1, wi1, work, dumma( 3 ) )
953 result( 10 ) = dumma( 3 )
954 IF( dumma( 4 ).GT.thresh )
THEN
955 WRITE( nounit, fmt = 9998 )
'Left',
'STREVC', dumma( 4 ),
962 CALL strevc(
'Left',
'Some',
SELECT, n, t1, lda, evectr,
963 $ ldu, dumma, ldu, n, in, work, iinfo )
964 IF( iinfo.NE.0 )
THEN
965 WRITE( nounit, fmt = 9999 )
'STREVC(L,S)', iinfo, n,
974 IF(
SELECT( j ) .AND. wi1( j ).EQ.zero )
THEN
976 IF( evectl( jj, j ).NE.evectr( jj, k ) )
THEN
982 ELSE IF(
SELECT( j ) .AND. wi1( j ).NE.zero )
THEN
984 IF( evectl( jj, j ).NE.evectr( jj, k ) .OR.
985 $ evectl( jj, j+1 ).NE.evectr( jj, k+1 ) )
THEN
995 $
WRITE( nounit, fmt = 9997 )
'Left',
'STREVC', n, jtype,
1001 result( 11 ) = ulpinv
1003 SELECT( j ) = .true.
1006 CALL shsein(
'Right',
'Qr',
'Ninitv',
SELECT, n, h, lda,
1007 $ wr3, wi3, dumma, ldu, evectx, ldu, n1, in,
1008 $ work, iwork, iwork, iinfo )
1009 IF( iinfo.NE.0 )
THEN
1010 WRITE( nounit, fmt = 9999 )
'SHSEIN(R)', iinfo, n, jtype,
1021 CALL sget22(
'N',
'N',
'N', n, h, lda, evectx, ldu, wr3,
1022 $ wi3, work, dumma( 1 ) )
1023 IF( dumma( 1 ).LT.ulpinv )
1024 $ result( 11 ) = dumma( 1 )*aninv
1025 IF( dumma( 2 ).GT.thresh )
THEN
1026 WRITE( nounit, fmt = 9998 )
'Right',
'SHSEIN',
1027 $ dumma( 2 ), n, jtype, ioldsd
1034 result( 12 ) = ulpinv
1036 SELECT( j ) = .true.
1039 CALL shsein(
'Left',
'Qr',
'Ninitv',
SELECT, n, h, lda, wr3,
1040 $ wi3, evecty, ldu, dumma, ldu, n1, in, work,
1041 $ iwork, iwork, iinfo )
1042 IF( iinfo.NE.0 )
THEN
1043 WRITE( nounit, fmt = 9999 )
'SHSEIN(L)', iinfo, n, jtype,
1054 CALL sget22(
'C',
'N',
'C', n, h, lda, evecty, ldu, wr3,
1055 $ wi3, work, dumma( 3 ) )
1056 IF( dumma( 3 ).LT.ulpinv )
1057 $ result( 12 ) = dumma( 3 )*aninv
1058 IF( dumma( 4 ).GT.thresh )
THEN
1059 WRITE( nounit, fmt = 9998 )
'Left',
'SHSEIN',
1060 $ dumma( 4 ), n, jtype, ioldsd
1067 result( 13 ) = ulpinv
1069 CALL sormhr(
'Left',
'No transpose', n, n, ilo, ihi, uu,
1070 $ ldu, tau, evectx, ldu, work, nwork, iinfo )
1071 IF( iinfo.NE.0 )
THEN
1072 WRITE( nounit, fmt = 9999 )
'SORMHR(R)', iinfo, n, jtype,
1083 CALL sget22(
'N',
'N',
'N', n, a, lda, evectx, ldu, wr3,
1084 $ wi3, work, dumma( 1 ) )
1085 IF( dumma( 1 ).LT.ulpinv )
1086 $ result( 13 ) = dumma( 1 )*aninv
1092 result( 14 ) = ulpinv
1094 CALL sormhr(
'Left',
'No transpose', n, n, ilo, ihi, uu,
1095 $ ldu, tau, evecty, ldu, work, nwork, iinfo )
1096 IF( iinfo.NE.0 )
THEN
1097 WRITE( nounit, fmt = 9999 )
'SORMHR(L)', iinfo, n, jtype,
1108 CALL sget22(
'C',
'N',
'C', n, a, lda, evecty, ldu, wr3,
1109 $ wi3, work, dumma( 3 ) )
1110 IF( dumma( 3 ).LT.ulpinv )
1111 $ result( 14 ) = dumma( 3 )*aninv
1118 ntestt = ntestt + ntest
1119 CALL slafts(
'SHS', n, n, jtype, ntest, result, ioldsd,
1120 $ thresh, nounit, nerrs )
1127 CALL slasum(
'SHS', nounit, nerrs, ntestt )
1131 9999
FORMAT(
' SCHKHS: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
1132 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
1133 9998
FORMAT(
' SCHKHS: ', a,
' Eigenvectors from ', a,
' incorrectly ',
1134 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
1135 $
'N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
1137 9997
FORMAT(
' SCHKHS: Selected ', a,
' Eigenvectors from ', a,
1138 $
' do not match other eigenvectors ', 9x,
'N=', i6,
1139 $
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )