numeric-linalg

Educational material on the SciPy implementation of numerical linear algebra algorithms

NameSizeMode
..
lapack/TESTING/EIG/ddrvsg2stg.f 49124B -rw-r--r--
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446
0447
0448
0449
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459
0460
0461
0462
0463
0464
0465
0466
0467
0468
0469
0470
0471
0472
0473
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495
0496
0497
0498
0499
0500
0501
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522
0523
0524
0525
0526
0527
0528
0529
0530
0531
0532
0533
0534
0535
0536
0537
0538
0539
0540
0541
0542
0543
0544
0545
0546
0547
0548
0549
0550
0551
0552
0553
0554
0555
0556
0557
0558
0559
0560
0561
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599
0600
0601
0602
0603
0604
0605
0606
0607
0608
0609
0610
0611
0612
0613
0614
0615
0616
0617
0618
0619
0620
0621
0622
0623
0624
0625
0626
0627
0628
0629
0630
0631
0632
0633
0634
0635
0636
0637
0638
0639
0640
0641
0642
0643
0644
0645
0646
0647
0648
0649
0650
0651
0652
0653
0654
0655
0656
0657
0658
0659
0660
0661
0662
0663
0664
0665
0666
0667
0668
0669
0670
0671
0672
0673
0674
0675
0676
0677
0678
0679
0680
0681
0682
0683
0684
0685
0686
0687
0688
0689
0690
0691
0692
0693
0694
0695
0696
0697
0698
0699
0700
0701
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715
0716
0717
0718
0719
0720
0721
0722
0723
0724
0725
0726
0727
0728
0729
0730
0731
0732
0733
0734
0735
0736
0737
0738
0739
0740
0741
0742
0743
0744
0745
0746
0747
0748
0749
0750
0751
0752
0753
0754
0755
0756
0757
0758
0759
0760
0761
0762
0763
0764
0765
0766
0767
0768
0769
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782
0783
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802
0803
0804
0805
0806
0807
0808
0809
0810
0811
0812
0813
0814
0815
0816
0817
0818
0819
0820
0821
0822
0823
0824
0825
0826
0827
0828
0829
0830
0831
0832
0833
0834
0835
0836
0837
0838
0839
0840
0841
0842
0843
0844
0845
0846
0847
0848
0849
0850
0851
0852
0853
0854
0855
0856
0857
0858
0859
0860
0861
0862
0863
0864
0865
0866
0867
0868
0869
0870
0871
0872
0873
0874
0875
0876
0877
0878
0879
0880
0881
0882
0883
0884
0885
0886
0887
0888
0889
0890
0891
0892
0893
0894
0895
0896
0897
0898
0899
0900
0901
0902
0903
0904
0905
0906
0907
0908
0909
0910
0911
0912
0913
0914
0915
0916
0917
0918
0919
0920
0921
0922
0923
0924
0925
0926
0927
0928
0929
0930
0931
0932
0933
0934
0935
0936
0937
0938
0939
0940
0941
0942
0943
0944
0945
0946
0947
0948
0949
0950
0951
0952
0953
0954
0955
0956
0957
0958
0959
0960
0961
0962
0963
0964
0965
0966
0967
0968
0969
0970
0971
0972
0973
0974
0975
0976
0977
0978
0979
0980
0981
0982
0983
0984
0985
0986
0987
0988
0989
0990
0991
0992
0993
0994
0995
0996
0997
0998
0999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
*> \brief \b DDRVSG2STG
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DDRVSG2STG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
*                              NOUNIT, A, LDA, B, LDB, D, D2, Z, LDZ, AB,
*                              BB, AP, BP, WORK, NWORK, IWORK, LIWORK, 
*                              RESULT, INFO )
*
*       IMPLICIT NONE
*       .. Scalar Arguments ..
*       INTEGER            INFO, LDA, LDB, LDZ, LIWORK, NOUNIT, NSIZES,
*      $                   NTYPES, NWORK
*       DOUBLE PRECISION   THRESH
*       ..
*       .. Array Arguments ..
*       LOGICAL            DOTYPE( * )
*       INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
*       DOUBLE PRECISION   A( LDA, * ), AB( LDA, * ), AP( * ),
*      $                   B( LDB, * ), BB( LDB, * ), BP( * ), D( * ),
*      $                   RESULT( * ), WORK( * ), Z( LDZ, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>      DDRVSG2STG checks the real symmetric generalized eigenproblem
*>      drivers.
*>
*>              DSYGV computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem.
*>
*>              DSYGVD computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem using a divide and conquer algorithm.
*>
*>              DSYGVX computes selected eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem.
*>
*>              DSPGV computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem in packed storage.
*>
*>              DSPGVD computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem in packed storage using a divide and
*>              conquer algorithm.
*>
*>              DSPGVX computes selected eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem in packed storage.
*>
*>              DSBGV computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite banded
*>              generalized eigenproblem.
*>
*>              DSBGVD computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite banded
*>              generalized eigenproblem using a divide and conquer
*>              algorithm.
*>
*>              DSBGVX computes selected eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite banded
*>              generalized eigenproblem.
*>
*>      When DDRVSG2STG is called, a number of matrix "sizes" ("n's") and a
*>      number of matrix "types" are specified.  For each size ("n")
*>      and each type of matrix, one matrix A of the given type will be
*>      generated; a random well-conditioned matrix B is also generated
*>      and the pair (A,B) is used to test the drivers.
*>
*>      For each pair (A,B), the following tests are performed:
*>
*>      (1) DSYGV with ITYPE = 1 and UPLO ='U':
*>
*>              | A Z - B Z D | / ( |A| |Z| n ulp )
*>              | D - D2 | / ( |D| ulp )   where D is computed by
*>                                         DSYGV and  D2 is computed by
*>                                         DSYGV_2STAGE. This test is
*>                                         only performed for DSYGV
*>
*>      (2) as (1) but calling DSPGV
*>      (3) as (1) but calling DSBGV
*>      (4) as (1) but with UPLO = 'L'
*>      (5) as (4) but calling DSPGV
*>      (6) as (4) but calling DSBGV
*>
*>      (7) DSYGV with ITYPE = 2 and UPLO ='U':
*>
*>              | A B Z - Z D | / ( |A| |Z| n ulp )
*>
*>      (8) as (7) but calling DSPGV
*>      (9) as (7) but with UPLO = 'L'
*>      (10) as (9) but calling DSPGV
*>
*>      (11) DSYGV with ITYPE = 3 and UPLO ='U':
*>
*>              | B A Z - Z D | / ( |A| |Z| n ulp )
*>
*>      (12) as (11) but calling DSPGV
*>      (13) as (11) but with UPLO = 'L'
*>      (14) as (13) but calling DSPGV
*>
*>      DSYGVD, DSPGVD and DSBGVD performed the same 14 tests.
*>
*>      DSYGVX, DSPGVX and DSBGVX performed the above 14 tests with
*>      the parameter RANGE = 'A', 'N' and 'I', respectively.
*>
*>      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.
*>      This type is used for the matrix A which has half-bandwidth KA.
*>      B is generated as a well-conditioned positive definite matrix
*>      with half-bandwidth KB (<= KA).
*>      Currently, the list of possible types for A 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 orthogonal 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 orthogonal 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 orthogonal 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) symmetric 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)
*>
*>      (16) Same as (8), but with KA = 1 and KB = 1
*>      (17) Same as (8), but with KA = 2 and KB = 1
*>      (18) Same as (8), but with KA = 2 and KB = 2
*>      (19) Same as (8), but with KA = 3 and KB = 1
*>      (20) Same as (8), but with KA = 3 and KB = 2
*>      (21) Same as (8), but with KA = 3 and KB = 3
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \verbatim
*>  NSIZES  INTEGER
*>          The number of sizes of matrices to use.  If it is zero,
*>          DDRVSG2STG does nothing.  It must be at least zero.
*>          Not modified.
*>
*>  NN      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.
*>          Not modified.
*>
*>  NTYPES  INTEGER
*>          The number of elements in DOTYPE.   If it is zero, DDRVSG2STG
*>          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. .
*>          Not modified.
*>
*>  DOTYPE  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.
*>          Not modified.
*>
*>  ISEED   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 DDRVSG2STG to continue the same random number
*>          sequence.
*>          Modified.
*>
*>  THRESH  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.
*>          Not modified.
*>
*>  NOUNIT  INTEGER
*>          The FORTRAN unit number for printing out error messages
*>          (e.g., if a routine returns IINFO not equal to 0.)
*>          Not modified.
*>
*>  A       DOUBLE PRECISION array, dimension (LDA , max(NN))
*>          Used to hold the matrix whose eigenvalues are to be
*>          computed.  On exit, A contains the last matrix actually
*>          used.
*>          Modified.
*>
*>  LDA     INTEGER
*>          The leading dimension of A and AB.  It must be at
*>          least 1 and at least max( NN ).
*>          Not modified.
*>
*>  B       DOUBLE PRECISION array, dimension (LDB , max(NN))
*>          Used to hold the symmetric positive definite matrix for
*>          the generalized problem.
*>          On exit, B contains the last matrix actually
*>          used.
*>          Modified.
*>
*>  LDB     INTEGER
*>          The leading dimension of B and BB.  It must be at
*>          least 1 and at least max( NN ).
*>          Not modified.
*>
*>  D       DOUBLE PRECISION array, dimension (max(NN))
*>          The eigenvalues of A. On exit, the eigenvalues in D
*>          correspond with the matrix in A.
*>          Modified.
*>
*>  Z       DOUBLE PRECISION array, dimension (LDZ, max(NN))
*>          The matrix of eigenvectors.
*>          Modified.
*>
*>  LDZ     INTEGER
*>          The leading dimension of Z.  It must be at least 1 and
*>          at least max( NN ).
*>          Not modified.
*>
*>  AB      DOUBLE PRECISION array, dimension (LDA, max(NN))
*>          Workspace.
*>          Modified.
*>
*>  BB      DOUBLE PRECISION array, dimension (LDB, max(NN))
*>          Workspace.
*>          Modified.
*>
*>  AP      DOUBLE PRECISION array, dimension (max(NN)**2)
*>          Workspace.
*>          Modified.
*>
*>  BP      DOUBLE PRECISION array, dimension (max(NN)**2)
*>          Workspace.
*>          Modified.
*>
*>  WORK    DOUBLE PRECISION array, dimension (NWORK)
*>          Workspace.
*>          Modified.
*>
*>  NWORK   INTEGER
*>          The number of entries in WORK.  This must be at least
*>          1+5*N+2*N*lg(N)+3*N**2 where N = max( NN(j) ) and
*>          lg( N ) = smallest integer k such that 2**k >= N.
*>          Not modified.
*>
*>  IWORK   INTEGER array, dimension (LIWORK)
*>          Workspace.
*>          Modified.
*>
*>  LIWORK  INTEGER
*>          The number of entries in WORK.  This must be at least 6*N.
*>          Not modified.
*>
*>  RESULT  DOUBLE PRECISION array, dimension (70)
*>          The values computed by the 70 tests described above.
*>          Modified.
*>
*>  INFO    INTEGER
*>          If 0, then everything ran OK.
*>           -1: NSIZES < 0
*>           -2: Some NN(j) < 0
*>           -3: NTYPES < 0
*>           -5: THRESH < 0
*>           -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
*>          -16: LDZ < 1 or LDZ < NMAX.
*>          -21: NWORK too small.
*>          -23: LIWORK too small.
*>          If  DLATMR, SLATMS, DSYGV, DSPGV, DSBGV, SSYGVD, SSPGVD,
*>              DSBGVD, DSYGVX, DSPGVX or SSBGVX returns an error code,
*>              the absolute value of it is returned.
*>          Modified.
*>
*> ----------------------------------------------------------------------
*>
*>       Some Local Variables and Parameters:
*>       ---- ----- --------- --- ----------
*>       ZERO, ONE       Real 0 and 1.
*>       MAXTYP          The number of types defined.
*>       NTEST           The number of tests that have been run
*>                       on this matrix.
*>       NTESTT          The total number of tests for this call.
*>       NMAX            Largest value in NN.
*>       NMATS           The number of matrices generated so far.
*>       NERRS           The number of tests which have exceeded THRESH
*>                       so far (computed by DLAFTS).
*>       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) )
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup double_eig
*
*  =====================================================================
      SUBROUTINE DDRVSG2STG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
     $                       NOUNIT, A, LDA, B, LDB, D, D2, Z, LDZ, AB,
     $                       BB, AP, BP, WORK, NWORK, IWORK, LIWORK, 
     $                       RESULT, INFO )
*
      IMPLICIT NONE
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDB, LDZ, LIWORK, NOUNIT, NSIZES,
     $                   NTYPES, NWORK
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            DOTYPE( * )
      INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
      DOUBLE PRECISION   A( LDA, * ), AB( LDA, * ), AP( * ),
     $                   B( LDB, * ), BB( LDB, * ), BP( * ), D( * ),
     $                   D2( * ), RESULT( * ), WORK( * ), Z( LDZ, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TEN
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 )
      INTEGER            MAXTYP
      PARAMETER          ( MAXTYP = 21 )
*     ..
*     .. Local Scalars ..
      LOGICAL            BADNN
      CHARACTER          UPLO
      INTEGER            I, IBTYPE, IBUPLO, IINFO, IJ, IL, IMODE, ITEMP,
     $                   ITYPE, IU, J, JCOL, JSIZE, JTYPE, KA, KA9, KB,
     $                   KB9, M, MTYPES, N, NERRS, NMATS, NMAX, NTEST,
     $                   NTESTT
      DOUBLE PRECISION   ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
     $                   RTUNFL, ULP, ULPINV, UNFL, VL, VU, TEMP1, TEMP2
*     ..
*     .. Local Arrays ..
      INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
     $                   KMAGN( MAXTYP ), KMODE( MAXTYP ),
     $                   KTYPE( MAXTYP )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH, DLARND
      EXTERNAL           LSAME, DLAMCH, DLARND
*     ..
*     .. External Subroutines ..
      EXTERNAL           DLACPY, DLAFTS, DLASET, DLASUM, DLATMR,
     $                   DLATMS, DSBGV, DSBGVD, DSBGVX, DSGT01, DSPGV,
     $                   DSPGVD, DSPGVX, DSYGV, DSYGVD, DSYGVX, XERBLA,
     $                   DSYGV_2STAGE
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, DBLE, MAX, MIN, SQRT
*     ..
*     .. Data statements ..
      DATA               KTYPE / 1, 2, 5*4, 5*5, 3*8, 6*9 /
      DATA               KMAGN / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
     $                   2, 3, 6*1 /
      DATA               KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
     $                   0, 0, 6*4 /
*     ..
*     .. Executable Statements ..
*
*     1)      Check for errors
*
      NTESTT = 0
      INFO = 0
*
      BADNN = .FALSE.
      NMAX = 0
      DO 10 J = 1, NSIZES
         NMAX = MAX( NMAX, NN( J ) )
         IF( NN( J ).LT.0 )
     $      BADNN = .TRUE.
   10 CONTINUE
*
*     Check for errors
*
      IF( NSIZES.LT.0 ) THEN
         INFO = -1
      ELSE IF( BADNN ) THEN
         INFO = -2
      ELSE IF( NTYPES.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
         INFO = -9
      ELSE IF( LDZ.LE.1 .OR. LDZ.LT.NMAX ) THEN
         INFO = -16
      ELSE IF( 2*MAX( NMAX, 3 )**2.GT.NWORK ) THEN
         INFO = -21
      ELSE IF( 2*MAX( NMAX, 3 )**2.GT.LIWORK ) THEN
         INFO = -23
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DDRVSG2STG', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
     $   RETURN
*
*     More Important constants
*
      UNFL = DLAMCH( 'Safe minimum' )
      OVFL = DLAMCH( 'Overflow' )
      ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
      ULPINV = ONE / ULP
      RTUNFL = SQRT( UNFL )
      RTOVFL = SQRT( OVFL )
*
      DO 20 I = 1, 4
         ISEED2( I ) = ISEED( I )
   20 CONTINUE
*
*     Loop over sizes, types
*
      NERRS = 0
      NMATS = 0
*
      DO 650 JSIZE = 1, NSIZES
         N = NN( JSIZE )
         ANINV = ONE / DBLE( MAX( 1, N ) )
*
         IF( NSIZES.NE.1 ) THEN
            MTYPES = MIN( MAXTYP, NTYPES )
         ELSE
            MTYPES = MIN( MAXTYP+1, NTYPES )
         END IF
*
         KA9 = 0
         KB9 = 0
         DO 640 JTYPE = 1, MTYPES
            IF( .NOT.DOTYPE( JTYPE ) )
     $         GO TO 640
            NMATS = NMATS + 1
            NTEST = 0
*
            DO 30 J = 1, 4
               IOLDSD( J ) = ISEED( J )
   30       CONTINUE
*
*           2)      Compute "A"
*
*                   Control parameters:
*
*               KMAGN  KMODE        KTYPE
*           =1  O(1)   clustered 1  zero
*           =2  large  clustered 2  identity
*           =3  small  exponential  (none)
*           =4         arithmetic   diagonal, w/ eigenvalues
*           =5         random log   hermitian, w/ eigenvalues
*           =6         random       (none)
*           =7                      random diagonal
*           =8                      random hermitian
*           =9                      banded, w/ eigenvalues
*
            IF( MTYPES.GT.MAXTYP )
     $         GO TO 90
*
            ITYPE = KTYPE( JTYPE )
            IMODE = KMODE( JTYPE )
*
*           Compute norm
*
            GO TO ( 40, 50, 60 )KMAGN( JTYPE )
*
   40       CONTINUE
            ANORM = ONE
            GO TO 70
*
   50       CONTINUE
            ANORM = ( RTOVFL*ULP )*ANINV
            GO TO 70
*
   60       CONTINUE
            ANORM = RTUNFL*N*ULPINV
            GO TO 70
*
   70       CONTINUE
*
            IINFO = 0
            COND = ULPINV
*
*           Special Matrices -- Identity & Jordan block
*
            IF( ITYPE.EQ.1 ) THEN
*
*              Zero
*
               KA = 0
               KB = 0
               CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
*
            ELSE IF( ITYPE.EQ.2 ) THEN
*
*              Identity
*
               KA = 0
               KB = 0
               CALL DLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
               DO 80 JCOL = 1, N
                  A( JCOL, JCOL ) = ANORM
   80          CONTINUE
*
            ELSE IF( ITYPE.EQ.4 ) THEN
*
*              Diagonal Matrix, [Eigen]values Specified
*
               KA = 0
               KB = 0
               CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
     $                      ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
     $                      IINFO )
*
            ELSE IF( ITYPE.EQ.5 ) THEN
*
*              symmetric, eigenvalues specified
*
               KA = MAX( 0, N-1 )
               KB = KA
               CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
     $                      ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
     $                      IINFO )
*
            ELSE IF( ITYPE.EQ.7 ) THEN
*
*              Diagonal, random eigenvalues
*
               KA = 0
               KB = 0
               CALL DLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
     $                      'T', 'N', WORK( N+1 ), 1, ONE,
     $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
     $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
*
            ELSE IF( ITYPE.EQ.8 ) THEN
*
*              symmetric, random eigenvalues
*
               KA = MAX( 0, N-1 )
               KB = KA
               CALL DLATMR( N, N, 'S', ISEED, 'H', WORK, 6, ONE, ONE,
     $                      'T', 'N', WORK( N+1 ), 1, ONE,
     $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
     $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
*
            ELSE IF( ITYPE.EQ.9 ) THEN
*
*              symmetric banded, eigenvalues specified
*
*              The following values are used for the half-bandwidths:
*
*                ka = 1   kb = 1
*                ka = 2   kb = 1
*                ka = 2   kb = 2
*                ka = 3   kb = 1
*                ka = 3   kb = 2
*                ka = 3   kb = 3
*
               KB9 = KB9 + 1
               IF( KB9.GT.KA9 ) THEN
                  KA9 = KA9 + 1
                  KB9 = 1
               END IF
               KA = MAX( 0, MIN( N-1, KA9 ) )
               KB = MAX( 0, MIN( N-1, KB9 ) )
               CALL DLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
     $                      ANORM, KA, KA, 'N', A, LDA, WORK( N+1 ),
     $                      IINFO )
*
            ELSE
*
               IINFO = 1
            END IF
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               RETURN
            END IF
*
   90       CONTINUE
*
            ABSTOL = UNFL + UNFL
            IF( N.LE.1 ) THEN
               IL = 1
               IU = N
            ELSE
               IL = 1 + INT( ( N-1 )*DLARND( 1, ISEED2 ) )
               IU = 1 + INT( ( N-1 )*DLARND( 1, ISEED2 ) )
               IF( IL.GT.IU ) THEN
                  ITEMP = IL
                  IL = IU
                  IU = ITEMP
               END IF
            END IF
*
*           3) Call DSYGV, DSPGV, DSBGV, SSYGVD, SSPGVD, SSBGVD,
*              DSYGVX, DSPGVX, and DSBGVX, do tests.
*
*           loop over the three generalized problems
*                 IBTYPE = 1: A*x = (lambda)*B*x
*                 IBTYPE = 2: A*B*x = (lambda)*x
*                 IBTYPE = 3: B*A*x = (lambda)*x
*
            DO 630 IBTYPE = 1, 3
*
*              loop over the setting UPLO
*
               DO 620 IBUPLO = 1, 2
                  IF( IBUPLO.EQ.1 )
     $               UPLO = 'U'
                  IF( IBUPLO.EQ.2 )
     $               UPLO = 'L'
*
*                 Generate random well-conditioned positive definite
*                 matrix B, of bandwidth not greater than that of A.
*
                  CALL DLATMS( N, N, 'U', ISEED, 'P', WORK, 5, TEN, ONE,
     $                         KB, KB, UPLO, B, LDB, WORK( N+1 ),
     $                         IINFO )
*
*                 Test DSYGV
*
                  NTEST = NTEST + 1
*
                  CALL DLACPY( ' ', N, N, A, LDA, Z, LDZ )
                  CALL DLACPY( UPLO, N, N, B, LDB, BB, LDB )
*
                  CALL DSYGV( IBTYPE, 'V', UPLO, N, Z, LDZ, BB, LDB, D,
     $                        WORK, NWORK, IINFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )'DSYGV(V,' // UPLO //
     $                  ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL DSGT01( IBTYPE, UPLO, N, N, A, LDA, B, LDB, Z,
     $                         LDZ, D, WORK, RESULT( NTEST ) )
*
*                 Test DSYGV_2STAGE
*
                  NTEST = NTEST + 1
*
                  CALL DLACPY( ' ', N, N, A, LDA, Z, LDZ )
                  CALL DLACPY( UPLO, N, N, B, LDB, BB, LDB )
*
                  CALL DSYGV_2STAGE( IBTYPE, 'N', UPLO, N, Z, LDZ,
     $                               BB, LDB, D2, WORK, NWORK, IINFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )
     $                  'DSYGV_2STAGE(V,' // UPLO //
     $                  ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
C                  CALL DSGT01( IBTYPE, UPLO, N, N, A, LDA, B, LDB, Z,
C     $                         LDZ, D, WORK, RESULT( NTEST ) )
*           
*                 Do Tests | D1 - D2 | / ( |D1| ulp )
*                 D1 computed using the standard 1-stage reduction as reference
*                 D2 computed using the 2-stage reduction
*           
                  TEMP1 = ZERO
                  TEMP2 = ZERO
                  DO 151 J = 1, N
                     TEMP1 = MAX( TEMP1, ABS( D( J ) ), 
     $                                   ABS( D2( J ) ) )
                     TEMP2 = MAX( TEMP2, ABS( D( J )-D2( J ) ) )
  151             CONTINUE
*           
                  RESULT( NTEST ) = TEMP2 / 
     $                              MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
*
*                 Test DSYGVD
*
                  NTEST = NTEST + 1
*
                  CALL DLACPY( ' ', N, N, A, LDA, Z, LDZ )
                  CALL DLACPY( UPLO, N, N, B, LDB, BB, LDB )
*
                  CALL DSYGVD( IBTYPE, 'V', UPLO, N, Z, LDZ, BB, LDB, D,
     $                         WORK, NWORK, IWORK, LIWORK, IINFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )'DSYGVD(V,' // UPLO //
     $                  ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL DSGT01( IBTYPE, UPLO, N, N, A, LDA, B, LDB, Z,
     $                         LDZ, D, WORK, RESULT( NTEST ) )
*
*                 Test DSYGVX
*
                  NTEST = NTEST + 1
*
                  CALL DLACPY( ' ', N, N, A, LDA, AB, LDA )
                  CALL DLACPY( UPLO, N, N, B, LDB, BB, LDB )
*
                  CALL DSYGVX( IBTYPE, 'V', 'A', UPLO, N, AB, LDA, BB,
     $                         LDB, VL, VU, IL, IU, ABSTOL, M, D, Z,
     $                         LDZ, WORK, NWORK, IWORK( N+1 ), IWORK,
     $                         IINFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )'DSYGVX(V,A' // UPLO //
     $                  ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL DSGT01( IBTYPE, UPLO, N, N, A, LDA, B, LDB, Z,
     $                         LDZ, D, WORK, RESULT( NTEST ) )
*
                  NTEST = NTEST + 1
*
                  CALL DLACPY( ' ', N, N, A, LDA, AB, LDA )
                  CALL DLACPY( UPLO, N, N, B, LDB, BB, LDB )
*
*                 since we do not know the exact eigenvalues of this
*                 eigenpair, we just set VL and VU as constants.
*                 It is quite possible that there are no eigenvalues
*                 in this interval.
*
                  VL = ZERO
                  VU = ANORM
                  CALL DSYGVX( IBTYPE, 'V', 'V', UPLO, N, AB, LDA, BB,
     $                         LDB, VL, VU, IL, IU, ABSTOL, M, D, Z,
     $                         LDZ, WORK, NWORK, IWORK( N+1 ), IWORK,
     $                         IINFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )'DSYGVX(V,V,' //
     $                  UPLO // ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
     $                         LDZ, D, WORK, RESULT( NTEST ) )
*
                  NTEST = NTEST + 1
*
                  CALL DLACPY( ' ', N, N, A, LDA, AB, LDA )
                  CALL DLACPY( UPLO, N, N, B, LDB, BB, LDB )
*
                  CALL DSYGVX( IBTYPE, 'V', 'I', UPLO, N, AB, LDA, BB,
     $                         LDB, VL, VU, IL, IU, ABSTOL, M, D, Z,
     $                         LDZ, WORK, NWORK, IWORK( N+1 ), IWORK,
     $                         IINFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )'DSYGVX(V,I,' //
     $                  UPLO // ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
     $                         LDZ, D, WORK, RESULT( NTEST ) )
*
  100             CONTINUE
*
*                 Test DSPGV
*
                  NTEST = NTEST + 1
*
*                 Copy the matrices into packed storage.
*
                  IF( LSAME( UPLO, 'U' ) ) THEN
                     IJ = 1
                     DO 120 J = 1, N
                        DO 110 I = 1, J
                           AP( IJ ) = A( I, J )
                           BP( IJ ) = B( I, J )
                           IJ = IJ + 1
  110                   CONTINUE
  120                CONTINUE
                  ELSE
                     IJ = 1
                     DO 140 J = 1, N
                        DO 130 I = J, N
                           AP( IJ ) = A( I, J )
                           BP( IJ ) = B( I, J )
                           IJ = IJ + 1
  130                   CONTINUE
  140                CONTINUE
                  END IF
*
                  CALL DSPGV( IBTYPE, 'V', UPLO, N, AP, BP, D, Z, LDZ,
     $                        WORK, IINFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )'DSPGV(V,' // UPLO //
     $                  ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 310
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL DSGT01( IBTYPE, UPLO, N, N, A, LDA, B, LDB, Z,
     $                         LDZ, D, WORK, RESULT( NTEST ) )
*
*                 Test DSPGVD
*
                  NTEST = NTEST + 1
*
*                 Copy the matrices into packed storage.
*
                  IF( LSAME( UPLO, 'U' ) ) THEN
                     IJ = 1
                     DO 160 J = 1, N
                        DO 150 I = 1, J
                           AP( IJ ) = A( I, J )
                           BP( IJ ) = B( I, J )
                           IJ = IJ + 1
  150                   CONTINUE
  160                CONTINUE
                  ELSE
                     IJ = 1
                     DO 180 J = 1, N
                        DO 170 I = J, N
                           AP( IJ ) = A( I, J )
                           BP( IJ ) = B( I, J )
                           IJ = IJ + 1
  170                   CONTINUE
  180                CONTINUE
                  END IF
*
                  CALL DSPGVD( IBTYPE, 'V', UPLO, N, AP, BP, D, Z, LDZ,
     $                         WORK, NWORK, IWORK, LIWORK, IINFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )'DSPGVD(V,' // UPLO //
     $                  ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 310
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL DSGT01( IBTYPE, UPLO, N, N, A, LDA, B, LDB, Z,
     $                         LDZ, D, WORK, RESULT( NTEST ) )
*
*                 Test DSPGVX
*
                  NTEST = NTEST + 1
*
*                 Copy the matrices into packed storage.
*
                  IF( LSAME( UPLO, 'U' ) ) THEN
                     IJ = 1
                     DO 200 J = 1, N
                        DO 190 I = 1, J
                           AP( IJ ) = A( I, J )
                           BP( IJ ) = B( I, J )
                           IJ = IJ + 1
  190                   CONTINUE
  200                CONTINUE
                  ELSE
                     IJ = 1
                     DO 220 J = 1, N
                        DO 210 I = J, N
                           AP( IJ ) = A( I, J )
                           BP( IJ ) = B( I, J )
                           IJ = IJ + 1
  210                   CONTINUE
  220                CONTINUE
                  END IF
*
                  CALL DSPGVX( IBTYPE, 'V', 'A', UPLO, N, AP, BP, VL,
     $                         VU, IL, IU, ABSTOL, M, D, Z, LDZ, WORK,
     $                         IWORK( N+1 ), IWORK, INFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )'DSPGVX(V,A' // UPLO //
     $                  ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 310
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
     $                         LDZ, D, WORK, RESULT( NTEST ) )
*
                  NTEST = NTEST + 1
*
*                 Copy the matrices into packed storage.
*
                  IF( LSAME( UPLO, 'U' ) ) THEN
                     IJ = 1
                     DO 240 J = 1, N
                        DO 230 I = 1, J
                           AP( IJ ) = A( I, J )
                           BP( IJ ) = B( I, J )
                           IJ = IJ + 1
  230                   CONTINUE
  240                CONTINUE
                  ELSE
                     IJ = 1
                     DO 260 J = 1, N
                        DO 250 I = J, N
                           AP( IJ ) = A( I, J )
                           BP( IJ ) = B( I, J )
                           IJ = IJ + 1
  250                   CONTINUE
  260                CONTINUE
                  END IF
*
                  VL = ZERO
                  VU = ANORM
                  CALL DSPGVX( IBTYPE, 'V', 'V', UPLO, N, AP, BP, VL,
     $                         VU, IL, IU, ABSTOL, M, D, Z, LDZ, WORK,
     $                         IWORK( N+1 ), IWORK, INFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )'DSPGVX(V,V' // UPLO //
     $                  ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 310
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
     $                         LDZ, D, WORK, RESULT( NTEST ) )
*
                  NTEST = NTEST + 1
*
*                 Copy the matrices into packed storage.
*
                  IF( LSAME( UPLO, 'U' ) ) THEN
                     IJ = 1
                     DO 280 J = 1, N
                        DO 270 I = 1, J
                           AP( IJ ) = A( I, J )
                           BP( IJ ) = B( I, J )
                           IJ = IJ + 1
  270                   CONTINUE
  280                CONTINUE
                  ELSE
                     IJ = 1
                     DO 300 J = 1, N
                        DO 290 I = J, N
                           AP( IJ ) = A( I, J )
                           BP( IJ ) = B( I, J )
                           IJ = IJ + 1
  290                   CONTINUE
  300                CONTINUE
                  END IF
