1
        2
        3
        4
        5
        6
        7
        8
        9
       10
       11
       12
       13
       14
       15
       16
       17
       18
       19
       20
       21
       22
       23
       24
       25
       26
       27
       28
       29
       30
       31
       32
       33
       34
       35
       36
       37
       38
       39
       40
       41
       42
       43
       44
       45
       46
       47
       48
       49
       50
       51
       52
       53
       54
       55
       56
       57
       58
       59
       60
       61
       62
       63
       64
       65
       66
       67
       68
       69
       70
       71
       72
       73
       74
       75
       76
       77
       78
       79
       80
       81
       82
       83
       84
       85
       86
       87
       88
       89
       90
       91
       92
       93
       94
       95
       96
       97
       98
       99
      100
      101
      102
      103
      104
      105
      106
      107
      108
      109
      110
      111
      112
      113
      114
      115
      116
      117
      118
      119
      120
      121
      122
      123
      124
      125
      126
      127
      128
      129
      130
      131
      132
      133
      134
      135
      136
      137
      138
      139
      140
      141
      142
      143
      144
      145
      146
      147
      148
      149
      150
      151
      152
      153
      154
      155
      156
      157
      158
      159
      160
      161
      162
      163
      164
      165
      166
      167
      168
      169
      170
      171
      172
      173
      174
      175
      176
      177
      178
      179
      180
      181
      182
      183
      184
      185
      186
      187
      188
      189
      190
      191
      192
      193
      194
      195
      196
      197
      198
      199
      200
      201
      202
      203
      204
      205
      206
      207
      208
      209
      210
      211
      212
      213
      214
      215
      216
      217
      218
      219
      220
      221
      222
      223
      224
      225
      226
      227
      228
      229
      230
      231
      232
      233
      234
      235
      236
      237
      238
      239
      240
      241
      242
      243
      244
      245
      246
      247
      248
      249
      250
      251
      252
      253
      254
      255
      256
      257
      258
      259
      260
      261
      262
      263
      264
      265
      266
      267
      268
      269
      270
      271
      272
      273
      274
      275
      276
      277
      278
      279
      280
      281
      282
      283
      284
      285
      286
      287
      288
      289
      290
      291
      292
      293
      294
      295
      296
      297
      298
      299
      300
      301
      302
      303
      304
      305
      306
      307
      308
      309
      310
      311
      312
      313
      314
      315
      316
      317
      318
      319
      320
      321
      322
      323
      324
      325
      326
      327
      328
      329
      330
      331
      332
      333
      334
      335
      336
      337
      338
      339
      340
      341
      342
      343
      344
      345
      346
      347
      348
      349
      350
      351
      352
      353
      354
      355
      356
      357
      358
      359
      360
      361
      362
      363
      364
      365
      366
      367
      368
      369
      370
      371
      372
      373
      374
      375
      376
      377
      378
      379
      380
      381
      382
      383
      384
      385
      386
      387
      388
      389
      390
      391
      392
      393
      394
      395
      396
      397
      398
      399
      400
      401
      402
      403
      404
      405
      406
      407
      408
      409
      410
      411
      412
      413
      414
      415
      416
      417
      418
      419
      420
      421
      422
      423
      424
      425
      426
      427
      428
      429
      430
      431
      432
      433
      434
      435
      436
      437
      438
      439
      440
      441
      442
      443
      444
      445
      446
      447
      448
      449
      450
      451
      452
      453
      454
      455
      456
      457
      458
      459
      460
      461
      462
      463
      464
      465
      466
      467
      468
      469
      470
      471
      472
      473
      474
      475
      476
      477
      478
      479
      480
      481
      482
      483
      484
      485
      486
      487
      488
      489
      490
      491
      492
      493
      494
      495
      496
      497
      498
      499
      500
      501
      502
      503
      504
      505
      506
      507
      508
      509
      510
      511
      512
      513
      514
      515
      516
      517
      518
      519
      520
      521
      522
      523
      524
      525
      526
      527
      528
      529
      530
      531
      532
      533
      534
      535
      536
      537
      538
      539
      540
      541
      542
      543
      544
      545
      546
      547
      548
      549
      550
      551
      552
      553
      554
      555
      556
      557
      558
      559
      560
      561
      562
      563
      564
      565
      566
      567
      568
      569
      570
      571
      572
      573
      574
      575
      576
      577
      578
      579
      580
      581
      582
      583
      584
      585
      586
      587
      588
      589
      590
      591
      592
      593
      594
      595
      596
      597
      598
      599
      600
      601
      602
      603
      604
      605
      606
      607
      608
      609
      610
      611
      612
      613
      614
      615
      616
      617
      618
      619
      620
      621
      622
      623
      624
      625
      626
      627
      628
      629
      630
      631
      632
      633
      634
      635
      636
      637
      638
      639
      640
      641
      642
      643
      644
      645
      646
      647
      648
      649
      650
      651
      652
      653
      654
      655
      656
      657
      658
      659
      660
      661
      662
      663
      664
      665
      666
      667
      668
      669
      670
      671
      672
      673
      674
      675
      676
      677
      678
      679
      680
      681
      682
      683
      684
      685
      686
      687
      688
      689
      690
      691
      692
      693
      694
      695
      696
      697
      698
      699
      700
      701
      702
      703
      704
      705
      706
      707
      708
      709
      710
      711
      712
      713
      714
      715
      716
      717
      718
      719
      720
      721
      722
      723
      724
      725
      726
      727
      728
      729
      730
      731
      732
      733
      734
      735
      736
      737
      738
      739
      740
      741
      742
      743
      744
      745
      746
      747
      748
      749
      750
      751
      752
      753
      754
      755
      756
      757
      758
      759
      760
      761
      762
      763
      764
      765
      766
      767
      768
      769
      770
      771
      772
      773
      774
      775
      776
      777
      778
      779
      780
      781
      782
      783
      784
      785
      786
      787
      788
      789
      790
      791
      792
      793
      794
      795
      796
      797
      798
      799
      800
      801
      802
      803
      804
      805
      806
      807
      808
      809
      810
      811
      812
      813
      814
      815
      816
      817
      818
      819
      820
      821
      822
      823
      824
      825
      826
      827
      828
      829
      830
      831
      832
      833
      834
      835
      836
      837
      838
      839
      840
      841
      842
      843
      844
      845
      846
      847
      848
      849
      850
      851
      852
      853
      854
      855
      856
      857
      858
      859
      860
      861
      862
      863
      864
      865
      866
      867
      868
      869
      870
      871
      872
      873
      874
      875
      876
      877
      878
      879
      880
      881
      882
      883
      884
      885
      886
      887
      888
      889
      890
      891
      892
      893
      894
      895
      896
      897
      898
      899
      900
      901
      902
      903
      904
      905
      906
      907
      908
      909
      910
      911
      912
      913
      914
      915
      916
      917
      918
      919
      920
      921
      922
      923
      924
      925
      926
      927
      928
      929
      930
      931
      932
      933
      934
      935
      936
      937
      938
      939
      940
      941
      942
      943
      944
      945
      946
      947
      948
      949
      950
      951
      952
      953
      954
      955
      956
      957
      958
      959
      960
      961
      962
      963
      964
      965
      966
      967
      968
      969
      970
      971
      972
      973
      974
      975
      976
      977
      978
      979
      980
      981
      982
      983
      984
      985
      986
      987
      988
      989
      990
      991
      992
      993
      994
      995
      996
      997
      998
      999
     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
      SUBROUTINE SSBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
     $                   LDX, WORK, INFO )
