358 SUBROUTINE ddrvsg2stg( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
359 $ NOUNIT, A, LDA, B, LDB, D, D2, Z, LDZ, AB,
360 $ BB, AP, BP, WORK, NWORK, IWORK, LIWORK,
370 INTEGER INFO, LDA, LDB, LDZ, LIWORK, NOUNIT, NSIZES,
372 DOUBLE PRECISION THRESH
376 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
377 DOUBLE PRECISION A( LDA, * ), AB( LDA, * ), AP( * ),
378 $ b( ldb, * ), bb( ldb, * ), bp( * ), d( * ),
379 $ d2( * ), result( * ), work( * ), z( ldz, * )
385 DOUBLE PRECISION ZERO, ONE, TEN
386 PARAMETER ( ZERO = 0.0d0, one = 1.0d0, ten = 10.0d0 )
388 parameter( maxtyp = 21 )
393 INTEGER I, IBTYPE, IBUPLO, IINFO, IJ, IL, IMODE, ITEMP,
394 $ itype, iu, j, jcol, jsize, jtype, ka, ka9, kb,
395 $ kb9, m, mtypes, n, nerrs, nmats, nmax, ntest,
397 DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
398 $ RTUNFL, ULP, ULPINV, UNFL, VL, VU, TEMP1, TEMP2
401 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
402 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
407 DOUBLE PRECISION DLAMCH, DLARND
408 EXTERNAL LSAME, DLAMCH, DLARND
417 INTRINSIC abs, dble, max, min, sqrt
420 DATA ktype / 1, 2, 5*4, 5*5, 3*8, 6*9 /
421 DATA kmagn / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
423 DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
436 nmax = max( nmax, nn( j ) )
443 IF( nsizes.LT.0 )
THEN
445 ELSE IF( badnn )
THEN
447 ELSE IF( ntypes.LT.0 )
THEN
449 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
451 ELSE IF( ldz.LE.1 .OR. ldz.LT.nmax )
THEN
453 ELSE IF( 2*max( nmax, 3 )**2.GT.nwork )
THEN
455 ELSE IF( 2*max( nmax, 3 )**2.GT.liwork )
THEN
460 CALL xerbla(
'DDRVSG2STG', -info )
466 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
471 unfl = dlamch(
'Safe minimum' )
472 ovfl = dlamch(
'Overflow' )
473 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
475 rtunfl = sqrt( unfl )
476 rtovfl = sqrt( ovfl )
479 iseed2( i ) = iseed( i )
487 DO 650 jsize = 1, nsizes
489 aninv = one / dble( max( 1, n ) )
491 IF( nsizes.NE.1 )
THEN
492 mtypes = min( maxtyp, ntypes )
494 mtypes = min( maxtyp+1, ntypes )
499 DO 640 jtype = 1, mtypes
500 IF( .NOT.dotype( jtype ) )
506 ioldsd( j ) = iseed( j )
524 IF( mtypes.GT.maxtyp )
527 itype = ktype( jtype )
528 imode = kmode( jtype )
532 GO TO ( 40, 50, 60 )kmagn( jtype )
539 anorm = ( rtovfl*ulp )*aninv
543 anorm = rtunfl*n*ulpinv
553 IF( itype.EQ.1 )
THEN
559 CALL dlaset(
'Full', lda, n, zero, zero, a, lda )
561 ELSE IF( itype.EQ.2 )
THEN
567 CALL dlaset(
'Full', lda, n, zero, zero, a, lda )
569 a( jcol, jcol ) = anorm
572 ELSE IF( itype.EQ.4 )
THEN
578 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
579 $ anorm, 0, 0,
'N', a, lda, work( n+1 ),
582 ELSE IF( itype.EQ.5 )
THEN
588 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
589 $ anorm, n, n,
'N', a, lda, work( n+1 ),
592 ELSE IF( itype.EQ.7 )
THEN
598 CALL dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
599 $
'T',
'N', work( n+1 ), 1, one,
600 $ work( 2*n+1 ), 1, one,
'N', idumma, 0, 0,
601 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
603 ELSE IF( itype.EQ.8 )
THEN
609 CALL dlatmr( n, n,
'S', iseed,
'H', work, 6, one, one,
610 $
'T',
'N', work( n+1 ), 1, one,
611 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
612 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
614 ELSE IF( itype.EQ.9 )
THEN
628 IF( kb9.GT.ka9 )
THEN
632 ka = max( 0, min( n-1, ka9 ) )
633 kb = max( 0, min( n-1, kb9 ) )
634 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
635 $ anorm, ka, ka,
'N', a, lda, work( n+1 ),
643 IF( iinfo.NE.0 )
THEN
644 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
657 il = 1 + int( ( n-1 )*dlarnd( 1, iseed2 ) )
658 iu = 1 + int( ( n-1 )*dlarnd( 1, iseed2 ) )
687 CALL dlatms( n, n,
'U', iseed,
'P', work, 5, ten, one,
688 $ kb, kb, uplo, b, ldb, work( n+1 ),
695 CALL dlacpy(
' ', n, n, a, lda, z, ldz )
696 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
698 CALL dsygv( ibtype,
'V', uplo, n, z, ldz, bb, ldb, d,
699 $ work, nwork, iinfo )
700 IF( iinfo.NE.0 )
THEN
701 WRITE( nounit, fmt = 9999 )
'DSYGV(V,' // uplo //
702 $
')', iinfo, n, jtype, ioldsd
704 IF( iinfo.LT.0 )
THEN
707 result( ntest ) = ulpinv
714 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
715 $ ldz, d, work, result( ntest ) )
721 CALL dlacpy(
' ', n, n, a, lda, z, ldz )
722 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
725 $ bb, ldb, d2, work, nwork, iinfo )
726 IF( iinfo.NE.0 )
THEN
727 WRITE( nounit, fmt = 9999 )
728 $
'DSYGV_2STAGE(V,' // uplo //
729 $
')', iinfo, n, jtype, ioldsd
731 IF( iinfo.LT.0 )
THEN
734 result( ntest ) = ulpinv
751 temp1 = max( temp1, abs( d( j ) ),
753 temp2 = max( temp2, abs( d( j )-d2( j ) ) )
756 result( ntest ) = temp2 /
757 $ max( unfl, ulp*max( temp1, temp2 ) )
763 CALL dlacpy(
' ', n, n, a, lda, z, ldz )
764 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
766 CALL dsygvd( ibtype,
'V', uplo, n, z, ldz, bb, ldb, d,
767 $ work, nwork, iwork, liwork, iinfo )
768 IF( iinfo.NE.0 )
THEN
769 WRITE( nounit, fmt = 9999 )
'DSYGVD(V,' // uplo //
770 $
')', iinfo, n, jtype, ioldsd
772 IF( iinfo.LT.0 )
THEN
775 result( ntest ) = ulpinv
782 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
783 $ ldz, d, work, result( ntest ) )
789 CALL dlacpy(
' ', n, n, a, lda, ab, lda )
790 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
792 CALL dsygvx( ibtype,
'V',
'A', uplo, n, ab, lda, bb,
793 $ ldb, vl, vu, il, iu, abstol, m, d, z,
794 $ ldz, work, nwork, iwork( n+1 ), iwork,
796 IF( iinfo.NE.0 )
THEN
797 WRITE( nounit, fmt = 9999 )
'DSYGVX(V,A' // uplo //
798 $
')', iinfo, n, jtype, ioldsd
800 IF( iinfo.LT.0 )
THEN
803 result( ntest ) = ulpinv
810 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
811 $ ldz, d, work, result( ntest ) )
815 CALL dlacpy(
' ', n, n, a, lda, ab, lda )
816 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
825 CALL dsygvx( ibtype,
'V',
'V', uplo, n, ab, lda, bb,
826 $ ldb, vl, vu, il, iu, abstol, m, d, z,
827 $ ldz, work, nwork, iwork( n+1 ), iwork,
829 IF( iinfo.NE.0 )
THEN
830 WRITE( nounit, fmt = 9999 )
'DSYGVX(V,V,' //
831 $ uplo //
')', iinfo, n, jtype, ioldsd
833 IF( iinfo.LT.0 )
THEN
836 result( ntest ) = ulpinv
843 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
844 $ ldz, d, work, result( ntest ) )
848 CALL dlacpy(
' ', n, n, a, lda, ab, lda )
849 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
851 CALL dsygvx( ibtype,
'V',
'I', uplo, n, ab, lda, bb,
852 $ ldb, vl, vu, il, iu, abstol, m, d, z,
853 $ ldz, work, nwork, iwork( n+1 ), iwork,
855 IF( iinfo.NE.0 )
THEN
856 WRITE( nounit, fmt = 9999 )
'DSYGVX(V,I,' //
857 $ uplo //
')', iinfo, n, jtype, ioldsd
859 IF( iinfo.LT.0 )
THEN
862 result( ntest ) = ulpinv
869 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
870 $ ldz, d, work, result( ntest ) )
880 IF( lsame( uplo,
'U' ) )
THEN
900 CALL dspgv( ibtype,
'V', uplo, n, ap, bp, d, z, ldz,
902 IF( iinfo.NE.0 )
THEN
903 WRITE( nounit, fmt = 9999 )
'DSPGV(V,' // uplo //
904 $
')', iinfo, n, jtype, ioldsd
906 IF( iinfo.LT.0 )
THEN
909 result( ntest ) = ulpinv
916 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
917 $ ldz, d, work, result( ntest ) )
925 IF( lsame( uplo,
'U' ) )
THEN
945 CALL dspgvd( ibtype,
'V', uplo, n, ap, bp, d, z, ldz,
946 $ work, nwork, iwork, liwork, iinfo )
947 IF( iinfo.NE.0 )
THEN
948 WRITE( nounit, fmt = 9999 )
'DSPGVD(V,' // uplo //
949 $
')', iinfo, n, jtype, ioldsd
951 IF( iinfo.LT.0 )
THEN
954 result( ntest ) = ulpinv
961 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
962 $ ldz, d, work, result( ntest ) )
970 IF( lsame( uplo,
'U' ) )
THEN
990 CALL dspgvx( ibtype,
'V',
'A', uplo, n, ap, bp, vl,
991 $ vu, il, iu, abstol, m, d, z, ldz, work,
992 $ iwork( n+1 ), iwork, info )
993 IF( iinfo.NE.0 )
THEN
994 WRITE( nounit, fmt = 9999 )
'DSPGVX(V,A' // uplo //
995 $
')', iinfo, n, jtype, ioldsd
997 IF( iinfo.LT.0 )
THEN
1000 result( ntest ) = ulpinv
1007 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
1008 $ ldz, d, work, result( ntest ) )
1014 IF( lsame( uplo,
'U' ) )
THEN
1018 ap( ij ) = a( i, j )
1019 bp( ij ) = b( i, j )
1027 ap( ij ) = a( i, j )
1028 bp( ij ) = b( i, j )
1036 CALL dspgvx( ibtype,
'V',
'V', uplo, n, ap, bp, vl,
1037 $ vu, il, iu, abstol, m, d, z, ldz, work,
1038 $ iwork( n+1 ), iwork, info )
1039 IF( iinfo.NE.0 )
THEN
1040 WRITE( nounit, fmt = 9999 )
'DSPGVX(V,V' // uplo //
1041 $
')', iinfo, n, jtype, ioldsd
1043 IF( iinfo.LT.0 )
THEN
1046 result( ntest ) = ulpinv
1053 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
1054 $ ldz, d, work, result( ntest ) )
1060 IF( lsame( uplo,
'U' ) )
THEN
1064 ap( ij ) = a( i, j )
1065 bp( ij ) = b( i, j )
1073 ap( ij ) = a( i, j )
1074 bp( ij ) = b( i, j )
1080 CALL dspgvx( ibtype,
'V',
'I', uplo, n, ap, bp, vl,
1081 $ vu, il, iu, abstol, m, d, z, ldz, work,
1082 $ iwork( n+1 ), iwork, info )
1083 IF( iinfo.NE.0 )
THEN
1084 WRITE( nounit, fmt = 9999 )
'DSPGVX(V,I' // uplo //
1085 $
')', iinfo, n, jtype, ioldsd
1087 IF( iinfo.LT.0 )
THEN
1090 result( ntest ) = ulpinv
1097 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
1098 $ ldz, d, work, result( ntest ) )
1102 IF( ibtype.EQ.1 )
THEN
1110 IF( lsame( uplo,
'U' ) )
THEN
1112 DO 320 i = max( 1, j-ka ), j
1113 ab( ka+1+i-j, j ) = a( i, j )
1115 DO 330 i = max( 1, j-kb ), j
1116 bb( kb+1+i-j, j ) = b( i, j )
1121 DO 350 i = j, min( n, j+ka )
1122 ab( 1+i-j, j ) = a( i, j )
1124 DO 360 i = j, min( n, j+kb )
1125 bb( 1+i-j, j ) = b( i, j )
1130 CALL dsbgv(
'V', uplo, n, ka, kb, ab, lda, bb, ldb,
1131 $ d, z, ldz, work, iinfo )
1132 IF( iinfo.NE.0 )
THEN
1133 WRITE( nounit, fmt = 9999 )
'DSBGV(V,' //
1134 $ uplo //
')', iinfo, n, jtype, ioldsd
1136 IF( iinfo.LT.0 )
THEN
1139 result( ntest ) = ulpinv
1146 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
1147 $ ldz, d, work, result( ntest ) )
1155 IF( lsame( uplo,
'U' ) )
THEN
1157 DO 380 i = max( 1, j-ka ), j
1158 ab( ka+1+i-j, j ) = a( i, j )
1160 DO 390 i = max( 1, j-kb ), j
1161 bb( kb+1+i-j, j ) = b( i, j )
1166 DO 410 i = j, min( n, j+ka )
1167 ab( 1+i-j, j ) = a( i, j )
1169 DO 420 i = j, min( n, j+kb )
1170 bb( 1+i-j, j ) = b( i, j )
1175 CALL dsbgvd(
'V', uplo, n, ka, kb, ab, lda, bb,
1176 $ ldb, d, z, ldz, work, nwork, iwork,
1178 IF( iinfo.NE.0 )
THEN
1179 WRITE( nounit, fmt = 9999 )
'DSBGVD(V,' //
1180 $ uplo //
')', iinfo, n, jtype, ioldsd
1182 IF( iinfo.LT.0 )
THEN
1185 result( ntest ) = ulpinv
1192 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
1193 $ ldz, d, work, result( ntest ) )
1201 IF( lsame( uplo,
'U' ) )
THEN
1203 DO 440 i = max( 1, j-ka ), j
1204 ab( ka+1+i-j, j ) = a( i, j )
1206 DO 450 i = max( 1, j-kb ), j
1207 bb( kb+1+i-j, j ) = b( i, j )
1212 DO 470 i = j, min( n, j+ka )
1213 ab( 1+i-j, j ) = a( i, j )
1215 DO 480 i = j, min( n, j+kb )
1216 bb( 1+i-j, j ) = b( i, j )
1221 CALL dsbgvx(
'V',
'A', uplo, n, ka, kb, ab, lda,
1222 $ bb, ldb, bp, max( 1, n ), vl, vu, il,
1223 $ iu, abstol, m, d, z, ldz, work,
1224 $ iwork( n+1 ), iwork, iinfo )
1225 IF( iinfo.NE.0 )
THEN
1226 WRITE( nounit, fmt = 9999 )
'DSBGVX(V,A' //
1227 $ uplo //
')', iinfo, n, jtype, ioldsd
1229 IF( iinfo.LT.0 )
THEN
1232 result( ntest ) = ulpinv
1239 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
1240 $ ldz, d, work, result( ntest ) )
1247 IF( lsame( uplo,
'U' ) )
THEN
1249 DO 500 i = max( 1, j-ka ), j
1250 ab( ka+1+i-j, j ) = a( i, j )
1252 DO 510 i = max( 1, j-kb ), j
1253 bb( kb+1+i-j, j ) = b( i, j )
1258 DO 530 i = j, min( n, j+ka )
1259 ab( 1+i-j, j ) = a( i, j )
1261 DO 540 i = j, min( n, j+kb )
1262 bb( 1+i-j, j ) = b( i, j )
1269 CALL dsbgvx(
'V',
'V', uplo, n, ka, kb, ab, lda,
1270 $ bb, ldb, bp, max( 1, n ), vl, vu, il,
1271 $ iu, abstol, m, d, z, ldz, work,
1272 $ iwork( n+1 ), iwork, iinfo )
1273 IF( iinfo.NE.0 )
THEN
1274 WRITE( nounit, fmt = 9999 )
'DSBGVX(V,V' //
1275 $ uplo //
')', iinfo, n, jtype, ioldsd
1277 IF( iinfo.LT.0 )
THEN
1280 result( ntest ) = ulpinv
1287 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
1288 $ ldz, d, work, result( ntest ) )
1294 IF( lsame( uplo,
'U' ) )
THEN
1296 DO 560 i = max( 1, j-ka ), j
1297 ab( ka+1+i-j, j ) = a( i, j )
1299 DO 570 i = max( 1, j-kb ), j
1300 bb( kb+1+i-j, j ) = b( i, j )
1305 DO 590 i = j, min( n, j+ka )
1306 ab( 1+i-j, j ) = a( i, j )
1308 DO 600 i = j, min( n, j+kb )
1309 bb( 1+i-j, j ) = b( i, j )
1314 CALL dsbgvx(
'V',
'I', uplo, n, ka, kb, ab, lda,
1315 $ bb, ldb, bp, max( 1, n ), vl, vu, il,
1316 $ iu, abstol, m, d, z, ldz, work,
1317 $ iwork( n+1 ), iwork, iinfo )
1318 IF( iinfo.NE.0 )
THEN
1319 WRITE( nounit, fmt = 9999 )
'DSBGVX(V,I' //
1320 $ uplo //
')', iinfo, n, jtype, ioldsd
1322 IF( iinfo.LT.0 )
THEN
1325 result( ntest ) = ulpinv
1332 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
1333 $ ldz, d, work, result( ntest ) )
1342 ntestt = ntestt + ntest
1343 CALL dlafts(
'DSG', n, n, jtype, ntest, result, ioldsd,
1344 $ thresh, nounit, nerrs )
1350 CALL dlasum(
'DSG', nounit, nerrs, ntestt )
1356 9999
FORMAT(
' DDRVSG2STG: ', a,
' returned INFO=', i6,
'.', / 9x,
1357 $
'N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )