LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ ssbgst()

subroutine ssbgst ( character  VECT,
character  UPLO,
integer  N,
integer  KA,
integer  KB,
real, dimension( ldab, * )  AB,
integer  LDAB,
real, dimension( ldbb, * )  BB,
integer  LDBB,
real, dimension( ldx, * )  X,
integer  LDX,
real, dimension( * )  WORK,
integer  INFO 
)

SSBGST

Download SSBGST + dependencies [TGZ] [ZIP] [TXT]

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.
Parameters
[in]VECT
          VECT is CHARACTER*1
          = 'N':  do not form the transformation matrix X;
          = 'V':  form X.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrices A and B.  N >= 0.
[in]KA
          KA is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
[in]KB
          KB is INTEGER
          The number of superdiagonals of the matrix B if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KA >= KB >= 0.
[in,out]AB
          AB is 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.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KA+1.
[in]BB
          BB is 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.
[in]LDBB
          LDBB is INTEGER
          The leading dimension of the array BB.  LDBB >= KB+1.
[out]X
          X is REAL array, dimension (LDX,N)
          If VECT = 'V', the n-by-n matrix X.
          If VECT = 'N', the array X is not referenced.
[in]LDX
          LDX is INTEGER
          The leading dimension of the array X.
          LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
[out]WORK
          WORK is REAL array, dimension (2*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 157 of file ssbgst.f.

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