*
*  -- LAPACK routine (version 3.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2006
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO, VECT
      INTEGER            INFO, KA, KB, LDAB, LDBB, LDX, N
*     ..
*     .. Array Arguments ..
      REAL               AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
     $                   X( LDX, * )
*     ..
*
*  Purpose
*  =======
*
*  SSBGST reduces a real symmetric-definite banded generalized
*  eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,
*  such that C has the same bandwidth as A.
*
*  B must have been previously factorized as S**T*S by SPBSTF, using a
*  split Cholesky factorization. A is overwritten by C = X**T*A*X, where
*  X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
*  bandwidth of A.
*
*  Arguments
*  =========
*
*  VECT    (input) CHARACTER*1
*          = 'N':  do not form the transformation matrix X;
*          = 'V':  form X.
*
*  UPLO    (input) CHARACTER*1
*          = 'U':  Upper triangle of A is stored;
*          = 'L':  Lower triangle of A is stored.
*
*  N       (input) INTEGER
*          The order of the matrices A and B.  N >= 0.
*
*  KA      (input) INTEGER
*          The number of superdiagonals of the matrix A if UPLO = 'U',
*          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
*
*  KB      (input) INTEGER
*          The number of superdiagonals of the matrix B if UPLO = 'U',
*          or the number of subdiagonals if UPLO = 'L'.  KA >= KB >= 0.
*
*  AB      (input/output) REAL array, dimension (LDAB,N)
*          On entry, the upper or lower triangle of the symmetric band
*          matrix A, stored in the first ka+1 rows of the array.  The
*          j-th column of A is stored in the j-th column of the array AB
*          as follows:
*          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
*          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
*
*          On exit, the transformed matrix X**T*A*X, stored in the same
*          format as A.
*
*  LDAB    (input) INTEGER
*          The leading dimension of the array AB.  LDAB >= KA+1.
*
*  BB      (input) REAL array, dimension (LDBB,N)
*          The banded factor S from the split Cholesky factorization of
*          B, as returned by SPBSTF, stored in the first KB+1 rows of
*          the array.
*
*  LDBB    (input) INTEGER
*          The leading dimension of the array BB.  LDBB >= KB+1.
*
*  X       (output) REAL array, dimension (LDX,N)
*          If VECT = 'V', the n-by-n matrix X.
*          If VECT = 'N', the array X is not referenced.
*
*  LDX     (input) INTEGER
*          The leading dimension of the array X.
*          LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
*
*  WORK    (workspace) REAL array, dimension (2*N)
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            UPDATE, UPPER, WANTX
      INTEGER            I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
     $                   KA1, KB1, KBT, L, M, NR, NRT, NX
      REAL               BII, RA, RA1, T
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGER, SLAR2V, SLARGV, SLARTG, SLARTV, SLASET,
     $                   SROT, SSCAL, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAXMIN
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters
*
      WANTX = LSAME( VECT, 'V' )
      UPPER = LSAME( UPLO, 'U' )
      KA1 = KA + 1
      KB1 = KB + 1
      INFO = 0
      IF.NOT.WANTX .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
         INFO = -1
      ELSE IF.NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( KA.LT.0 ) THEN
         INFO = -4
      ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
         INFO = -5
      ELSE IF( LDAB.LT.KA+1 ) THEN
         INFO = -7
      ELSE IF( LDBB.LT.KB+1 ) THEN
         INFO = -9
      ELSE IF( LDX.LT.1 .OR. WANTX .AND. LDX.LT.MAX1, N ) ) THEN
         INFO = -11
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SSBGST'-INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( N.EQ.0 )
     $   RETURN
*
      INCA = LDAB*KA1
*
*     Initialize X to the unit matrix, if needed
*
      IF( WANTX )
     $   CALL SLASET( 'Full', N, N, ZERO, ONE, X, LDX )
*
*     Set M to the splitting point m. It must be the same value as is
*     used in SPBSTF. The chosen value allows the arrays WORK and RWORK
*     to be of dimension (N).
*
      M = ( N+KB ) / 2