*
                  CALL DSPGVX( IBTYPE, 'V', 'I', UPLO, N, AP, BP, VL,
     $                         VU, IL, IU, ABSTOL, M, D, Z, LDZ, WORK,
     $                         IWORK( N+1 ), IWORK, INFO )
                  IF( IINFO.NE.0 ) THEN
                     WRITE( NOUNIT, FMT = 9999 )'DSPGVX(V,I' // UPLO //
     $                  ')', IINFO, N, JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     IF( IINFO.LT.0 ) THEN
                        RETURN
                     ELSE
                        RESULT( NTEST ) = ULPINV
                        GO TO 310
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
     $                         LDZ, D, WORK, RESULT( NTEST ) )
*
  310             CONTINUE
*
                  IF( IBTYPE.EQ.1 ) THEN
*
*                    TEST DSBGV
*
                     NTEST = NTEST + 1
*
*                    Copy the matrices into band storage.
*
                     IF( LSAME( UPLO, 'U' ) ) THEN
                        DO 340 J = 1, N
                           DO 320 I = MAX( 1, J-KA ), J
                              AB( KA+1+I-J, J ) = A( I, J )
  320                      CONTINUE
                           DO 330 I = MAX( 1, J-KB ), J
                              BB( KB+1+I-J, J ) = B( I, J )
  330                      CONTINUE
  340                   CONTINUE
                     ELSE
                        DO 370 J = 1, N
                           DO 350 I = J, MIN( N, J+KA )
                              AB( 1+I-J, J ) = A( I, J )
  350                      CONTINUE
                           DO 360 I = J, MIN( N, J+KB )
                              BB( 1+I-J, J ) = B( I, J )
  360                      CONTINUE
  370                   CONTINUE
                     END IF
*
                     CALL DSBGV( 'V', UPLO, N, KA, KB, AB, LDA, BB, LDB,
     $                           D, Z, LDZ, WORK, IINFO )
                     IF( IINFO.NE.0 ) THEN
                        WRITE( NOUNIT, FMT = 9999 )'DSBGV(V,' //
     $                     UPLO // ')', IINFO, N, JTYPE, IOLDSD
                        INFO = ABS( IINFO )
                        IF( IINFO.LT.0 ) THEN
                           RETURN
                        ELSE
                           RESULT( NTEST ) = ULPINV
                           GO TO 620
                        END IF
                     END IF
*
*                    Do Test
*
                     CALL DSGT01( IBTYPE, UPLO, N, N, A, LDA, B, LDB, Z,
     $                            LDZ, D, WORK, RESULT( NTEST ) )
*
*                    TEST DSBGVD
*
                     NTEST = NTEST + 1
*
*                    Copy the matrices into band storage.
*
                     IF( LSAME( UPLO, 'U' ) ) THEN
                        DO 400 J = 1, N
                           DO 380 I = MAX( 1, J-KA ), J
                              AB( KA+1+I-J, J ) = A( I, J )
  380                      CONTINUE
                           DO 390 I = MAX( 1, J-KB ), J
                              BB( KB+1+I-J, J ) = B( I, J )
  390                      CONTINUE
  400                   CONTINUE
                     ELSE
                        DO 430 J = 1, N
                           DO 410 I = J, MIN( N, J+KA )
                              AB( 1+I-J, J ) = A( I, J )
  410                      CONTINUE
                           DO 420 I = J, MIN( N, J+KB )
                              BB( 1+I-J, J ) = B( I, J )
  420                      CONTINUE
  430                   CONTINUE
                     END IF
*
                     CALL DSBGVD( 'V', UPLO, N, KA, KB, AB, LDA, BB,
     $                            LDB, D, Z, LDZ, WORK, NWORK, IWORK,
     $                            LIWORK, IINFO )
                     IF( IINFO.NE.0 ) THEN
                        WRITE( NOUNIT, FMT = 9999 )'DSBGVD(V,' //
     $                     UPLO // ')', IINFO, N, JTYPE, IOLDSD
                        INFO = ABS( IINFO )
                        IF( IINFO.LT.0 ) THEN
                           RETURN
                        ELSE
                           RESULT( NTEST ) = ULPINV
                           GO TO 620
                        END IF
                     END IF
*
*                    Do Test
*
                     CALL DSGT01( IBTYPE, UPLO, N, N, A, LDA, B, LDB, Z,
     $                            LDZ, D, WORK, RESULT( NTEST ) )
*
*                    Test DSBGVX
*
                     NTEST = NTEST + 1
*
*                    Copy the matrices into band storage.
*
                     IF( LSAME( UPLO, 'U' ) ) THEN
                        DO 460 J = 1, N
                           DO 440 I = MAX( 1, J-KA ), J
                              AB( KA+1+I-J, J ) = A( I, J )
  440                      CONTINUE
                           DO 450 I = MAX( 1, J-KB ), J
                              BB( KB+1+I-J, J ) = B( I, J )
  450                      CONTINUE
  460                   CONTINUE
                     ELSE
                        DO 490 J = 1, N
                           DO 470 I = J, MIN( N, J+KA )
                              AB( 1+I-J, J ) = A( I, J )
  470                      CONTINUE
                           DO 480 I = J, MIN( N, J+KB )
                              BB( 1+I-J, J ) = B( I, J )
  480                      CONTINUE
  490                   CONTINUE
                     END IF
