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, , 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' )
474 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
476 rtunfl = sqrt( unfl )
477 rtovfl = sqrt( ovfl )
480 iseed2( i ) = iseed( i )
488 DO 650 jsize = 1, nsizes
490 aninv = one / dble(
max( 1, n ) )
492 IF( nsizes.NE.1 )
THEN
493 mtypes =
min( maxtyp, ntypes )
495 mtypes =
min( maxtyp+1, ntypes )
500 DO 640 jtype = 1, mtypes
501 IF( .NOT.dotype( jtype ) )
507 ioldsd( j ) = iseed( j )
525 IF( mtypes.GT.maxtyp )
528 itype = ktype( jtype )
529 imode = kmode( jtype )
533 GO TO ( 40, 50, 60 )kmagn( jtype )
540 anorm = ( rtovfl*ulp )*aninv
544 anorm = rtunfl*n*ulpinv
554 IF( itype.EQ.1 )
THEN
560 CALL dlaset(
'Full', lda, n, zero, zero, a, lda )
562 ELSE IF( itype.EQ.2 )
THEN
568 CALL dlaset(
'Full', lda, n, zero, zero, a, lda )
570 a( jcol, jcol ) = anorm
573 ELSE IF( itype.EQ.4 )
THEN
579 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
580 $ anorm, 0, 0,
'N', a, lda, work( n+1 ),
583 ELSE IF( itype.EQ.5 )
THEN
589 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
590 $ anorm, n, n,
'N', a, lda, work( n+1 ),
593 ELSE IF( itype.EQ.7 )
THEN
599 CALL dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
600 $
'T',
'N', work( n+1 ), 1, one,
601 $ work( 2*n+1 ), 1, one,
'N', idumma, 0, 0,
602 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
604 ELSE IF( itype.EQ.8 )
THEN
610 CALL dlatmr( n, n,
'S', iseed,
'H', work, 6, one, one,
611 $
'T',
'N', work( n+1 ), 1, one,
612 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
613 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
615 ELSE IF( itype.EQ.9 )
THEN
629 IF( kb9.GT.ka9 )
THEN
633 ka =
max( 0,
min( n-1, ka9 ) )
634 kb =
max( 0,
min( n-1, kb9 ) )
635 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
636 $ anorm, ka, ka,
'N', a, lda, work( n+1 ),
644 IF( iinfo.NE.0 )
THEN
645 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
658 il = 1 + int( ( n-1 )*dlarnd( 1, iseed2 ) )
659 iu = 1 + int( ( n-1 )*dlarnd( 1, iseed2 ) )
688 CALL dlatms( n, n,
'U', iseed,
'P', work, 5, ten, one,
689 $ kb, kb, uplo, b, ldb, work( n+1 ),
696 CALL dlacpy(
' ', n, n, a, lda, z, ldz )
697 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
699 CALL dsygv( ibtype,
'V', uplo, n, z, ldz, bb, ldb, d,
700 $ work, nwork, iinfo )
701 IF( iinfo.NE.0 )
THEN
702 WRITE( nounit, fmt = 9999 )
'DSYGV(V,' // uplo //
703 $
')', iinfo, n, jtype, ioldsd
705 IF( iinfo.LT.0 )
THEN
708 result( ntest ) = ulpinv
715 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
716 $ ldz, d, work, result( ntest ) )
722 CALL dlacpy(
' ', n, n, a, lda, z, ldz )
723 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
726 $ bb, ldb, d2, work, nwork, iinfo )
727 IF( iinfo.NE.0 )
THEN
728 WRITE( nounit, fmt = 9999 )
729 $
'DSYGV_2STAGE(V,' // uplo //
730 $
')', iinfo, n, jtype, ioldsd
732 IF( iinfo.LT.0 )
THEN
735 result( ntest ) = ulpinv
752 temp1 =
max( temp1, abs( d( j ) ),
754 temp2 =
max( temp2, abs( d( j )-d2( j ) ) )
757 result( ntest ) = temp2 /
758 $
max( unfl, ulp*
max( temp1, temp2 ) )
764 CALL dlacpy(
' ', n, n, a, lda, z, ldz )
765 CALL dlacpy( uplo, n, n, b
767 CALL dsygvd( ibtype,
'V', uplo, n, z, ldz, bb, ldb, d,
768 $ work, nwork, iwork, liwork, iinfo )
770 WRITE( nounit, fmt = 9999 )
'DSYGVD(V,' // uplo //
771 $
')', iinfo, n, jtype, ioldsd
773 IF( iinfo.LT.0 )
THEN
776 result( ntest ) = ulpinv
783 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
784 $ ldz, d, work, result( ntest ) )
790 CALL dlacpy(
' ', n, n, a, lda, ab, lda )
791 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
793 CALL dsygvx( ibtype,
'V',
'A', uplo, n, ab, lda, bb,
794 $ ldb, vl, vu, il, iu, abstol, m, d, z,
795 $ ldz, work, nwork, iwork( n+1 ), iwork,
797 IF( iinfo.NE.0 )
THEN
798 WRITE( nounit, fmt = 9999 )
'DSYGVX(V,A' // uplo //
799 $
')', iinfo, n, jtype, ioldsd
801 IF( iinfo.LT.0 )
THEN
804 result( ntest ) = ulpinv
811 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
812 $ ldz, d, work, result( ntest ) )
816 CALL dlacpy(
' ', n, n, a, lda, ab, lda )
817 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
826 CALL dsygvx( ibtype,
'V',
'V', uplo, n, ab, lda, bb,
827 $ ldb, vl, vu, il, iu, abstol, m, d, z,
828 $ ldz, work, nwork, iwork( n+1 ), iwork,
830 IF( iinfo.NE.0 )
THEN
831 WRITE( nounit, fmt = 9999 )
'DSYGVX(V,V,' //
832 $ uplo //
')', iinfo, n, jtype, ioldsd
834 IF( iinfo.LT.0 )
THEN
837 result( ntest ) = ulpinv
844 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
845 $ ldz, d, work, result( ntest ) )
849 CALL dlacpy(
' ', n, n, a, lda, ab, lda )
850 CALL dlacpy( uplo, n, n, b, ldb, bb, ldb )
852 CALL dsygvx( ibtype,
'V',
'I', uplo, n, ab, lda, bb,
853 $ ldb, vl, vu, il, iu, abstol, m, d, z,
854 $ ldz, work, nwork, iwork
856 IF( iinfo.NE.0 )
THEN
857 WRITE( nounit, fmt = 9999 )
'DSYGVX(V,I,' //
858 $ uplo //
')', iinfo, n, jtype, ioldsd
860 IF( iinfo.LT.0 )
THEN
863 result( ntest ) = ulpinv
870 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
871 $ ldz, d, work, result(
881 IF( lsame( uplo,
'U' ) )
THEN
901 CALL dspgv( ibtype,
'V', uplo, n, ap, bp, d, z, ldz,
903 IF( iinfo.NE.0 )
THEN
904 WRITE( nounit, fmt = 9999 )
'DSPGV(V,' // uplo //
905 $
')', iinfo, n, jtype, ioldsd
907 IF( iinfo.LT.0 )
THEN
910 result( ntest ) = ulpinv
917 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
918 $ ldz, d, work, result( ntest ) )
926 IF( lsame( uplo,
'U' ) )
THEN
946 CALL dspgvd( ibtype,
'V', uplo, n, ap, bp, d, z, ldz,
947 $ work, nwork, iwork, liwork, iinfo )
948 IF( iinfo.NE.0 )
THEN
949 WRITE( nounit, fmt = 9999 )
'DSPGVD(V,' // uplo //
950 $
')', iinfo, n, jtype, ioldsd
952 IF( iinfo.LT.0 )
THEN
955 result( ntest ) = ulpinv
962 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
963 $ ldz, d, work, result( ntest ) )
971 IF( lsame( uplo,
'U' ) )
THEN
991 CALL dspgvx( ibtype,
'V',
'A', uplo, n, ap, bp, vl,
992 $ vu, il, iu, abstol, m, d, z, ldz, work,
993 $ iwork( n+1 ), iwork, info )
994 IF( iinfo.NE.0 )
THEN
995 WRITE( nounit, fmt = 9999 )
'DSPGVX(V,A' // uplo //
996 $
')', iinfo, n, jtype, ioldsd
998 IF( iinfo.LT.0 )
THEN
1001 result( ntest ) = ulpinv
1008 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
1009 $ ldz, d, work, result( ntest ) )
1015 IF( lsame( uplo
'U' ) )
THEN
1019 ap( ij ) = a( i, j )
1020 bp( ij ) = b( i, j )
1028 ap( ij ) = a( i, j )
1029 bp( ij ) = b( i, j )
1037 CALL dspgvx( ibtype,
'V',
'V', uplo, n, ap, bp, vl,
1038 $ vu, il, iu, abstol, m, d, z, ldz, work,
1039 $ iwork( n+1 ), iwork, info )
1040 IF( iinfo.NE.0 )
THEN
1041 WRITE( nounit, fmt = 9999 )
'DSPGVX(V,V' // uplo //
1042 $
')', iinfo, n, jtype, ioldsd
1044 IF( iinfo.LT.0 )
THEN
1047 result( ntest ) = ulpinv
1054 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
1055 $ ldz, d, work, result( ntest ) )
1061 IF( lsame( uplo,
'U' ) )
THEN
1065 ap( ij ) = a( i, j )
1066 bp( ij ) = b( i, j )
1074 ap( ij ) = a( i, j )
1075 bp( ij ) = b( i, j )
1081 CALL dspgvx( ibtype,
'V',
'I', uplo, n, ap, bp, vl,
1082 $ vu, il, iu, abstol, m, d, z, ldz, work,
1083 $ iwork( n+1 ), iwork, info )
1084 IF( iinfo.NE.0 )
THEN
1085 WRITE( nounit, fmt = 9999 )
'DSPGVX(V,I' // uplo //
1086 $
')', iinfo, n, jtype, ioldsd
1088 IF( iinfo.LT.0 )
THEN
1091 result( ntest ) = ulpinv
1098 CALL dsgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
1099 $ ldz, d, work, result( ntest ) )
1103 IF( ibtype.EQ.1 )
THEN
1111 IF( lsame( uplo,
'U' ) )
THEN
1113 DO 320 i =
max( 1, j-ka ), j
1114 ab( ka+1+i-j, j ) = a( i, j )
1116 DO 330 i =
max( 1, j-kb ), j
1117 bb( kb+1+i-j, j ) = b( i, j )
1122 DO 350 i = j,
min( n, j+ka )
1123 ab( 1+i-j, j ) = a( i, j )
1125 DO 360 i = j,
min( n, j+kb )
1126 bb( 1+i-j, j ) = b( i, j )
1131 CALL dsbgv(
'V', uplo, n, ka, kb, ab, lda, bb, ldb,
1132 $ d, z, ldz, work, iinfo )
1133 IF( iinfo.NE.0 )
THEN
1134 WRITE( nounit, fmt = 9999 )
'DSBGV(V,' //
1135 $ uplo //
')', iinfo, n, jtype, ioldsd
1137 IF( iinfo.LT.0 )
THEN
1140 result( ntest ) = ulpinv
1147 CALL dsgt01( ibtype, uplo, n
1148 $ ldz, d, work, result
1156 IF( lsame( uplo,
'U' ) )
THEN
1158 DO 380 i =
max( 1, j-ka ), j
1159 ab( ka+1+i-j, j ) = a( i, j )
1161 DO 390 i =
max( 1, j-kb ), j
1162 bb( kb+1+i-j, j ) = b( i, j )
1167 DO 410 i = j,
min( n, j+ka )
1168 ab( 1+i-j, j ) = a( i, j )
1170 DO 420 i = j,
min( n, j+kb )
1171 bb( 1+i-j, j ) = b( i, j )
1176 CALL dsbgvd(
'V', uplo, n, ka, kb, ab, lda, bb,
1177 $ ldb, d, z, ldz, work, nwork, iwork,
1179 IF( iinfo.NE.0 )
THEN
1180 WRITE( nounit, fmt = 9999 )
'DSBGVD(V,'
1181 $ uplo //
')', iinfo, n, jtype, ioldsd
1183 IF( iinfo.LT.0 )
THEN
1186 result( ntest ) = ulpinv
1193 CALL dsgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
1194 $ ldz, d, work, result( ntest ) )
1202 IF( lsame( uplo, 'u
' ) ) THEN
1204 DO 440 I = MAX( 1, J-KA ), J
1205 AB( KA+1+I-J, J ) = A( I, J )
1207 DO 450 I = MAX( 1, J-KB ), J
1208 BB( KB+1+I-J, J ) = B( I, J )
1213 DO 470 I = J, MIN( N, J+KA )
1214 AB( 1+I-J, J ) = A( I, J )
1216 DO 480 I = J, MIN( N, J+KB )
1217 BB( 1+I-J, J ) = B( I, J )
1222 CALL DSBGVX( 'v
', 'a
', UPLO, N, KA, KB, AB, LDA,
1223 $ BB, LDB, BP, MAX( 1, N ), VL, VU, IL,
1224 $ IU, ABSTOL, M, D, Z, LDZ, WORK,
1225 $ IWORK( N+1 ), IWORK, IINFO )
1226.NE.
IF( IINFO0 ) THEN
1227 WRITE( NOUNIT, FMT = 9999 )'dsbgvx(v,a
' //
1228 $ UPLO // ')
', IINFO, N, JTYPE, IOLDSD
1230.LT.
IF( IINFO0 ) THEN
1233 RESULT( NTEST ) = ULPINV
1240 CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
1241 $ LDZ, D, WORK, RESULT( NTEST ) )
1248 IF( LSAME( UPLO, 'u
' ) ) THEN
1250 DO 500 I = MAX( 1, J-KA ), J
1251 AB( KA+1+I-J, J ) = A( I, J )
1253 DO 510 I = MAX( 1, J-KB ), J
1254 BB( KB+1+I-J, J ) = B( I, J )
1259 DO 530 I = J, MIN( N, J+KA )
1260 AB( 1+I-J, J ) = A( I, J )
1262 DO 540 I = J, MIN( N, J+KB )
1263 BB( 1+I-J, J ) = B( I, J )
1270 CALL DSBGVX( 'v
', 'v
', UPLO, N, KA, KB, AB, LDA,
1271 $ BB, LDB, BP, MAX( 1, N ), VL, VU, IL,
1272 $ IU, ABSTOL, M, D, Z, LDZ, WORK,
1273 $ IWORK( N+1 ), IWORK, IINFO )
1274.NE.
IF( IINFO0 ) THEN
1275 WRITE( NOUNIT, FMT = 9999 )'dsbgvx(v,v
' //
1276 $ UPLO // ')
', IINFO, N, JTYPE, IOLDSD
1278.LT.
IF( IINFO0 ) THEN
1281 RESULT( NTEST ) = ULPINV
1288 CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
1289 $ LDZ, D, WORK, RESULT( NTEST ) )
1295 IF( LSAME( UPLO, 'u
' ) ) THEN
1297 DO 560 I = MAX( 1, J-KA ), J
1298 AB( KA+1+I-J, J ) = A( I, J )
1300 DO 570 I = MAX( 1, J-KB ), J
1301 BB( KB+1+I-J, J ) = B( I, J )
1306 DO 590 I = J, MIN( N, J+KA )
1307 AB( 1+I-J, J ) = A( I, J )
1309 DO 600 I = J, MIN( N, J+KB )
1310 BB( 1+I-J, J ) = B( I, J )
1315 CALL DSBGVX( 'v
', 'i
', UPLO, N, KA, KB, AB, LDA,
1316 $ BB, LDB, BP, MAX( 1, N ), VL, VU, IL,
1317 $ IU, ABSTOL, M, D, Z, LDZ, WORK,
1318 $ IWORK( N+1 ), IWORK, IINFO )
1319.NE.
IF( IINFO0 ) THEN
1320 WRITE( NOUNIT, FMT = 9999 )'dsbgvx(v,i
' //
1321 $ UPLO // ')
', IINFO, N, JTYPE, IOLDSD
1323.LT.
IF( IINFO0 ) THEN
1326 RESULT( NTEST ) = ULPINV
1333 CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
1334 $ LDZ, D, WORK, RESULT( NTEST ) )
1343 NTESTT = NTESTT + NTEST
1344 CALL DLAFTS( 'dsg
', N, N, JTYPE, NTEST, RESULT, IOLDSD,
1345 $ THRESH, NOUNIT, NERRS )
1351 CALL DLASUM( 'dsg
', NOUNIT, NERRS, NTESTT )
1357 9999 FORMAT( ' ddrvsg2stg:
', A, ' returned info=
', I6, '.
', / 9X,
1358 $ 'n=
', I6, ', jtype=
', I6, ', iseed=(
', 3( I5, ',
' ), I5, ')
' )