*
*     The routine works in two phases, corresponding to the two halves
*     of the split Cholesky factorization of B as S**T*S where
*
*     S = ( U    )
*         ( M  L )
*
*     with U upper triangular of order m, and L lower triangular of
*     order n-m. S has the same bandwidth as B.
*
*     S is treated as a product of elementary matrices:
*
*     S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n)
*
*     where S(i) is determined by the i-th row of S.
*
*     In phase 1, the index i takes the values n, n-1, ... , m+1;
*     in phase 2, it takes the values 1, 2, ... , m.
*
*     For each value of i, the current matrix A is updated by forming
*     inv(S(i))**T*A*inv(S(i)). This creates a triangular bulge outside
*     the band of A. The bulge is then pushed down toward the bottom of
*     A in phase 1, and up toward the top of A in phase 2, by applying
*     plane rotations.
*
*     There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
*     of them are linearly independent, so annihilating a bulge requires
*     only 2*kb-1 plane rotations. The rotations are divided into a 1st
*     set of kb-1 rotations, and a 2nd set of kb rotations.
*
*     Wherever possible, rotations are generated and applied in vector
*     operations of length NR between the indices J1 and J2 (sometimes
*     replaced by modified values NRT, J1T or J2T).
*
*     The cosines and sines of the rotations are stored in the array
*     WORK. The cosines of the 1st set of rotations are stored in
*     elements n+2:n+m-kb-1 and the sines of the 1st set in elements
*     2:m-kb-1; the cosines of the 2nd set are stored in elements
*     n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n.
*
*     The bulges are not formed explicitly; nonzero elements outside the
*     band are created only when they are required for generating new
*     rotations; they are stored in the array WORK, in positions where
*     they are later overwritten by the sines of the rotations which
*     annihilate them.
*
*     **************************** Phase 1 *****************************
*
*     The logical structure of this phase is:
*
*     UPDATE = .TRUE.
*     DO I = N, M + 1, -1
*        use S(i) to update A and create a new bulge
*        apply rotations to push all bulges KA positions downward
*     END DO
*     UPDATE = .FALSE.
*     DO I = M + KA + 1, N - 1
*        apply rotations to push all bulges KA positions downward
*     END DO
*
*     To avoid duplicating code, the two loops are merged.
*
      UPDATE = .TRUE.
      I = N + 1
   10 CONTINUE
      IF( UPDATE ) THEN
         I = I - 1
         KBT = MIN( KB, I-1 )
         I0 = I - 1
         I1 = MIN( N, I+KA )
         I2 = I - KBT + KA1
         IF( I.LT.M+1 ) THEN
            UPDATE = .FALSE.
            I = I + 1
            I0 = M
            IF( KA.EQ.0 )
     $         GO TO 480
            GO TO 10
         END IF
      ELSE
         I = I + KA
         IF( I.GT.N-1 )
     $      GO TO 480
      END IF
*
      IF( UPPER ) THEN
*
*        Transform A, working with the upper triangle
*
         IF( UPDATE ) THEN
*
*           Form  inv(S(i))**T * A * inv(S(i))
*
            BII = BB( KB1, I )
            DO 20 J = I, I1
               AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
   20       CONTINUE
            DO 30 J = MAX1, I-KA ), I
               AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
   30       CONTINUE
            DO 60 K = I - KBT, I - 1
               DO 40 J = I - KBT, K
                  AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
     $                               BB( J-I+KB1, I )*AB( K-I+KA1, I ) -
     $                               BB( K-I+KB1, I )*AB( J-I+KA1, I ) +
     $                               AB( KA1, I )*BB( J-I+KB1, I )*
     $                               BB( K-I+KB1, I )
   40          CONTINUE
               DO 50 J = MAX1, I-KA ), I - KBT - 1
                  AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
     $                               BB( K-I+KB1, I )*AB( J-I+KA1, I )
   50          CONTINUE
   60       CONTINUE
            DO 80 J = I, I1
               DO 70 K = MAX( J-KA, I-KBT ), I - 1
                  AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
     $                               BB( K-I+KB1, I )*AB( I-J+KA1, J )
   70          CONTINUE
   80       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by inv(S(i))
*
               CALL SSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
               IF( KBT.GT.0 )
     $            CALL SGER( N-M, KBT, -ONE, X( M+1, I ), 1,
     $                       BB( KB1-KBT, I ), 1, X( M+1, I-KBT ), LDX )
            END IF
*
*           store a(i,i1) in RA1 for use in next loop over K
*
            RA1 = AB( I-I1+KA1, I1 )
         END IF
*
*        Generate and apply vectors of rotations to chase all the
*        existing bulges KA positions down toward the bottom of the
*        band
*
         DO 130 K = 1, KB - 1
            IF( UPDATE ) THEN
*
*              Determine the rotations which would annihilate the bulge
*              which has in theory just been created
*
               IF( I-K+KA.LT..AND. I-K.GT.1 ) THEN
*
*                 generate rotation to annihilate a(i,i-k+ka+1)
*
                  CALL SLARTG( AB( K+1, I-K+KA ), RA1,
     $                         WORK( N+I-K+KA-M ), WORK( I-K+KA-M ),
     $                         RA )
*
*                 create nonzero element a(i-k,i-k+ka+1) outside the
*                 band and store it in WORK(i-k)
*
                  T = -BB( KB1-K, I )*RA1
                  WORK( I-K ) = WORK( N+I-K+KA-M )*-
     $                          WORK( I-K+KA-M )*AB( 1, I-K+KA )
                  AB( 1, I-K+KA ) = WORK( I-K+KA-M )*+
     $                              WORK( N+I-K+KA-M )*AB( 1, I-K+KA )
                  RA1 = RA
               END IF
            END IF
            J2 = I - K - 1 + MAX1, K-I0+2 )*KA1
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            IF( UPDATE ) THEN
               J2T = MAX( J2, I+2*KA-K+1 )
            ELSE
               J2T = J2
            END IF
            NRT = ( N-J2T+KA ) / KA1
            DO 90 J = J2T, J1, KA1
*
*              create nonzero element a(j-ka,j+1) outside the band
*              and store it in WORK(j-m)
*
               WORK( J-M ) = WORK( J-M )*AB( 1, J+1 )
               AB( 1, J+1 ) = WORK( N+J-M )*AB( 1, J+1 )
   90       CONTINUE
*
*           generate rotations in 1st set to annihilate elements which
*           have been created outside the band
*
            IF( NRT.GT.0 )
     $         CALL SLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1,
     $                      WORK( N+J2T-M ), KA1 )
            IF( NR.GT.0 ) THEN
*
*              apply rotations in 1st set from the right
*
               DO 100 L = 1, KA - 1
                  CALL SLARTV( NR, AB( KA1-L, J2 ), INCA,
     $                         AB( KA-L, J2+1 ), INCA, WORK( N+J2-M ),
     $                         WORK( J2-M ), KA1 )
  100          CONTINUE
