LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dsbgst.f
Go to the documentation of this file.
1 *> \brief \b DSBGST
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DSBGST + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbgst.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbgst.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbgst.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DSBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
22 * LDX, WORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO, VECT
26 * INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N
27 * ..
28 * .. Array Arguments ..
29 * DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
30 * $ X( LDX, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DSBGST reduces a real symmetric-definite banded generalized
40 *> eigenproblem A*x = lambda*B*x to standard form C*y = lambda*y,
41 *> such that C has the same bandwidth as A.
42 *>
43 *> B must have been previously factorized as S**T*S by DPBSTF, using a
44 *> split Cholesky factorization. A is overwritten by C = X**T*A*X, where
45 *> X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
46 *> bandwidth of A.
47 *> \endverbatim
48 *
49 * Arguments:
50 * ==========
51 *
52 *> \param[in] VECT
53 *> \verbatim
54 *> VECT is CHARACTER*1
55 *> = 'N': do not form the transformation matrix X;
56 *> = 'V': form X.
57 *> \endverbatim
58 *>
59 *> \param[in] UPLO
60 *> \verbatim
61 *> UPLO is CHARACTER*1
62 *> = 'U': Upper triangle of A is stored;
63 *> = 'L': Lower triangle of A is stored.
64 *> \endverbatim
65 *>
66 *> \param[in] N
67 *> \verbatim
68 *> N is INTEGER
69 *> The order of the matrices A and B. N >= 0.
70 *> \endverbatim
71 *>
72 *> \param[in] KA
73 *> \verbatim
74 *> KA is INTEGER
75 *> The number of superdiagonals of the matrix A if UPLO = 'U',
76 *> or the number of subdiagonals if UPLO = 'L'. KA >= 0.
77 *> \endverbatim
78 *>
79 *> \param[in] KB
80 *> \verbatim
81 *> KB is INTEGER
82 *> The number of superdiagonals of the matrix B if UPLO = 'U',
83 *> or the number of subdiagonals if UPLO = 'L'. KA >= KB >= 0.
84 *> \endverbatim
85 *>
86 *> \param[in,out] AB
87 *> \verbatim
88 *> AB is DOUBLE PRECISION array, dimension (LDAB,N)
89 *> On entry, the upper or lower triangle of the symmetric band
90 *> matrix A, stored in the first ka+1 rows of the array. The
91 *> j-th column of A is stored in the j-th column of the array AB
92 *> as follows:
93 *> if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
94 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka).
95 *>
96 *> On exit, the transformed matrix X**T*A*X, stored in the same
97 *> format as A.
98 *> \endverbatim
99 *>
100 *> \param[in] LDAB
101 *> \verbatim
102 *> LDAB is INTEGER
103 *> The leading dimension of the array AB. LDAB >= KA+1.
104 *> \endverbatim
105 *>
106 *> \param[in] BB
107 *> \verbatim
108 *> BB is DOUBLE PRECISION array, dimension (LDBB,N)
109 *> The banded factor S from the split Cholesky factorization of
110 *> B, as returned by DPBSTF, stored in the first KB+1 rows of
111 *> the array.
112 *> \endverbatim
113 *>
114 *> \param[in] LDBB
115 *> \verbatim
116 *> LDBB is INTEGER
117 *> The leading dimension of the array BB. LDBB >= KB+1.
118 *> \endverbatim
119 *>
120 *> \param[out] X
121 *> \verbatim
122 *> X is DOUBLE PRECISION array, dimension (LDX,N)
123 *> If VECT = 'V', the n-by-n matrix X.
124 *> If VECT = 'N', the array X is not referenced.
125 *> \endverbatim
126 *>
127 *> \param[in] LDX
128 *> \verbatim
129 *> LDX is INTEGER
130 *> The leading dimension of the array X.
131 *> LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
132 *> \endverbatim
133 *>
134 *> \param[out] WORK
135 *> \verbatim
136 *> WORK is DOUBLE PRECISION array, dimension (2*N)
137 *> \endverbatim
138 *>
139 *> \param[out] INFO
140 *> \verbatim
141 *> INFO is INTEGER
142 *> = 0: successful exit
143 *> < 0: if INFO = -i, the i-th argument had an illegal value.
144 *> \endverbatim
145 *
146 * Authors:
147 * ========
148 *
149 *> \author Univ. of Tennessee
150 *> \author Univ. of California Berkeley
151 *> \author Univ. of Colorado Denver
152 *> \author NAG Ltd.
153 *
154 *> \date November 2011
155 *
156 *> \ingroup doubleOTHERcomputational
157 *
158 * =====================================================================
159  SUBROUTINE dsbgst( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
160  $ ldx, work, info )
161 *
162 * -- LAPACK computational routine (version 3.4.0) --
163 * -- LAPACK is a software package provided by Univ. of Tennessee, --
164 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165 * November 2011
166 *
167 * .. Scalar Arguments ..
168  CHARACTER uplo, vect
169  INTEGER info, ka, kb, ldab, ldbb, ldx, n
170 * ..
171 * .. Array Arguments ..
172  DOUBLE PRECISION ab( ldab, * ), bb( ldbb, * ), work( * ),
173  $ x( ldx, * )
174 * ..
175 *
176 * =====================================================================
177 *
178 * .. Parameters ..
179  DOUBLE PRECISION zero, one
180  parameter( zero = 0.0d+0, one = 1.0d+0 )
181 * ..
182 * .. Local Scalars ..
183  LOGICAL update, upper, wantx
184  INTEGER i, i0, i1, i2, inca, j, j1, j1t, j2, j2t, k,
185  $ ka1, kb1, kbt, l, m, nr, nrt, nx
186  DOUBLE PRECISION bii, ra, ra1, t
187 * ..
188 * .. External Functions ..
189  LOGICAL lsame
190  EXTERNAL lsame
191 * ..
192 * .. External Subroutines ..
193  EXTERNAL dger, dlar2v, dlargv, dlartg, dlartv, dlaset,
194  $ drot, dscal, xerbla
195 * ..
196 * .. Intrinsic Functions ..
197  INTRINSIC max, min
198 * ..
199 * .. Executable Statements ..
200 *
201 * Test the input parameters
202 *
203  wantx = lsame( vect, 'V' )
204  upper = lsame( uplo, 'U' )
205  ka1 = ka + 1
206  kb1 = kb + 1
207  info = 0
208  IF( .NOT.wantx .AND. .NOT.lsame( vect, 'N' ) ) THEN
209  info = -1
210  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
211  info = -2
212  ELSE IF( n.LT.0 ) THEN
213  info = -3
214  ELSE IF( ka.LT.0 ) THEN
215  info = -4
216  ELSE IF( kb.LT.0 .OR. kb.GT.ka ) THEN
217  info = -5
218  ELSE IF( ldab.LT.ka+1 ) THEN
219  info = -7
220  ELSE IF( ldbb.LT.kb+1 ) THEN
221  info = -9
222  ELSE IF( ldx.LT.1 .OR. wantx .AND. ldx.LT.max( 1, n ) ) THEN
223  info = -11
224  END IF
225  IF( info.NE.0 ) THEN
226  CALL xerbla( 'DSBGST', -info )
227  return
228  END IF
229 *
230 * Quick return if possible
231 *
232  IF( n.EQ.0 )
233  $ return
234 *
235  inca = ldab*ka1
236 *
237 * Initialize X to the unit matrix, if needed
238 *
239  IF( wantx )
240  $ CALL dlaset( 'Full', n, n, zero, one, x, ldx )
241 *
242 * Set M to the splitting point m. It must be the same value as is
243 * used in DPBSTF. The chosen value allows the arrays WORK and RWORK
244 * to be of dimension (N).
245 *
246  m = ( n+kb ) / 2
247 *
248 * The routine works in two phases, corresponding to the two halves
249 * of the split Cholesky factorization of B as S**T*S where
250 *
251 * S = ( U )
252 * ( M L )
253 *
254 * with U upper triangular of order m, and L lower triangular of
255 * order n-m. S has the same bandwidth as B.
256 *
257 * S is treated as a product of elementary matrices:
258 *
259 * S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n)
260 *
261 * where S(i) is determined by the i-th row of S.
262 *
263 * In phase 1, the index i takes the values n, n-1, ... , m+1;
264 * in phase 2, it takes the values 1, 2, ... , m.
265 *
266 * For each value of i, the current matrix A is updated by forming
267 * inv(S(i))**T*A*inv(S(i)). This creates a triangular bulge outside
268 * the band of A. The bulge is then pushed down toward the bottom of
269 * A in phase 1, and up toward the top of A in phase 2, by applying
270 * plane rotations.
271 *
272 * There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
273 * of them are linearly independent, so annihilating a bulge requires
274 * only 2*kb-1 plane rotations. The rotations are divided into a 1st
275 * set of kb-1 rotations, and a 2nd set of kb rotations.
276 *
277 * Wherever possible, rotations are generated and applied in vector
278 * operations of length NR between the indices J1 and J2 (sometimes
279 * replaced by modified values NRT, J1T or J2T).
280 *
281 * The cosines and sines of the rotations are stored in the array
282 * WORK. The cosines of the 1st set of rotations are stored in
283 * elements n+2:n+m-kb-1 and the sines of the 1st set in elements
284 * 2:m-kb-1; the cosines of the 2nd set are stored in elements
285 * n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n.
286 *
287 * The bulges are not formed explicitly; nonzero elements outside the
288 * band are created only when they are required for generating new
289 * rotations; they are stored in the array WORK, in positions where
290 * they are later overwritten by the sines of the rotations which
291 * annihilate them.
292 *
293 * **************************** Phase 1 *****************************
294 *
295 * The logical structure of this phase is:
296 *
297 * UPDATE = .TRUE.
298 * DO I = N, M + 1, -1
299 * use S(i) to update A and create a new bulge
300 * apply rotations to push all bulges KA positions downward
301 * END DO
302 * UPDATE = .FALSE.
303 * DO I = M + KA + 1, N - 1
304 * apply rotations to push all bulges KA positions downward
305 * END DO
306 *
307 * To avoid duplicating code, the two loops are merged.
308 *
309  update = .true.
310  i = n + 1
311  10 continue
312  IF( update ) THEN
313  i = i - 1
314  kbt = min( kb, i-1 )
315  i0 = i - 1
316  i1 = min( n, i+ka )
317  i2 = i - kbt + ka1
318  IF( i.LT.m+1 ) THEN
319  update = .false.
320  i = i + 1
321  i0 = m
322  IF( ka.EQ.0 )
323  $ go to 480
324  go to 10
325  END IF
326  ELSE
327  i = i + ka
328  IF( i.GT.n-1 )
329  $ go to 480
330  END IF
331 *
332  IF( upper ) THEN
333 *
334 * Transform A, working with the upper triangle
335 *
336  IF( update ) THEN
337 *
338 * Form inv(S(i))**T * A * inv(S(i))
339 *
340  bii = bb( kb1, i )
341  DO 20 j = i, i1
342  ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
343  20 continue
344  DO 30 j = max( 1, i-ka ), i
345  ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
346  30 continue
347  DO 60 k = i - kbt, i - 1
348  DO 40 j = i - kbt, k
349  ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
350  $ bb( j-i+kb1, i )*ab( k-i+ka1, i ) -
351  $ bb( k-i+kb1, i )*ab( j-i+ka1, i ) +
352  $ ab( ka1, i )*bb( j-i+kb1, i )*
353  $ bb( k-i+kb1, i )
354  40 continue
355  DO 50 j = max( 1, i-ka ), i - kbt - 1
356  ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
357  $ bb( k-i+kb1, i )*ab( j-i+ka1, i )
358  50 continue
359  60 continue
360  DO 80 j = i, i1
361  DO 70 k = max( j-ka, i-kbt ), i - 1
362  ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
363  $ bb( k-i+kb1, i )*ab( i-j+ka1, j )
364  70 continue
365  80 continue
366 *
367  IF( wantx ) THEN
368 *
369 * post-multiply X by inv(S(i))
370 *
371  CALL dscal( n-m, one / bii, x( m+1, i ), 1 )
372  IF( kbt.GT.0 )
373  $ CALL dger( n-m, kbt, -one, x( m+1, i ), 1,
374  $ bb( kb1-kbt, i ), 1, x( m+1, i-kbt ), ldx )
375  END IF
376 *
377 * store a(i,i1) in RA1 for use in next loop over K
378 *
379  ra1 = ab( i-i1+ka1, i1 )
380  END IF
381 *
382 * Generate and apply vectors of rotations to chase all the
383 * existing bulges KA positions down toward the bottom of the
384 * band
385 *
386  DO 130 k = 1, kb - 1
387  IF( update ) THEN
388 *
389 * Determine the rotations which would annihilate the bulge
390 * which has in theory just been created
391 *
392  IF( i-k+ka.LT.n .AND. i-k.GT.1 ) THEN
393 *
394 * generate rotation to annihilate a(i,i-k+ka+1)
395 *
396  CALL dlartg( ab( k+1, i-k+ka ), ra1,
397  $ work( n+i-k+ka-m ), work( i-k+ka-m ),
398  $ ra )
399 *
400 * create nonzero element a(i-k,i-k+ka+1) outside the
401 * band and store it in WORK(i-k)
402 *
403  t = -bb( kb1-k, i )*ra1
404  work( i-k ) = work( n+i-k+ka-m )*t -
405  $ work( i-k+ka-m )*ab( 1, i-k+ka )
406  ab( 1, i-k+ka ) = work( i-k+ka-m )*t +
407  $ work( n+i-k+ka-m )*ab( 1, i-k+ka )
408  ra1 = ra
409  END IF
410  END IF
411  j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
412  nr = ( n-j2+ka ) / ka1
413  j1 = j2 + ( nr-1 )*ka1
414  IF( update ) THEN
415  j2t = max( j2, i+2*ka-k+1 )
416  ELSE
417  j2t = j2
418  END IF
419  nrt = ( n-j2t+ka ) / ka1
420  DO 90 j = j2t, j1, ka1
421 *
422 * create nonzero element a(j-ka,j+1) outside the band
423 * and store it in WORK(j-m)
424 *
425  work( j-m ) = work( j-m )*ab( 1, j+1 )
426  ab( 1, j+1 ) = work( n+j-m )*ab( 1, j+1 )
427  90 continue
428 *
429 * generate rotations in 1st set to annihilate elements which
430 * have been created outside the band
431 *
432  IF( nrt.GT.0 )
433  $ CALL dlargv( nrt, ab( 1, j2t ), inca, work( j2t-m ), ka1,
434  $ work( n+j2t-m ), ka1 )
435  IF( nr.GT.0 ) THEN
436 *
437 * apply rotations in 1st set from the right
438 *
439  DO 100 l = 1, ka - 1
440  CALL dlartv( nr, ab( ka1-l, j2 ), inca,
441  $ ab( ka-l, j2+1 ), inca, work( n+j2-m ),
442  $ work( j2-m ), ka1 )
443  100 continue
444 *
445 * apply rotations in 1st set from both sides to diagonal
446 * blocks
447 *
448  CALL dlar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
449  $ ab( ka, j2+1 ), inca, work( n+j2-m ),
450  $ work( j2-m ), ka1 )
451 *
452  END IF
453 *
454 * start applying rotations in 1st set from the left
455 *
456  DO 110 l = ka - 1, kb - k + 1, -1
457  nrt = ( n-j2+l ) / ka1
458  IF( nrt.GT.0 )
459  $ CALL dlartv( nrt, ab( l, j2+ka1-l ), inca,
460  $ ab( l+1, j2+ka1-l ), inca,
461  $ work( n+j2-m ), work( j2-m ), ka1 )
462  110 continue
463 *
464  IF( wantx ) THEN
465 *
466 * post-multiply X by product of rotations in 1st set
467 *
468  DO 120 j = j2, j1, ka1
469  CALL drot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
470  $ work( n+j-m ), work( j-m ) )
471  120 continue
472  END IF
473  130 continue
474 *
475  IF( update ) THEN
476  IF( i2.LE.n .AND. kbt.GT.0 ) THEN
477 *
478 * create nonzero element a(i-kbt,i-kbt+ka+1) outside the
479 * band and store it in WORK(i-kbt)
480 *
481  work( i-kbt ) = -bb( kb1-kbt, i )*ra1
482  END IF
483  END IF
484 *
485  DO 170 k = kb, 1, -1
486  IF( update ) THEN
487  j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
488  ELSE
489  j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
490  END IF
491 *
492 * finish applying rotations in 2nd set from the left
493 *
494  DO 140 l = kb - k, 1, -1
495  nrt = ( n-j2+ka+l ) / ka1
496  IF( nrt.GT.0 )
497  $ CALL dlartv( nrt, ab( l, j2-l+1 ), inca,
498  $ ab( l+1, j2-l+1 ), inca, work( n+j2-ka ),
499  $ work( j2-ka ), ka1 )
500  140 continue
501  nr = ( n-j2+ka ) / ka1
502  j1 = j2 + ( nr-1 )*ka1
503  DO 150 j = j1, j2, -ka1
504  work( j ) = work( j-ka )
505  work( n+j ) = work( n+j-ka )
506  150 continue
507  DO 160 j = j2, j1, ka1
508 *
509 * create nonzero element a(j-ka,j+1) outside the band
510 * and store it in WORK(j)
511 *
512  work( j ) = work( j )*ab( 1, j+1 )
513  ab( 1, j+1 ) = work( n+j )*ab( 1, j+1 )
514  160 continue
515  IF( update ) THEN
516  IF( i-k.LT.n-ka .AND. k.LE.kbt )
517  $ work( i-k+ka ) = work( i-k )
518  END IF
519  170 continue
520 *
521  DO 210 k = kb, 1, -1
522  j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
523  nr = ( n-j2+ka ) / ka1
524  j1 = j2 + ( nr-1 )*ka1
525  IF( nr.GT.0 ) THEN
526 *
527 * generate rotations in 2nd set to annihilate elements
528 * which have been created outside the band
529 *
530  CALL dlargv( nr, ab( 1, j2 ), inca, work( j2 ), ka1,
531  $ work( n+j2 ), ka1 )
532 *
533 * apply rotations in 2nd set from the right
534 *
535  DO 180 l = 1, ka - 1
536  CALL dlartv( nr, ab( ka1-l, j2 ), inca,
537  $ ab( ka-l, j2+1 ), inca, work( n+j2 ),
538  $ work( j2 ), ka1 )
539  180 continue
540 *
541 * apply rotations in 2nd set from both sides to diagonal
542 * blocks
543 *
544  CALL dlar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
545  $ ab( ka, j2+1 ), inca, work( n+j2 ),
546  $ work( j2 ), ka1 )
547 *
548  END IF
549 *
550 * start applying rotations in 2nd set from the left
551 *
552  DO 190 l = ka - 1, kb - k + 1, -1
553  nrt = ( n-j2+l ) / ka1
554  IF( nrt.GT.0 )
555  $ CALL dlartv( nrt, ab( l, j2+ka1-l ), inca,
556  $ ab( l+1, j2+ka1-l ), inca, work( n+j2 ),
557  $ work( j2 ), ka1 )
558  190 continue
559 *
560  IF( wantx ) THEN
561 *
562 * post-multiply X by product of rotations in 2nd set
563 *
564  DO 200 j = j2, j1, ka1
565  CALL drot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
566  $ work( n+j ), work( j ) )
567  200 continue
568  END IF
569  210 continue
570 *
571  DO 230 k = 1, kb - 1
572  j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
573 *
574 * finish applying rotations in 1st set from the left
575 *
576  DO 220 l = kb - k, 1, -1
577  nrt = ( n-j2+l ) / ka1
578  IF( nrt.GT.0 )
579  $ CALL dlartv( nrt, ab( l, j2+ka1-l ), inca,
580  $ ab( l+1, j2+ka1-l ), inca,
581  $ work( n+j2-m ), work( j2-m ), ka1 )
582  220 continue
583  230 continue
584 *
585  IF( kb.GT.1 ) THEN
586  DO 240 j = n - 1, i - kb + 2*ka + 1, -1
587  work( n+j-m ) = work( n+j-ka-m )
588  work( j-m ) = work( j-ka-m )
589  240 continue
590  END IF
591 *
592  ELSE
593 *
594 * Transform A, working with the lower triangle
595 *
596  IF( update ) THEN
597 *
598 * Form inv(S(i))**T * A * inv(S(i))
599 *
600  bii = bb( 1, i )
601  DO 250 j = i, i1
602  ab( j-i+1, i ) = ab( j-i+1, i ) / bii
603  250 continue
604  DO 260 j = max( 1, i-ka ), i
605  ab( i-j+1, j ) = ab( i-j+1, j ) / bii
606  260 continue
607  DO 290 k = i - kbt, i - 1
608  DO 270 j = i - kbt, k
609  ab( k-j+1, j ) = ab( k-j+1, j ) -
610  $ bb( i-j+1, j )*ab( i-k+1, k ) -
611  $ bb( i-k+1, k )*ab( i-j+1, j ) +
612  $ ab( 1, i )*bb( i-j+1, j )*
613  $ bb( i-k+1, k )
614  270 continue
615  DO 280 j = max( 1, i-ka ), i - kbt - 1
616  ab( k-j+1, j ) = ab( k-j+1, j ) -
617  $ bb( i-k+1, k )*ab( i-j+1, j )
618  280 continue
619  290 continue
620  DO 310 j = i, i1
621  DO 300 k = max( j-ka, i-kbt ), i - 1
622  ab( j-k+1, k ) = ab( j-k+1, k ) -
623  $ bb( i-k+1, k )*ab( j-i+1, i )
624  300 continue
625  310 continue
626 *
627  IF( wantx ) THEN
628 *
629 * post-multiply X by inv(S(i))
630 *
631  CALL dscal( n-m, one / bii, x( m+1, i ), 1 )
632  IF( kbt.GT.0 )
633  $ CALL dger( n-m, kbt, -one, x( m+1, i ), 1,
634  $ bb( kbt+1, i-kbt ), ldbb-1,
635  $ x( m+1, i-kbt ), ldx )
636  END IF
637 *
638 * store a(i1,i) in RA1 for use in next loop over K
639 *
640  ra1 = ab( i1-i+1, i )
641  END IF
642 *
643 * Generate and apply vectors of rotations to chase all the
644 * existing bulges KA positions down toward the bottom of the
645 * band
646 *
647  DO 360 k = 1, kb - 1
648  IF( update ) THEN
649 *
650 * Determine the rotations which would annihilate the bulge
651 * which has in theory just been created
652 *
653  IF( i-k+ka.LT.n .AND. i-k.GT.1 ) THEN
654 *
655 * generate rotation to annihilate a(i-k+ka+1,i)
656 *
657  CALL dlartg( ab( ka1-k, i ), ra1, work( n+i-k+ka-m ),
658  $ work( i-k+ka-m ), ra )
659 *
660 * create nonzero element a(i-k+ka+1,i-k) outside the
661 * band and store it in WORK(i-k)
662 *
663  t = -bb( k+1, i-k )*ra1
664  work( i-k ) = work( n+i-k+ka-m )*t -
665  $ work( i-k+ka-m )*ab( ka1, i-k )
666  ab( ka1, i-k ) = work( i-k+ka-m )*t +
667  $ work( n+i-k+ka-m )*ab( ka1, i-k )
668  ra1 = ra
669  END IF
670  END IF
671  j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
672  nr = ( n-j2+ka ) / ka1
673  j1 = j2 + ( nr-1 )*ka1
674  IF( update ) THEN
675  j2t = max( j2, i+2*ka-k+1 )
676  ELSE
677  j2t = j2
678  END IF
679  nrt = ( n-j2t+ka ) / ka1
680  DO 320 j = j2t, j1, ka1
681 *
682 * create nonzero element a(j+1,j-ka) outside the band
683 * and store it in WORK(j-m)
684 *
685  work( j-m ) = work( j-m )*ab( ka1, j-ka+1 )
686  ab( ka1, j-ka+1 ) = work( n+j-m )*ab( ka1, j-ka+1 )
687  320 continue
688 *
689 * generate rotations in 1st set to annihilate elements which
690 * have been created outside the band
691 *
692  IF( nrt.GT.0 )
693  $ CALL dlargv( nrt, ab( ka1, j2t-ka ), inca, work( j2t-m ),
694  $ ka1, work( n+j2t-m ), ka1 )
695  IF( nr.GT.0 ) THEN
696 *
697 * apply rotations in 1st set from the left
698 *
699  DO 330 l = 1, ka - 1
700  CALL dlartv( nr, ab( l+1, j2-l ), inca,
701  $ ab( l+2, j2-l ), inca, work( n+j2-m ),
702  $ work( j2-m ), ka1 )
703  330 continue
704 *
705 * apply rotations in 1st set from both sides to diagonal
706 * blocks
707 *
708  CALL dlar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2, j2 ),
709  $ inca, work( n+j2-m ), work( j2-m ), ka1 )
710 *
711  END IF
712 *
713 * start applying rotations in 1st set from the right
714 *
715  DO 340 l = ka - 1, kb - k + 1, -1
716  nrt = ( n-j2+l ) / ka1
717  IF( nrt.GT.0 )
718  $ CALL dlartv( nrt, ab( ka1-l+1, j2 ), inca,
719  $ ab( ka1-l, j2+1 ), inca, work( n+j2-m ),
720  $ work( j2-m ), ka1 )
721  340 continue
722 *
723  IF( wantx ) THEN
724 *
725 * post-multiply X by product of rotations in 1st set
726 *
727  DO 350 j = j2, j1, ka1
728  CALL drot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
729  $ work( n+j-m ), work( j-m ) )
730  350 continue
731  END IF
732  360 continue
733 *
734  IF( update ) THEN
735  IF( i2.LE.n .AND. kbt.GT.0 ) THEN
736 *
737 * create nonzero element a(i-kbt+ka+1,i-kbt) outside the
738 * band and store it in WORK(i-kbt)
739 *
740  work( i-kbt ) = -bb( kbt+1, i-kbt )*ra1
741  END IF
742  END IF
743 *
744  DO 400 k = kb, 1, -1
745  IF( update ) THEN
746  j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
747  ELSE
748  j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
749  END IF
750 *
751 * finish applying rotations in 2nd set from the right
752 *
753  DO 370 l = kb - k, 1, -1
754  nrt = ( n-j2+ka+l ) / ka1
755  IF( nrt.GT.0 )
756  $ CALL dlartv( nrt, ab( ka1-l+1, j2-ka ), inca,
757  $ ab( ka1-l, j2-ka+1 ), inca,
758  $ work( n+j2-ka ), work( j2-ka ), ka1 )
759  370 continue
760  nr = ( n-j2+ka ) / ka1
761  j1 = j2 + ( nr-1 )*ka1
762  DO 380 j = j1, j2, -ka1
763  work( j ) = work( j-ka )
764  work( n+j ) = work( n+j-ka )
765  380 continue
766  DO 390 j = j2, j1, ka1
767 *
768 * create nonzero element a(j+1,j-ka) outside the band
769 * and store it in WORK(j)
770 *
771  work( j ) = work( j )*ab( ka1, j-ka+1 )
772  ab( ka1, j-ka+1 ) = work( n+j )*ab( ka1, j-ka+1 )
773  390 continue
774  IF( update ) THEN
775  IF( i-k.LT.n-ka .AND. k.LE.kbt )
776  $ work( i-k+ka ) = work( i-k )
777  END IF
778  400 continue
779 *
780  DO 440 k = kb, 1, -1
781  j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
782  nr = ( n-j2+ka ) / ka1
783  j1 = j2 + ( nr-1 )*ka1
784  IF( nr.GT.0 ) THEN
785 *
786 * generate rotations in 2nd set to annihilate elements
787 * which have been created outside the band
788 *
789  CALL dlargv( nr, ab( ka1, j2-ka ), inca, work( j2 ), ka1,
790  $ work( n+j2 ), ka1 )
791 *
792 * apply rotations in 2nd set from the left
793 *
794  DO 410 l = 1, ka - 1
795  CALL dlartv( nr, ab( l+1, j2-l ), inca,
796  $ ab( l+2, j2-l ), inca, work( n+j2 ),
797  $ work( j2 ), ka1 )
798  410 continue
799 *
800 * apply rotations in 2nd set from both sides to diagonal
801 * blocks
802 *
803  CALL dlar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2, j2 ),
804  $ inca, work( n+j2 ), work( j2 ), ka1 )
805 *
806  END IF
807 *
808 * start applying rotations in 2nd set from the right
809 *
810  DO 420 l = ka - 1, kb - k + 1, -1
811  nrt = ( n-j2+l ) / ka1
812  IF( nrt.GT.0 )
813  $ CALL dlartv( nrt, ab( ka1-l+1, j2 ), inca,
814  $ ab( ka1-l, j2+1 ), inca, work( n+j2 ),
815  $ work( j2 ), ka1 )
816  420 continue
817 *
818  IF( wantx ) THEN
819 *
820 * post-multiply X by product of rotations in 2nd set
821 *
822  DO 430 j = j2, j1, ka1
823  CALL drot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
824  $ work( n+j ), work( j ) )
825  430 continue
826  END IF
827  440 continue
828 *
829  DO 460 k = 1, kb - 1
830  j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
831 *
832 * finish applying rotations in 1st set from the right
833 *
834  DO 450 l = kb - k, 1, -1
835  nrt = ( n-j2+l ) / ka1
836  IF( nrt.GT.0 )
837  $ CALL dlartv( nrt, ab( ka1-l+1, j2 ), inca,
838  $ ab( ka1-l, j2+1 ), inca, work( n+j2-m ),
839  $ work( j2-m ), ka1 )
840  450 continue
841  460 continue
842 *
843  IF( kb.GT.1 ) THEN
844  DO 470 j = n - 1, i - kb + 2*ka + 1, -1
845  work( n+j-m ) = work( n+j-ka-m )
846  work( j-m ) = work( j-ka-m )
847  470 continue
848  END IF
849 *
850  END IF
851 *
852  go to 10
853 *
854  480 continue
855 *
856 * **************************** Phase 2 *****************************
857 *
858 * The logical structure of this phase is:
859 *
860 * UPDATE = .TRUE.
861 * DO I = 1, M
862 * use S(i) to update A and create a new bulge
863 * apply rotations to push all bulges KA positions upward
864 * END DO
865 * UPDATE = .FALSE.
866 * DO I = M - KA - 1, 2, -1
867 * apply rotations to push all bulges KA positions upward
868 * END DO
869 *
870 * To avoid duplicating code, the two loops are merged.
871 *
872  update = .true.
873  i = 0
874  490 continue
875  IF( update ) THEN
876  i = i + 1
877  kbt = min( kb, m-i )
878  i0 = i + 1
879  i1 = max( 1, i-ka )
880  i2 = i + kbt - ka1
881  IF( i.GT.m ) THEN
882  update = .false.
883  i = i - 1
884  i0 = m + 1
885  IF( ka.EQ.0 )
886  $ return
887  go to 490
888  END IF
889  ELSE
890  i = i - ka
891  IF( i.LT.2 )
892  $ return
893  END IF
894 *
895  IF( i.LT.m-kbt ) THEN
896  nx = m
897  ELSE
898  nx = n
899  END IF
900 *
901  IF( upper ) THEN
902 *
903 * Transform A, working with the upper triangle
904 *
905  IF( update ) THEN
906 *
907 * Form inv(S(i))**T * A * inv(S(i))
908 *
909  bii = bb( kb1, i )
910  DO 500 j = i1, i
911  ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
912  500 continue
913  DO 510 j = i, min( n, i+ka )
914  ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
915  510 continue
916  DO 540 k = i + 1, i + kbt
917  DO 520 j = k, i + kbt
918  ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
919  $ bb( i-j+kb1, j )*ab( i-k+ka1, k ) -
920  $ bb( i-k+kb1, k )*ab( i-j+ka1, j ) +
921  $ ab( ka1, i )*bb( i-j+kb1, j )*
922  $ bb( i-k+kb1, k )
923  520 continue
924  DO 530 j = i + kbt + 1, min( n, i+ka )
925  ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
926  $ bb( i-k+kb1, k )*ab( i-j+ka1, j )
927  530 continue
928  540 continue
929  DO 560 j = i1, i
930  DO 550 k = i + 1, min( j+ka, i+kbt )
931  ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
932  $ bb( i-k+kb1, k )*ab( j-i+ka1, i )
933  550 continue
934  560 continue
935 *
936  IF( wantx ) THEN
937 *
938 * post-multiply X by inv(S(i))
939 *
940  CALL dscal( nx, one / bii, x( 1, i ), 1 )
941  IF( kbt.GT.0 )
942  $ CALL dger( nx, kbt, -one, x( 1, i ), 1, bb( kb, i+1 ),
943  $ ldbb-1, x( 1, i+1 ), ldx )
944  END IF
945 *
946 * store a(i1,i) in RA1 for use in next loop over K
947 *
948  ra1 = ab( i1-i+ka1, i )
949  END IF
950 *
951 * Generate and apply vectors of rotations to chase all the
952 * existing bulges KA positions up toward the top of the band
953 *
954  DO 610 k = 1, kb - 1
955  IF( update ) THEN
956 *
957 * Determine the rotations which would annihilate the bulge
958 * which has in theory just been created
959 *
960  IF( i+k-ka1.GT.0 .AND. i+k.LT.m ) THEN
961 *
962 * generate rotation to annihilate a(i+k-ka-1,i)
963 *
964  CALL dlartg( ab( k+1, i ), ra1, work( n+i+k-ka ),
965  $ work( i+k-ka ), ra )
966 *
967 * create nonzero element a(i+k-ka-1,i+k) outside the
968 * band and store it in WORK(m-kb+i+k)
969 *
970  t = -bb( kb1-k, i+k )*ra1
971  work( m-kb+i+k ) = work( n+i+k-ka )*t -
972  $ work( i+k-ka )*ab( 1, i+k )
973  ab( 1, i+k ) = work( i+k-ka )*t +
974  $ work( n+i+k-ka )*ab( 1, i+k )
975  ra1 = ra
976  END IF
977  END IF
978  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
979  nr = ( j2+ka-1 ) / ka1
980  j1 = j2 - ( nr-1 )*ka1
981  IF( update ) THEN
982  j2t = min( j2, i-2*ka+k-1 )
983  ELSE
984  j2t = j2
985  END IF
986  nrt = ( j2t+ka-1 ) / ka1
987  DO 570 j = j1, j2t, ka1
988 *
989 * create nonzero element a(j-1,j+ka) outside the band
990 * and store it in WORK(j)
991 *
992  work( j ) = work( j )*ab( 1, j+ka-1 )
993  ab( 1, j+ka-1 ) = work( n+j )*ab( 1, j+ka-1 )
994  570 continue
995 *
996 * generate rotations in 1st set to annihilate elements which
997 * have been created outside the band
998 *
999  IF( nrt.GT.0 )
1000  $ CALL dlargv( nrt, ab( 1, j1+ka ), inca, work( j1 ), ka1,
1001  $ work( n+j1 ), ka1 )
1002  IF( nr.GT.0 ) THEN
1003 *
1004 * apply rotations in 1st set from the left
1005 *
1006  DO 580 l = 1, ka - 1
1007  CALL dlartv( nr, ab( ka1-l, j1+l ), inca,
1008  $ ab( ka-l, j1+l ), inca, work( n+j1 ),
1009  $ work( j1 ), ka1 )
1010  580 continue
1011 *
1012 * apply rotations in 1st set from both sides to diagonal
1013 * blocks
1014 *
1015  CALL dlar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1016  $ ab( ka, j1 ), inca, work( n+j1 ),
1017  $ work( j1 ), ka1 )
1018 *
1019  END IF
1020 *
1021 * start applying rotations in 1st set from the right
1022 *
1023  DO 590 l = ka - 1, kb - k + 1, -1
1024  nrt = ( j2+l-1 ) / ka1
1025  j1t = j2 - ( nrt-1 )*ka1
1026  IF( nrt.GT.0 )
1027  $ CALL dlartv( nrt, ab( l, j1t ), inca,
1028  $ ab( l+1, j1t-1 ), inca, work( n+j1t ),
1029  $ work( j1t ), ka1 )
1030  590 continue
1031 *
1032  IF( wantx ) THEN
1033 *
1034 * post-multiply X by product of rotations in 1st set
1035 *
1036  DO 600 j = j1, j2, ka1
1037  CALL drot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1038  $ work( n+j ), work( j ) )
1039  600 continue
1040  END IF
1041  610 continue
1042 *
1043  IF( update ) THEN
1044  IF( i2.GT.0 .AND. kbt.GT.0 ) THEN
1045 *
1046 * create nonzero element a(i+kbt-ka-1,i+kbt) outside the
1047 * band and store it in WORK(m-kb+i+kbt)
1048 *
1049  work( m-kb+i+kbt ) = -bb( kb1-kbt, i+kbt )*ra1
1050  END IF
1051  END IF
1052 *
1053  DO 650 k = kb, 1, -1
1054  IF( update ) THEN
1055  j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1056  ELSE
1057  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1058  END IF
1059 *
1060 * finish applying rotations in 2nd set from the right
1061 *
1062  DO 620 l = kb - k, 1, -1
1063  nrt = ( j2+ka+l-1 ) / ka1
1064  j1t = j2 - ( nrt-1 )*ka1
1065  IF( nrt.GT.0 )
1066  $ CALL dlartv( nrt, ab( l, j1t+ka ), inca,
1067  $ ab( l+1, j1t+ka-1 ), inca,
1068  $ work( n+m-kb+j1t+ka ),
1069  $ work( m-kb+j1t+ka ), ka1 )
1070  620 continue
1071  nr = ( j2+ka-1 ) / ka1
1072  j1 = j2 - ( nr-1 )*ka1
1073  DO 630 j = j1, j2, ka1
1074  work( m-kb+j ) = work( m-kb+j+ka )
1075  work( n+m-kb+j ) = work( n+m-kb+j+ka )
1076  630 continue
1077  DO 640 j = j1, j2, ka1
1078 *
1079 * create nonzero element a(j-1,j+ka) outside the band
1080 * and store it in WORK(m-kb+j)
1081 *
1082  work( m-kb+j ) = work( m-kb+j )*ab( 1, j+ka-1 )
1083  ab( 1, j+ka-1 ) = work( n+m-kb+j )*ab( 1, j+ka-1 )
1084  640 continue
1085  IF( update ) THEN
1086  IF( i+k.GT.ka1 .AND. k.LE.kbt )
1087  $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1088  END IF
1089  650 continue
1090 *
1091  DO 690 k = kb, 1, -1
1092  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1093  nr = ( j2+ka-1 ) / ka1
1094  j1 = j2 - ( nr-1 )*ka1
1095  IF( nr.GT.0 ) THEN
1096 *
1097 * generate rotations in 2nd set to annihilate elements
1098 * which have been created outside the band
1099 *
1100  CALL dlargv( nr, ab( 1, j1+ka ), inca, work( m-kb+j1 ),
1101  $ ka1, work( n+m-kb+j1 ), ka1 )
1102 *
1103 * apply rotations in 2nd set from the left
1104 *
1105  DO 660 l = 1, ka - 1
1106  CALL dlartv( nr, ab( ka1-l, j1+l ), inca,
1107  $ ab( ka-l, j1+l ), inca,
1108  $ work( n+m-kb+j1 ), work( m-kb+j1 ), ka1 )
1109  660 continue
1110 *
1111 * apply rotations in 2nd set from both sides to diagonal
1112 * blocks
1113 *
1114  CALL dlar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1115  $ ab( ka, j1 ), inca, work( n+m-kb+j1 ),
1116  $ work( m-kb+j1 ), ka1 )
1117 *
1118  END IF
1119 *
1120 * start applying rotations in 2nd set from the right
1121 *
1122  DO 670 l = ka - 1, kb - k + 1, -1
1123  nrt = ( j2+l-1 ) / ka1
1124  j1t = j2 - ( nrt-1 )*ka1
1125  IF( nrt.GT.0 )
1126  $ CALL dlartv( nrt, ab( l, j1t ), inca,
1127  $ ab( l+1, j1t-1 ), inca,
1128  $ work( n+m-kb+j1t ), work( m-kb+j1t ),
1129  $ ka1 )
1130  670 continue
1131 *
1132  IF( wantx ) THEN
1133 *
1134 * post-multiply X by product of rotations in 2nd set
1135 *
1136  DO 680 j = j1, j2, ka1
1137  CALL drot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1138  $ work( n+m-kb+j ), work( m-kb+j ) )
1139  680 continue
1140  END IF
1141  690 continue
1142 *
1143  DO 710 k = 1, kb - 1
1144  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1145 *
1146 * finish applying rotations in 1st set from the right
1147 *
1148  DO 700 l = kb - k, 1, -1
1149  nrt = ( j2+l-1 ) / ka1
1150  j1t = j2 - ( nrt-1 )*ka1
1151  IF( nrt.GT.0 )
1152  $ CALL dlartv( nrt, ab( l, j1t ), inca,
1153  $ ab( l+1, j1t-1 ), inca, work( n+j1t ),
1154  $ work( j1t ), ka1 )
1155  700 continue
1156  710 continue
1157 *
1158  IF( kb.GT.1 ) THEN
1159  DO 720 j = 2, min( i+kb, m ) - 2*ka - 1
1160  work( n+j ) = work( n+j+ka )
1161  work( j ) = work( j+ka )
1162  720 continue
1163  END IF
1164 *
1165  ELSE
1166 *
1167 * Transform A, working with the lower triangle
1168 *
1169  IF( update ) THEN
1170 *
1171 * Form inv(S(i))**T * A * inv(S(i))
1172 *
1173  bii = bb( 1, i )
1174  DO 730 j = i1, i
1175  ab( i-j+1, j ) = ab( i-j+1, j ) / bii
1176  730 continue
1177  DO 740 j = i, min( n, i+ka )
1178  ab( j-i+1, i ) = ab( j-i+1, i ) / bii
1179  740 continue
1180  DO 770 k = i + 1, i + kbt
1181  DO 750 j = k, i + kbt
1182  ab( j-k+1, k ) = ab( j-k+1, k ) -
1183  $ bb( j-i+1, i )*ab( k-i+1, i ) -
1184  $ bb( k-i+1, i )*ab( j-i+1, i ) +
1185  $ ab( 1, i )*bb( j-i+1, i )*
1186  $ bb( k-i+1, i )
1187  750 continue
1188  DO 760 j = i + kbt + 1, min( n, i+ka )
1189  ab( j-k+1, k ) = ab( j-k+1, k ) -
1190  $ bb( k-i+1, i )*ab( j-i+1, i )
1191  760 continue
1192  770 continue
1193  DO 790 j = i1, i
1194  DO 780 k = i + 1, min( j+ka, i+kbt )
1195  ab( k-j+1, j ) = ab( k-j+1, j ) -
1196  $ bb( k-i+1, i )*ab( i-j+1, j )
1197  780 continue
1198  790 continue
1199 *
1200  IF( wantx ) THEN
1201 *
1202 * post-multiply X by inv(S(i))
1203 *
1204  CALL dscal( nx, one / bii, x( 1, i ), 1 )
1205  IF( kbt.GT.0 )
1206  $ CALL dger( nx, kbt, -one, x( 1, i ), 1, bb( 2, i ), 1,
1207  $ x( 1, i+1 ), ldx )
1208  END IF
1209 *
1210 * store a(i,i1) in RA1 for use in next loop over K
1211 *
1212  ra1 = ab( i-i1+1, i1 )
1213  END IF
1214 *
1215 * Generate and apply vectors of rotations to chase all the
1216 * existing bulges KA positions up toward the top of the band
1217 *
1218  DO 840 k = 1, kb - 1
1219  IF( update ) THEN
1220 *
1221 * Determine the rotations which would annihilate the bulge
1222 * which has in theory just been created
1223 *
1224  IF( i+k-ka1.GT.0 .AND. i+k.LT.m ) THEN
1225 *
1226 * generate rotation to annihilate a(i,i+k-ka-1)
1227 *
1228  CALL dlartg( ab( ka1-k, i+k-ka ), ra1,
1229  $ work( n+i+k-ka ), work( i+k-ka ), ra )
1230 *
1231 * create nonzero element a(i+k,i+k-ka-1) outside the
1232 * band and store it in WORK(m-kb+i+k)
1233 *
1234  t = -bb( k+1, i )*ra1
1235  work( m-kb+i+k ) = work( n+i+k-ka )*t -
1236  $ work( i+k-ka )*ab( ka1, i+k-ka )
1237  ab( ka1, i+k-ka ) = work( i+k-ka )*t +
1238  $ work( n+i+k-ka )*ab( ka1, i+k-ka )
1239  ra1 = ra
1240  END IF
1241  END IF
1242  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1243  nr = ( j2+ka-1 ) / ka1
1244  j1 = j2 - ( nr-1 )*ka1
1245  IF( update ) THEN
1246  j2t = min( j2, i-2*ka+k-1 )
1247  ELSE
1248  j2t = j2
1249  END IF
1250  nrt = ( j2t+ka-1 ) / ka1
1251  DO 800 j = j1, j2t, ka1
1252 *
1253 * create nonzero element a(j+ka,j-1) outside the band
1254 * and store it in WORK(j)
1255 *
1256  work( j ) = work( j )*ab( ka1, j-1 )
1257  ab( ka1, j-1 ) = work( n+j )*ab( ka1, j-1 )
1258  800 continue
1259 *
1260 * generate rotations in 1st set to annihilate elements which
1261 * have been created outside the band
1262 *
1263  IF( nrt.GT.0 )
1264  $ CALL dlargv( nrt, ab( ka1, j1 ), inca, work( j1 ), ka1,
1265  $ work( n+j1 ), ka1 )
1266  IF( nr.GT.0 ) THEN
1267 *
1268 * apply rotations in 1st set from the right
1269 *
1270  DO 810 l = 1, ka - 1
1271  CALL dlartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),
1272  $ inca, work( n+j1 ), work( j1 ), ka1 )
1273  810 continue
1274 *
1275 * apply rotations in 1st set from both sides to diagonal
1276 * blocks
1277 *
1278  CALL dlar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1279  $ ab( 2, j1-1 ), inca, work( n+j1 ),
1280  $ work( j1 ), ka1 )
1281 *
1282  END IF
1283 *
1284 * start applying rotations in 1st set from the left
1285 *
1286  DO 820 l = ka - 1, kb - k + 1, -1
1287  nrt = ( j2+l-1 ) / ka1
1288  j1t = j2 - ( nrt-1 )*ka1
1289  IF( nrt.GT.0 )
1290  $ CALL dlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1291  $ ab( ka1-l, j1t-ka1+l ), inca,
1292  $ work( n+j1t ), work( j1t ), ka1 )
1293  820 continue
1294 *
1295  IF( wantx ) THEN
1296 *
1297 * post-multiply X by product of rotations in 1st set
1298 *
1299  DO 830 j = j1, j2, ka1
1300  CALL drot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1301  $ work( n+j ), work( j ) )
1302  830 continue
1303  END IF
1304  840 continue
1305 *
1306  IF( update ) THEN
1307  IF( i2.GT.0 .AND. kbt.GT.0 ) THEN
1308 *
1309 * create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1310 * band and store it in WORK(m-kb+i+kbt)
1311 *
1312  work( m-kb+i+kbt ) = -bb( kbt+1, i )*ra1
1313  END IF
1314  END IF
1315 *
1316  DO 880 k = kb, 1, -1
1317  IF( update ) THEN
1318  j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1319  ELSE
1320  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1321  END IF
1322 *
1323 * finish applying rotations in 2nd set from the left
1324 *
1325  DO 850 l = kb - k, 1, -1
1326  nrt = ( j2+ka+l-1 ) / ka1
1327  j1t = j2 - ( nrt-1 )*ka1
1328  IF( nrt.GT.0 )
1329  $ CALL dlartv( nrt, ab( ka1-l+1, j1t+l-1 ), inca,
1330  $ ab( ka1-l, j1t+l-1 ), inca,
1331  $ work( n+m-kb+j1t+ka ),
1332  $ work( m-kb+j1t+ka ), ka1 )
1333  850 continue
1334  nr = ( j2+ka-1 ) / ka1
1335  j1 = j2 - ( nr-1 )*ka1
1336  DO 860 j = j1, j2, ka1
1337  work( m-kb+j ) = work( m-kb+j+ka )
1338  work( n+m-kb+j ) = work( n+m-kb+j+ka )
1339  860 continue
1340  DO 870 j = j1, j2, ka1
1341 *
1342 * create nonzero element a(j+ka,j-1) outside the band
1343 * and store it in WORK(m-kb+j)
1344 *
1345  work( m-kb+j ) = work( m-kb+j )*ab( ka1, j-1 )
1346  ab( ka1, j-1 ) = work( n+m-kb+j )*ab( ka1, j-1 )
1347  870 continue
1348  IF( update ) THEN
1349  IF( i+k.GT.ka1 .AND. k.LE.kbt )
1350  $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1351  END IF
1352  880 continue
1353 *
1354  DO 920 k = kb, 1, -1
1355  j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1356  nr = ( j2+ka-1 ) / ka1
1357  j1 = j2 - ( nr-1 )*ka1
1358  IF( nr.GT.0 ) THEN
1359 *
1360 * generate rotations in 2nd set to annihilate elements
1361 * which have been created outside the band
1362 *
1363  CALL dlargv( nr, ab( ka1, j1 ), inca, work( m-kb+j1 ),
1364  $ ka1, work( n+m-kb+j1 ), ka1 )
1365 *
1366 * apply rotations in 2nd set from the right
1367 *
1368  DO 890 l = 1, ka - 1
1369  CALL dlartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),
1370  $ inca, work( n+m-kb+j1 ), work( m-kb+j1 ),
1371  $ ka1 )
1372  890 continue
1373 *
1374 * apply rotations in 2nd set from both sides to diagonal
1375 * blocks
1376 *
1377  CALL dlar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1378  $ ab( 2, j1-1 ), inca, work( n+m-kb+j1 ),
1379  $ work( m-kb+j1 ), ka1 )
1380 *
1381  END IF
1382 *
1383 * start applying rotations in 2nd set from the left
1384 *
1385  DO 900 l = ka - 1, kb - k + 1, -1
1386  nrt = ( j2+l-1 ) / ka1
1387  j1t = j2 - ( nrt-1 )*ka1
1388  IF( nrt.GT.0 )
1389  $ CALL dlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1390  $ ab( ka1-l, j1t-ka1+l ), inca,
1391  $ work( n+m-kb+j1t ), work( m-kb+j1t ),
1392  $ ka1 )
1393  900 continue
1394 *
1395  IF( wantx ) THEN
1396 *
1397 * post-multiply X by product of rotations in 2nd set
1398 *
1399  DO 910 j = j1, j2, ka1
1400  CALL drot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1401  $ work( n+m-kb+j ), work( m-kb+j ) )
1402  910 continue
1403  END IF
1404  920 continue
1405 *
1406  DO 940 k = 1, kb - 1
1407  j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1408 *
1409 * finish applying rotations in 1st set from the left
1410 *
1411  DO 930 l = kb - k, 1, -1
1412  nrt = ( j2+l-1 ) / ka1
1413  j1t = j2 - ( nrt-1 )*ka1
1414  IF( nrt.GT.0 )
1415  $ CALL dlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1416  $ ab( ka1-l, j1t-ka1+l ), inca,
1417  $ work( n+j1t ), work( j1t ), ka1 )
1418  930 continue
1419  940 continue
1420 *
1421  IF( kb.GT.1 ) THEN
1422  DO 950 j = 2, min( i+kb, m ) - 2*ka - 1
1423  work( n+j ) = work( n+j+ka )
1424  work( j ) = work( j+ka )
1425  950 continue
1426  END IF
1427 *
1428  END IF
1429 *
1430  go to 490
1431 *
1432 * End of DSBGST
1433 *
1434  END