*
                     CALL DSBGVX( 'V', 'A', UPLO, N, KA, KB, AB, LDA,
     $                            BB, LDB, BP, MAX( 1, N ), VL, VU, IL,
     $                            IU, ABSTOL, M, D, Z, LDZ, WORK,
     $                            IWORK( N+1 ), IWORK, IINFO )
                     IF( IINFO.NE.0 ) THEN
                        WRITE( NOUNIT, FMT = 9999 )'DSBGVX(V,A' //
     $                     UPLO // ')', IINFO, N, JTYPE, IOLDSD
                        INFO = ABS( IINFO )
                        IF( IINFO.LT.0 ) THEN
                           RETURN
                        ELSE
                           RESULT( NTEST ) = ULPINV
                           GO TO 620
                        END IF
                     END IF
*
*                    Do Test
*
                     CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
     $                            LDZ, D, WORK, RESULT( NTEST ) )
*
*
                     NTEST = NTEST + 1
*
*                    Copy the matrices into band storage.
*
                     IF( LSAME( UPLO, 'U' ) ) THEN
                        DO 520 J = 1, N
                           DO 500 I = MAX( 1, J-KA ), J
                              AB( KA+1+I-J, J ) = A( I, J )
  500                      CONTINUE
                           DO 510 I = MAX( 1, J-KB ), J
                              BB( KB+1+I-J, J ) = B( I, J )
  510                      CONTINUE
  520                   CONTINUE
                     ELSE
                        DO 550 J = 1, N
                           DO 530 I = J, MIN( N, J+KA )
                              AB( 1+I-J, J ) = A( I, J )
  530                      CONTINUE
                           DO 540 I = J, MIN( N, J+KB )
                              BB( 1+I-J, J ) = B( I, J )
  540                      CONTINUE
  550                   CONTINUE
                     END IF
*
                     VL = ZERO
                     VU = ANORM
                     CALL DSBGVX( 'V', 'V', UPLO, N, KA, KB, AB, LDA,
     $                            BB, LDB, BP, MAX( 1, N ), VL, VU, IL,
     $                            IU, ABSTOL, M, D, Z, LDZ, WORK,
     $                            IWORK( N+1 ), IWORK, IINFO )
                     IF( IINFO.NE.0 ) THEN
                        WRITE( NOUNIT, FMT = 9999 )'DSBGVX(V,V' //
     $                     UPLO // ')', IINFO, N, JTYPE, IOLDSD
                        INFO = ABS( IINFO )
                        IF( IINFO.LT.0 ) THEN
                           RETURN
                        ELSE
                           RESULT( NTEST ) = ULPINV
                           GO TO 620
                        END IF
                     END IF
*
*                    Do Test
*
                     CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
     $                            LDZ, D, WORK, RESULT( NTEST ) )
*
                     NTEST = NTEST + 1
*
*                    Copy the matrices into band storage.
*
                     IF( LSAME( UPLO, 'U' ) ) THEN
                        DO 580 J = 1, N
                           DO 560 I = MAX( 1, J-KA ), J
                              AB( KA+1+I-J, J ) = A( I, J )
  560                      CONTINUE
                           DO 570 I = MAX( 1, J-KB ), J
                              BB( KB+1+I-J, J ) = B( I, J )
  570                      CONTINUE
  580                   CONTINUE
                     ELSE
                        DO 610 J = 1, N
                           DO 590 I = J, MIN( N, J+KA )
                              AB( 1+I-J, J ) = A( I, J )
  590                      CONTINUE
                           DO 600 I = J, MIN( N, J+KB )
                              BB( 1+I-J, J ) = B( I, J )
  600                      CONTINUE
  610                   CONTINUE
                     END IF
*
                     CALL DSBGVX( 'V', 'I', UPLO, N, KA, KB, AB, LDA,
     $                            BB, LDB, BP, MAX( 1, N ), VL, VU, IL,
     $                            IU, ABSTOL, M, D, Z, LDZ, WORK,
     $                            IWORK( N+1 ), IWORK, IINFO )
                     IF( IINFO.NE.0 ) THEN
                        WRITE( NOUNIT, FMT = 9999 )'DSBGVX(V,I' //
     $                     UPLO // ')', IINFO, N, JTYPE, IOLDSD
                        INFO = ABS( IINFO )
                        IF( IINFO.LT.0 ) THEN
                           RETURN
                        ELSE
                           RESULT( NTEST ) = ULPINV
                           GO TO 620
                        END IF
                     END IF
*
*                    Do Test
*
                     CALL DSGT01( IBTYPE, UPLO, N, M, A, LDA, B, LDB, Z,
     $                            LDZ, D, WORK, RESULT( NTEST ) )
*
                  END IF
*
  620          CONTINUE
  630       CONTINUE
*
*           End of Loop -- Check for RESULT(j) > THRESH
*
            NTESTT = NTESTT + NTEST
            CALL DLAFTS( 'DSG', N, N, JTYPE, NTEST, RESULT, IOLDSD,
     $                   THRESH, NOUNIT, NERRS )
  640    CONTINUE
  650 CONTINUE
*
*     Summary
*
      CALL DLASUM( 'DSG', NOUNIT, NERRS, NTESTT )
*
      RETURN
*
*     End of DDRVSG2STG
*
 9999 FORMAT( ' DDRVSG2STG: ', A, ' returned INFO=', I6, '.', / 9X,
     $    'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
      END