*
*              apply rotations in 1st set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
     $                      AB( KA, J2+1 ), INCA, WORK( N+J2-M ),
     $                      WORK( J2-M ), KA1 )
*
            END IF
*
*           start applying rotations in 1st set from the left
*
            DO 110 L = KA - 1, KB - K + 1-1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J2+KA1-L ), INCA,
     $                         AB( L+1, J2+KA1-L ), INCA,
     $                         WORK( N+J2-M ), WORK( J2-M ), KA1 )
  110       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 1st set
*
               DO 120 J = J2, J1, KA1
                  CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       WORK( N+J-M ), WORK( J-M ) )
  120          CONTINUE
            END IF
  130    CONTINUE
*
         IF( UPDATE ) THEN
            IF( I2.LE..AND. KBT.GT.0 ) THEN
*
*              create nonzero element a(i-kbt,i-kbt+ka+1) outside the
*              band and store it in WORK(i-kbt)
*
               WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1
            END IF
         END IF
*
         DO 170 K = KB, 1-1
            IF( UPDATE ) THEN
               J2 = I - K - 1 + MAX2, K-I0+1 )*KA1
            ELSE
               J2 = I - K - 1 + MAX1, K-I0+1 )*KA1
            END IF
*
*           finish applying rotations in 2nd set from the left
*
            DO 140 L = KB - K, 1-1
               NRT = ( N-J2+KA+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J2-L+1 ), INCA,
     $                         AB( L+1, J2-L+1 ), INCA, WORK( N+J2-KA ),
     $                         WORK( J2-KA ), KA1 )
  140       CONTINUE
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            DO 150 J = J1, J2, -KA1
               WORK( J ) = WORK( J-KA )
               WORK( N+J ) = WORK( N+J-KA )
  150       CONTINUE
            DO 160 J = J2, J1, KA1
*
*              create nonzero element a(j-ka,j+1) outside the band
*              and store it in WORK(j)
*
               WORK( J ) = WORK( J )*AB( 1, J+1 )
               AB( 1, J+1 ) = WORK( N+J )*AB( 1, J+1 )
  160       CONTINUE
            IF( UPDATE ) THEN
               IF( I-K.LT.N-KA .AND. K.LE.KBT )
     $            WORK( I-K+KA ) = WORK( I-K )
            END IF
  170    CONTINUE
*
         DO 210 K = KB, 1-1
            J2 = I - K - 1 + MAX1, K-I0+1 )*KA1
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            IF( NR.GT.0 ) THEN
*
*              generate rotations in 2nd set to annihilate elements
*              which have been created outside the band
*
               CALL SLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1,
     $                      WORK( N+J2 ), KA1 )
*
*              apply rotations in 2nd set from the right
*
               DO 180 L = 1, KA - 1
                  CALL SLARTV( NR, AB( KA1-L, J2 ), INCA,
     $                         AB( KA-L, J2+1 ), INCA, WORK( N+J2 ),
     $                         WORK( J2 ), KA1 )
  180          CONTINUE
*
*              apply rotations in 2nd set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
     $                      AB( KA, J2+1 ), INCA, WORK( N+J2 ),
     $                      WORK( J2 ), KA1 )
*
            END IF
*
*           start applying rotations in 2nd set from the left
*
            DO 190 L = KA - 1, KB - K + 1-1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J2+KA1-L ), INCA,
     $                         AB( L+1, J2+KA1-L ), INCA, WORK( N+J2 ),
     $                         WORK( J2 ), KA1 )
  190       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 2nd set
*
               DO 200 J = J2, J1, KA1
                  CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       WORK( N+J ), WORK( J ) )
  200          CONTINUE
            END IF
  210    CONTINUE
*
         DO 230 K = 1, KB - 1
            J2 = I - K - 1 + MAX1, K-I0+2 )*KA1
*
*           finish applying rotations in 1st set from the left
*
            DO 220 L = KB - K, 1-1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J2+KA1-L ), INCA,
     $                         AB( L+1, J2+KA1-L ), INCA,
     $                         WORK( N+J2-M ), WORK( J2-M ), KA1 )
  220       CONTINUE
  230    CONTINUE
*
         IF( KB.GT.1 ) THEN
            DO 240 J = N - 1, I - KB + 2*KA + 1-1
               WORK( N+J-M ) = WORK( N+J-KA-M )
               WORK( J-M ) = WORK( J-KA-M )
  240       CONTINUE
         END IF
*
      ELSE
*
*        Transform A, working with the lower triangle
*
         IF( UPDATE ) THEN
*
*           Form  inv(S(i))**T * A * inv(S(i))
*
            BII = BB( 1, I )
            DO 250 J = I, I1
               AB( J-I+1, I ) = AB( J-I+1, I ) / BII
  250       CONTINUE
            DO 260 J = MAX1, I-KA ), I
               AB( I-J+1, J ) = AB( I-J+1, J ) / BII
  260       CONTINUE
            DO 290 K = I - KBT, I - 1
               DO 270 J = I - KBT, K
                  AB( K-J+1, J ) = AB( K-J+1, J ) -
     $                             BB( I-J+1, J )*AB( I-K+1, K ) -
     $                             BB( I-K+1, K )*AB( I-J+1, J ) +
     $                             AB( 1, I )*BB( I-J+1, J )*
     $                             BB( I-K+1, K )
  270          CONTINUE
               DO 280 J = MAX1, I-KA ), I - KBT - 1
                  AB( K-J+1, J ) = AB( K-J+1, J ) -
     $                             BB( I-K+1, K )*AB( I-J+1, J )
  280          CONTINUE
  290       CONTINUE
            DO 310 J = I, I1
               DO 300 K = MAX( J-KA, I-KBT ), I - 1
                  AB( J-K+1, K ) = AB( J-K+1, K ) -
     $                             BB( I-K+1, K )*AB( J-I+1, I )
  300          CONTINUE
  310       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by inv(S(i))
*
               CALL SSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
               IF( KBT.GT.0 )
     $            CALL SGER( N-M, KBT, -ONE, X( M+1, I ), 1,
     $                       BB( KBT+1, I-KBT ), LDBB-1,
     $                       X( M+1, I-KBT ), LDX )
            END IF
*
*           store a(i1,i) in RA1 for use in next loop over K
*
            RA1 = AB( I1-I+1, I )
         END IF
*
*        Generate and apply vectors of rotations to chase all the
*        existing bulges KA positions down toward the bottom of the
*        band
*
         DO 360 K = 1, KB - 1
            IF( UPDATE ) THEN
*
*              Determine the rotations which would annihilate the bulge
*              which has in theory just been created
*
               IF( I-K+KA.LT..AND. I-K.GT.1 ) THEN
*
*                 generate rotation to annihilate a(i-k+ka+1,i)
*
                  CALL SLARTG( AB( KA1-K, I ), RA1, WORK( N+I-K+KA-M ),
     $                         WORK( I-K+KA-M ), RA )
*
*                 create nonzero element a(i-k+ka+1,i-k) outside the
*                 band and store it in WORK(i-k)
*
                  T = -BB( K+1, I-K )*RA1
                  WORK( I-K ) = WORK( N+I-K+KA-M )*-
     $                          WORK( I-K+KA-M )*AB( KA1, I-K )
                  AB( KA1, I-K ) = WORK( I-K+KA-M )*+
     $                             WORK( N+I-K+KA-M )*AB( KA1, I-K )
                  RA1 = RA
               END IF
            END IF
            J2 = I - K - 1 + MAX1, K-I0+2 )*KA1
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            IF( UPDATE ) THEN
               J2T = MAX( J2, I+2*KA-K+1 )
            ELSE
               J2T = J2
            END IF
            NRT = ( N-J2T+KA ) / KA1
            DO 320 J = J2T, J1, KA1
*
*              create nonzero element a(j+1,j-ka) outside the band
*              and store it in WORK(j-m)
*
               WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
               AB( KA1, J-KA+1 ) = WORK( N+J-M )*AB( KA1, J-KA+1 )
  320       CONTINUE
*
*           generate rotations in 1st set to annihilate elements which
*           have been created outside the band
*
            IF( NRT.GT.0 )
     $         CALL SLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
     $                      KA1, WORK( N+J2T-M ), KA1 )
            IF( NR.GT.0 ) THEN
*
*              apply rotations in 1st set from the left
*
               DO 330 L = 1, KA - 1
                  CALL SLARTV( NR, AB( L+1, J2-L ), INCA,
     $                         AB( L+2, J2-L ), INCA, WORK( N+J2-M ),
     $                         WORK( J2-M ), KA1 )
  330          CONTINUE
*
*              apply rotations in 1st set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
     $                      INCA, WORK( N+J2-M ), WORK( J2-M ), KA1 )
*
            END IF
*
*           start applying rotations in 1st set from the right
*
            DO 340 L = KA - 1, KB - K + 1-1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
     $                         AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
     $                         WORK( J2-M ), KA1 )
  340       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 1st set
*
               DO 350 J = J2, J1, KA1
                  CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       WORK( N+J-M ), WORK( J-M ) )
  350          CONTINUE
            END IF
  360    CONTINUE
*
         IF( UPDATE ) THEN
            IF( I2.LE..AND. KBT.GT.0 ) THEN
*
*              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
*              band and store it in WORK(i-kbt)
*
               WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
            END IF
         END IF
*
         DO 400 K = KB, 1-1
            IF( UPDATE ) THEN
               J2 = I - K - 1 + MAX2, K-I0+1 )*KA1
            ELSE
               J2 = I - K - 1 + MAX1, K-I0+1 )*KA1
            END IF
*
*           finish applying rotations in 2nd set from the right
*
            DO 370 L = KB - K, 1-1
               NRT = ( N-J2+KA+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
     $                         AB( KA1-L, J2-KA+1 ), INCA,
     $                         WORK( N+J2-KA ), WORK( J2-KA ), KA1 )
  370       CONTINUE
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            DO 380 J = J1, J2, -KA1
               WORK( J ) = WORK( J-KA )
               WORK( N+J ) = WORK( N+J-KA )
  380       CONTINUE
            DO 390 J = J2, J1, KA1
*
*              create nonzero element a(j+1,j-ka) outside the band
*              and store it in WORK(j)
*
               WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
               AB( KA1, J-KA+1 ) = WORK( N+J )*AB( KA1, J-KA+1 )
  390       CONTINUE
            IF( UPDATE ) THEN
               IF( I-K.LT.N-KA .AND. K.LE.KBT )
     $            WORK( I-K+KA ) = WORK( I-K )
            END IF
  400    CONTINUE
*
         DO 440 K = KB, 1-1
            J2 = I - K - 1 + MAX1, K-I0+1 )*KA1
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            IF( NR.GT.0 ) THEN
*
*              generate rotations in 2nd set to annihilate elements
*              which have been created outside the band
*
               CALL SLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
     $                      WORK( N+J2 ), KA1 )
*
*              apply rotations in 2nd set from the left
*
               DO 410 L = 1, KA - 1
                  CALL SLARTV( NR, AB( L+1, J2-L ), INCA,
     $                         AB( L+2, J2-L ), INCA, WORK( N+J2 ),
     $                         WORK( J2 ), KA1 )
  410          CONTINUE
*
*              apply rotations in 2nd set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
     $                      INCA, WORK( N+J2 ), WORK( J2 ), KA1 )
*
            END IF
*
*           start applying rotations in 2nd set from the right
*
            DO 420 L = KA - 1, KB - K + 1-1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
     $                         AB( KA1-L, J2+1 ), INCA, WORK( N+J2 ),
     $                         WORK( J2 ), KA1 )
  420       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 2nd set
*
               DO 430 J = J2, J1, KA1
                  CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       WORK( N+J ), WORK( J ) )
  430          CONTINUE
            END IF
  440    CONTINUE
*
         DO 460 K = 1, KB - 1
            J2 = I - K - 1 + MAX1, K-I0+2 )*KA1
*
*           finish applying rotations in 1st set from the right
*
            DO 450 L = KB - K, 1-1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
     $                         AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
     $                         WORK( J2-M ), KA1 )
  450       CONTINUE
  460    CONTINUE
*
         IF( KB.GT.1 ) THEN
            DO 470 J = N - 1, I - KB + 2*KA + 1-1
               WORK( N+J-M ) = WORK( N+J-KA-M )
               WORK( J-M ) = WORK( J-KA-M )
  470       CONTINUE
         END IF
*
      END IF
*
      GO TO 10
*
  480 CONTINUE
*
*     **************************** Phase 2 *****************************
*
*     The logical structure of this phase is:
*
*     UPDATE = .TRUE.
*     DO I = 1, M
*        use S(i) to update A and create a new bulge
*        apply rotations to push all bulges KA positions upward
*     END DO
*     UPDATE = .FALSE.
*     DO I = M - KA - 1, 2, -1
*        apply rotations to push all bulges KA positions upward
*     END DO
*
*     To avoid duplicating code, the two loops are merged.
*
      UPDATE = .TRUE.
      I = 0
  490 CONTINUE
      IF( UPDATE ) THEN
         I = I + 1
         KBT = MIN( KB, M-I )
         I0 = I + 1
         I1 = MAX1, I-KA )
         I2 = I + KBT - KA1
         IF( I.GT.M ) THEN
            UPDATE = .FALSE.
            I = I - 1
            I0 = M + 1
            IF( KA.EQ.0 )
     $         RETURN
            GO TO 490
         END IF
      ELSE
         I = I - KA
         IF( I.LT.2 )
     $      RETURN
      END IF
*
      IF( I.LT.M-KBT ) THEN
         NX = M
      ELSE
         NX = N
      END IF
*
      IF( UPPER ) THEN
*
*        Transform A, working with the upper triangle
*
         IF( UPDATE ) THEN
*
*           Form  inv(S(i))**T * A * inv(S(i))
*
            BII = BB( KB1, I )
            DO 500 J = I1, I
               AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
  500       CONTINUE
            DO 510 J = I, MIN( N, I+KA )
               AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
  510       CONTINUE
            DO 540 K = I + 1, I + KBT
               DO 520 J = K, I + KBT
                  AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
     $                               BB( I-J+KB1, J )*AB( I-K+KA1, K ) -
     $                               BB( I-K+KB1, K )*AB( I-J+KA1, J ) +
     $                               AB( KA1, I )*BB( I-J+KB1, J )*
     $                               BB( I-K+KB1, K )
  520          CONTINUE
               DO 530 J = I + KBT + 1MIN( N, I+KA )
                  AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
     $                               BB( I-K+KB1, K )*AB( I-J+KA1, J )
  530          CONTINUE
  540       CONTINUE
            DO 560 J = I1, I
               DO 550 K = I + 1MIN( J+KA, I+KBT )
                  AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
     $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
  550          CONTINUE
  560       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by inv(S(i))
*
               CALL SSCAL( NX, ONE / BII, X( 1, I ), 1 )
               IF( KBT.GT.0 )
     $            CALL SGER( NX, KBT, -ONE, X( 1, I ), 1, BB( KB, I+1 ),
     $                       LDBB-1, X( 1, I+1 ), LDX )
            END IF
*
*           store a(i1,i) in RA1 for use in next loop over K
*
            RA1 = AB( I1-I+KA1, I )
         END IF
*
*        Generate and apply vectors of rotations to chase all the
*        existing bulges KA positions up toward the top of the band
*
         DO 610 K = 1, KB - 1
            IF( UPDATE ) THEN
*
*              Determine the rotations which would annihilate the bulge
*              which has in theory just been created
*
               IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
*
*                 generate rotation to annihilate a(i+k-ka-1,i)
*
                  CALL SLARTG( AB( K+1, I ), RA1, WORK( N+I+K-KA ),
     $                         WORK( I+K-KA ), RA )
*
*                 create nonzero element a(i+k-ka-1,i+k) outside the
*                 band and store it in WORK(m-kb+i+k)
*
                  T = -BB( KB1-K, I+K )*RA1
                  WORK( M-KB+I+K ) = WORK( N+I+K-KA )*-
     $                               WORK( I+K-KA )*AB( 1, I+K )
                  AB( 1, I+K ) = WORK( I+K-KA )*+
     $                           WORK( N+I+K-KA )*AB( 1, I+K )
                  RA1 = RA
               END IF
            END IF
            J2 = I + K + 1 - MAX1, K+I0-M+1 )*KA1
            NR = ( J2+KA-1 ) / KA1
            J1 = J2 - ( NR-1 )*KA1
            IF( UPDATE ) THEN
               J2T = MIN( J2, I-2*KA+K-1 )
            ELSE
               J2T = J2
            END IF
            NRT = ( J2T+KA-1 ) / KA1
            DO 570 J = J1, J2T, KA1
*
*              create nonzero element a(j-1,j+ka) outside the band
*              and store it in WORK(j)
*
               WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
               AB( 1, J+KA-1 ) = WORK( N+J )*AB( 1, J+KA-1 )
  570       CONTINUE
*
*           generate rotations in 1st set to annihilate elements which
*           have been created outside the band
*
            IF( NRT.GT.0 )
     $         CALL SLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
     $                      WORK( N+J1 ), KA1 )
            IF( NR.GT.0 ) THEN
*
*              apply rotations in 1st set from the left
*
               DO 580 L = 1, KA - 1
                  CALL SLARTV( NR, AB( KA1-L, J1+L ), INCA,
     $                         AB( KA-L, J1+L ), INCA, WORK( N+J1 ),
     $                         WORK( J1 ), KA1 )
  580          CONTINUE
*
*              apply rotations in 1st set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
     $                      AB( KA, J1 ), INCA, WORK( N+J1 ),
     $                      WORK( J1 ), KA1 )
*
            END IF
*
*           start applying rotations in 1st set from the right
*
            DO 590 L = KA - 1, KB - K + 1-1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J1T ), INCA,
     $                         AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
     $                         WORK( J1T ), KA1 )
  590       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 1st set
*
               DO 600 J = J1, J2, KA1
                  CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
     $                       WORK( N+J ), WORK( J ) )
  600          CONTINUE
            END IF
  610    CONTINUE
*
         IF( UPDATE ) THEN
            IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
*
*              create nonzero element a(i+kbt-ka-1,i+kbt) outside the
*              band and store it in WORK(m-kb+i+kbt)
*
               WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
            END IF
         END IF
*
         DO 650 K = KB, 1-1
            IF( UPDATE ) THEN
               J2 = I + K + 1 - MAX2, K+I0-M )*KA1
            ELSE
               J2 = I + K + 1 - MAX1, K+I0-M )*KA1
            END IF
*
*           finish applying rotations in 2nd set from the right
*
            DO 620 L = KB - K, 1-1
               NRT = ( J2+KA+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J1T+KA ), INCA,
     $                         AB( L+1, J1T+KA-1 ), INCA,
     $                         WORK( N+M-KB+J1T+KA ),
     $                         WORK( M-KB+J1T+KA ), KA1 )
  620       CONTINUE
            NR = ( J2+KA-1 ) / KA1
            J1 = J2 - ( NR-1 )*KA1
            DO 630 J = J1, J2, KA1
               WORK( M-KB+J ) = WORK( M-KB+J+KA )
               WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
  630       CONTINUE
            DO 640 J = J1, J2, KA1
*
*              create nonzero element a(j-1,j+ka) outside the band
*              and store it in WORK(m-kb+j)
*
               WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
               AB( 1, J+KA-1 ) = WORK( N+M-KB+J )*AB( 1, J+KA-1 )
  640       CONTINUE
            IF( UPDATE ) THEN
               IF( I+K.GT.KA1 .AND. K.LE.KBT )
     $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
            END IF
  650    CONTINUE
*
         DO 690 K = KB, 1-1
            J2 = I + K + 1 - MAX1, K+I0-M )*KA1
            NR = ( J2+KA-1 ) / KA1
            J1 = J2 - ( NR-1 )*KA1
            IF( NR.GT.0 ) THEN
*
*              generate rotations in 2nd set to annihilate elements
*              which have been created outside the band
*
               CALL SLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
     $                      KA1, WORK( N+M-KB+J1 ), KA1 )
*
*              apply rotations in 2nd set from the left
*
               DO 660 L = 1, KA - 1
                  CALL SLARTV( NR, AB( KA1-L, J1+L ), INCA,
     $                         AB( KA-L, J1+L ), INCA,
     $                         WORK( N+M-KB+J1 ), WORK( M-KB+J1 ), KA1 )
  660          CONTINUE
*
*              apply rotations in 2nd set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
     $                      AB( KA, J1 ), INCA, WORK( N+M-KB+J1 ),
     $                      WORK( M-KB+J1 ), KA1 )
*
            END IF
*
*           start applying rotations in 2nd set from the right
*
            DO 670 L = KA - 1, KB - K + 1-1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J1T ), INCA,
     $                         AB( L+1, J1T-1 ), INCA,
     $                         WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
     $                         KA1 )
  670       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 2nd set
*
               DO 680 J = J1, J2, KA1
                  CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
     $                       WORK( N+M-KB+J ), WORK( M-KB+J ) )
  680          CONTINUE
            END IF
  690    CONTINUE
*
         DO 710 K = 1, KB - 1
            J2 = I + K + 1 - MAX1, K+I0-M+1 )*KA1
*
*           finish applying rotations in 1st set from the right
*
            DO 700 L = KB - K, 1-1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J1T ), INCA,
     $                         AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
     $                         WORK( J1T ), KA1 )
  700       CONTINUE
  710    CONTINUE
*
         IF( KB.GT.1 ) THEN
            DO 720 J = 2MIN( I+KB, M ) - 2*KA - 1
               WORK( N+J ) = WORK( N+J+KA )
               WORK( J ) = WORK( J+KA )
  720       CONTINUE
         END IF
*
      ELSE
*
*        Transform A, working with the lower triangle
*
         IF( UPDATE ) THEN
*
*           Form  inv(S(i))**T * A * inv(S(i))
*
            BII = BB( 1, I )
            DO 730 J = I1, I
               AB( I-J+1, J ) = AB( I-J+1, J ) / BII
  730       CONTINUE
            DO 740 J = I, MIN( N, I+KA )
               AB( J-I+1, I ) = AB( J-I+1, I ) / BII
  740       CONTINUE
            DO 770 K = I + 1, I + KBT
               DO 750 J = K, I + KBT
                  AB( J-K+1, K ) = AB( J-K+1, K ) -
     $                             BB( J-I+1, I )*AB( K-I+1, I ) -
     $                             BB( K-I+1, I )*AB( J-I+1, I ) +
     $                             AB( 1, I )*BB( J-I+1, I )*
     $                             BB( K-I+1, I )
  750          CONTINUE
               DO 760 J = I + KBT + 1MIN( N, I+KA )
                  AB( J-K+1, K ) = AB( J-K+1, K ) -
     $                             BB( K-I+1, I )*AB( J-I+1, I )
  760          CONTINUE
  770       CONTINUE
            DO 790 J = I1, I
               DO 780 K = I + 1MIN( J+KA, I+KBT )
                  AB( K-J+1, J ) = AB( K-J+1, J ) -
     $                             BB( K-I+1, I )*AB( I-J+1, J )
  780          CONTINUE
  790       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by inv(S(i))
*
               CALL SSCAL( NX, ONE / BII, X( 1, I ), 1 )
               IF( KBT.GT.0 )
     $            CALL SGER( NX, KBT, -ONE, X( 1, I ), 1, BB( 2, I ), 1,
     $                       X( 1, I+1 ), LDX )
            END IF
*
*           store a(i,i1) in RA1 for use in next loop over K
*
            RA1 = AB( I-I1+1, I1 )
         END IF
*
*        Generate and apply vectors of rotations to chase all the
*        existing bulges KA positions up toward the top of the band
*
         DO 840 K = 1, KB - 1
            IF( UPDATE ) THEN
*
*              Determine the rotations which would annihilate the bulge
*              which has in theory just been created
*
               IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
*
*                 generate rotation to annihilate a(i,i+k-ka-1)
*
                  CALL SLARTG( AB( KA1-K, I+K-KA ), RA1,
     $                         WORK( N+I+K-KA ), WORK( I+K-KA ), RA )
*
*                 create nonzero element a(i+k,i+k-ka-1) outside the
*                 band and store it in WORK(m-kb+i+k)
*
                  T = -BB( K+1, I )*RA1
                  WORK( M-KB+I+K ) = WORK( N+I+K-KA )*-
     $                               WORK( I+K-KA )*AB( KA1, I+K-KA )
                  AB( KA1, I+K-KA ) = WORK( I+K-KA )*+
     $                                WORK( N+I+K-KA )*AB( KA1, I+K-KA )
                  RA1 = RA
               END IF
            END IF
            J2 = I + K + 1 - MAX1, K+I0-M+1 )*KA1
            NR = ( J2+KA-1 ) / KA1
            J1 = J2 - ( NR-1 )*KA1
            IF( UPDATE ) THEN
               J2T = MIN( J2, I-2*KA+K-1 )
            ELSE
               J2T = J2
            END IF
            NRT = ( J2T+KA-1 ) / KA1
            DO 800 J = J1, J2T, KA1
*
*              create nonzero element a(j+ka,j-1) outside the band
*              and store it in WORK(j)
*
               WORK( J ) = WORK( J )*AB( KA1, J-1 )
               AB( KA1, J-1 ) = WORK( N+J )*AB( KA1, J-1 )
  800       CONTINUE
*
*           generate rotations in 1st set to annihilate elements which
*           have been created outside the band
*
            IF( NRT.GT.0 )
     $         CALL SLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
     $                      WORK( N+J1 ), KA1 )
            IF( NR.GT.0 ) THEN
*
*              apply rotations in 1st set from the right
*
               DO 810 L = 1, KA - 1
                  CALL SLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
     $                         INCA, WORK( N+J1 ), WORK( J1 ), KA1 )
  810          CONTINUE
*
*              apply rotations in 1st set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
     $                      AB( 2, J1-1 ), INCA, WORK( N+J1 ),
     $                      WORK( J1 ), KA1 )
*
            END IF
*
*           start applying rotations in 1st set from the left
*
            DO 820 L = KA - 1, KB - K + 1-1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
     $                         AB( KA1-L, J1T-KA1+L ), INCA,
     $                         WORK( N+J1T ), WORK( J1T ), KA1 )
  820       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 1st set
*
               DO 830 J = J1, J2, KA1
                  CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
     $                       WORK( N+J ), WORK( J ) )
  830          CONTINUE
            END IF
  840    CONTINUE
*
         IF( UPDATE ) THEN
            IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
*
*              create nonzero element a(i+kbt,i+kbt-ka-1) outside the
*              band and store it in WORK(m-kb+i+kbt)
*
               WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
            END IF
         END IF
*
         DO 880 K = KB, 1-1
            IF( UPDATE ) THEN
               J2 = I + K + 1 - MAX2, K+I0-M )*KA1
            ELSE
               J2 = I + K + 1 - MAX1, K+I0-M )*KA1
            END IF
*
*           finish applying rotations in 2nd set from the left
*
            DO 850 L = KB - K, 1-1
               NRT = ( J2+KA+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
     $                         AB( KA1-L, J1T+L-1 ), INCA,
     $                         WORK( N+M-KB+J1T+KA ),
     $                         WORK( M-KB+J1T+KA ), KA1 )
  850       CONTINUE
            NR = ( J2+KA-1 ) / KA1
            J1 = J2 - ( NR-1 )*KA1
            DO 860 J = J1, J2, KA1
               WORK( M-KB+J ) = WORK( M-KB+J+KA )
               WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
  860       CONTINUE
            DO 870 J = J1, J2, KA1
*
*              create nonzero element a(j+ka,j-1) outside the band
*              and store it in WORK(m-kb+j)
*
               WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
               AB( KA1, J-1 ) = WORK( N+M-KB+J )*AB( KA1, J-1 )
  870       CONTINUE
            IF( UPDATE ) THEN
               IF( I+K.GT.KA1 .AND. K.LE.KBT )
     $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
            END IF
  880    CONTINUE
*
         DO 920 K = KB, 1-1
            J2 = I + K + 1 - MAX1, K+I0-M )*KA1
            NR = ( J2+KA-1 ) / KA1
            J1 = J2 - ( NR-1 )*KA1
            IF( NR.GT.0 ) THEN
*
*              generate rotations in 2nd set to annihilate elements
*              which have been created outside the band
*
               CALL SLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
     $                      KA1, WORK( N+M-KB+J1 ), KA1 )
*
*              apply rotations in 2nd set from the right
*
               DO 890 L = 1, KA - 1
                  CALL SLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
     $                         INCA, WORK( N+M-KB+J1 ), WORK( M-KB+J1 ),
     $                         KA1 )
  890          CONTINUE
*
*              apply rotations in 2nd set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
     $                      AB( 2, J1-1 ), INCA, WORK( N+M-KB+J1 ),
     $                      WORK( M-KB+J1 ), KA1 )
*
            END IF
*
*           start applying rotations in 2nd set from the left
*
            DO 900 L = KA - 1, KB - K + 1-1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
     $                         AB( KA1-L, J1T-KA1+L ), INCA,
     $                         WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
     $                         KA1 )
  900       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 2nd set
*
               DO 910 J = J1, J2, KA1
                  CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
     $                       WORK( N+M-KB+J ), WORK( M-KB+J ) )
  910          CONTINUE
            END IF
  920    CONTINUE
*
         DO 940 K = 1, KB - 1
            J2 = I + K + 1 - MAX1, K+I0-M+1 )*KA1
*
*           finish applying rotations in 1st set from the left
*
            DO 930 L = KB - K, 1-1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
     $                         AB( KA1-L, J1T-KA1+L ), INCA,
     $                         WORK( N+J1T ), WORK( J1T ), KA1 )
  930       CONTINUE
  940    CONTINUE
*
         IF( KB.GT.1 ) THEN
            DO 950 J = 2MIN( I+KB, M ) - 2*KA - 1
               WORK( N+J ) = WORK( N+J+KA )
               WORK( J ) = WORK( J+KA )
  950       CONTINUE
         END IF
*
      END IF
*
      GO TO 490
*
*     End of SSBGST
*
      END