LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> Download DSBGST + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbgst.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbgst.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbgst.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DSBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
20* LDX, WORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER UPLO, VECT
24* INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N
25* ..
26* .. Array Arguments ..
27* DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
28* $ X( LDX, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> DSBGST reduces a real symmetric-definite banded generalized
38*> eigenproblem A*x = lambda*B*x to standard form C*y = lambda*y,
39*> such that C has the same bandwidth as A.
40*>
41*> B must have been previously factorized as S**T*S by DPBSTF, using a
42*> split Cholesky factorization. A is overwritten by C = X**T*A*X, where
43*> X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
44*> bandwidth of A.
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] VECT
51*> \verbatim
52*> VECT is CHARACTER*1
53*> = 'N': do not form the transformation matrix X;
54*> = 'V': form X.
55*> \endverbatim
56*>
57*> \param[in] UPLO
58*> \verbatim
59*> UPLO is CHARACTER*1
60*> = 'U': Upper triangle of A is stored;
61*> = 'L': Lower triangle of A is stored.
62*> \endverbatim
63*>
64*> \param[in] N
65*> \verbatim
66*> N is INTEGER
67*> The order of the matrices A and B. N >= 0.
68*> \endverbatim
69*>
70*> \param[in] KA
71*> \verbatim
72*> KA is INTEGER
73*> The number of superdiagonals of the matrix A if UPLO = 'U',
74*> or the number of subdiagonals if UPLO = 'L'. KA >= 0.
75*> \endverbatim
76*>
77*> \param[in] KB
78*> \verbatim
79*> KB is INTEGER
80*> The number of superdiagonals of the matrix B if UPLO = 'U',
81*> or the number of subdiagonals if UPLO = 'L'. KA >= KB >= 0.
82*> \endverbatim
83*>
84*> \param[in,out] AB
85*> \verbatim
86*> AB is DOUBLE PRECISION array, dimension (LDAB,N)
87*> On entry, the upper or lower triangle of the symmetric band
88*> matrix A, stored in the first ka+1 rows of the array. The
89*> j-th column of A is stored in the j-th column of the array AB
90*> as follows:
91*> if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
92*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka).
93*>
94*> On exit, the transformed matrix X**T*A*X, stored in the same
95*> format as A.
96*> \endverbatim
97*>
98*> \param[in] LDAB
99*> \verbatim
100*> LDAB is INTEGER
101*> The leading dimension of the array AB. LDAB >= KA+1.
102*> \endverbatim
103*>
104*> \param[in] BB
105*> \verbatim
106*> BB is DOUBLE PRECISION array, dimension (LDBB,N)
107*> The banded factor S from the split Cholesky factorization of
108*> B, as returned by DPBSTF, stored in the first KB+1 rows of
109*> the array.
110*> \endverbatim
111*>
112*> \param[in] LDBB
113*> \verbatim
114*> LDBB is INTEGER
115*> The leading dimension of the array BB. LDBB >= KB+1.
116*> \endverbatim
117*>
118*> \param[out] X
119*> \verbatim
120*> X is DOUBLE PRECISION array, dimension (LDX,N)
121*> If VECT = 'V', the n-by-n matrix X.
122*> If VECT = 'N', the array X is not referenced.
123*> \endverbatim
124*>
125*> \param[in] LDX
126*> \verbatim
127*> LDX is INTEGER
128*> The leading dimension of the array X.
129*> LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
130*> \endverbatim
131*>
132*> \param[out] WORK
133*> \verbatim
134*> WORK is DOUBLE PRECISION array, dimension (2*N)
135*> \endverbatim
136*>
137*> \param[out] INFO
138*> \verbatim
139*> INFO is INTEGER
140*> = 0: successful exit
141*> < 0: if INFO = -i, the i-th argument had an illegal value.
142*> \endverbatim
143*
144* Authors:
145* ========
146*
147*> \author Univ. of Tennessee
148*> \author Univ. of California Berkeley
149*> \author Univ. of Colorado Denver
150*> \author NAG Ltd.
151*
152*> \ingroup hbgst
153*
154* =====================================================================
155 SUBROUTINE dsbgst( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB,
156 $ X,
157 $ LDX, WORK, INFO )
158*
159* -- LAPACK computational routine --
160* -- LAPACK is a software package provided by Univ. of Tennessee, --
161* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162*
163* .. Scalar Arguments ..
164 CHARACTER UPLO, VECT
165 INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N
166* ..
167* .. Array Arguments ..
168 DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
169 $ X( LDX, * )
170* ..
171*
172* =====================================================================
173*
174* .. Parameters ..
175 DOUBLE PRECISION ZERO, ONE
176 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
177* ..
178* .. Local Scalars ..
179 LOGICAL UPDATE, UPPER, WANTX
180 INTEGER I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
181 $ ka1, kb1, kbt, l, m, nr, nrt, nx
182 DOUBLE PRECISION BII, RA, RA1, T
183* ..
184* .. External Functions ..
185 LOGICAL LSAME
186 EXTERNAL LSAME
187* ..
188* .. External Subroutines ..
189 EXTERNAL dger, dlar2v, dlargv, dlartg, dlartv,
190 $ dlaset,
191 $ drot, dscal, 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( 'DSBGST', -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 dlaset( '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 DPBSTF. 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 dscal( n-m, one / bii, x( m+1, i ), 1 )
369 IF( kbt.GT.0 )
370 $ CALL dger( 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 dlartg( 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 dlargv( nrt, ab( 1, j2t ), inca, work( j2t-m ),
431 $ ka1,
432 $ work( n+j2t-m ), ka1 )
433 IF( nr.GT.0 ) THEN
434*
435* apply rotations in 1st set from the right
436*
437 DO 100 l = 1, ka - 1
438 CALL dlartv( nr, ab( ka1-l, j2 ), inca,
439 $ ab( ka-l, j2+1 ), inca, work( n+j2-m ),
440 $ work( j2-m ), ka1 )
441 100 CONTINUE
442*
443* apply rotations in 1st set from both sides to diagonal
444* blocks
445*
446 CALL dlar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
447 $ ab( ka, j2+1 ), inca, work( n+j2-m ),
448 $ work( j2-m ), ka1 )
449*
450 END IF
451*
452* start applying rotations in 1st set from the left
453*
454 DO 110 l = ka - 1, kb - k + 1, -1
455 nrt = ( n-j2+l ) / ka1
456 IF( nrt.GT.0 )
457 $ CALL dlartv( nrt, ab( l, j2+ka1-l ), inca,
458 $ ab( l+1, j2+ka1-l ), inca,
459 $ work( n+j2-m ), work( j2-m ), ka1 )
460 110 CONTINUE
461*
462 IF( wantx ) THEN
463*
464* post-multiply X by product of rotations in 1st set
465*
466 DO 120 j = j2, j1, ka1
467 CALL drot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
468 $ work( n+j-m ), work( j-m ) )
469 120 CONTINUE
470 END IF
471 130 CONTINUE
472*
473 IF( update ) THEN
474 IF( i2.LE.n .AND. kbt.GT.0 ) THEN
475*
476* create nonzero element a(i-kbt,i-kbt+ka+1) outside the
477* band and store it in WORK(i-kbt)
478*
479 work( i-kbt ) = -bb( kb1-kbt, i )*ra1
480 END IF
481 END IF
482*
483 DO 170 k = kb, 1, -1
484 IF( update ) THEN
485 j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
486 ELSE
487 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
488 END IF
489*
490* finish applying rotations in 2nd set from the left
491*
492 DO 140 l = kb - k, 1, -1
493 nrt = ( n-j2+ka+l ) / ka1
494 IF( nrt.GT.0 )
495 $ CALL dlartv( nrt, ab( l, j2-l+1 ), inca,
496 $ ab( l+1, j2-l+1 ), inca, work( n+j2-ka ),
497 $ work( j2-ka ), ka1 )
498 140 CONTINUE
499 nr = ( n-j2+ka ) / ka1
500 j1 = j2 + ( nr-1 )*ka1
501 DO 150 j = j1, j2, -ka1
502 work( j ) = work( j-ka )
503 work( n+j ) = work( n+j-ka )
504 150 CONTINUE
505 DO 160 j = j2, j1, ka1
506*
507* create nonzero element a(j-ka,j+1) outside the band
508* and store it in WORK(j)
509*
510 work( j ) = work( j )*ab( 1, j+1 )
511 ab( 1, j+1 ) = work( n+j )*ab( 1, j+1 )
512 160 CONTINUE
513 IF( update ) THEN
514 IF( i-k.LT.n-ka .AND. k.LE.kbt )
515 $ work( i-k+ka ) = work( i-k )
516 END IF
517 170 CONTINUE
518*
519 DO 210 k = kb, 1, -1
520 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
521 nr = ( n-j2+ka ) / ka1
522 j1 = j2 + ( nr-1 )*ka1
523 IF( nr.GT.0 ) THEN
524*
525* generate rotations in 2nd set to annihilate elements
526* which have been created outside the band
527*
528 CALL dlargv( nr, ab( 1, j2 ), inca, work( j2 ), ka1,
529 $ work( n+j2 ), ka1 )
530*
531* apply rotations in 2nd set from the right
532*
533 DO 180 l = 1, ka - 1
534 CALL dlartv( nr, ab( ka1-l, j2 ), inca,
535 $ ab( ka-l, j2+1 ), inca, work( n+j2 ),
536 $ work( j2 ), ka1 )
537 180 CONTINUE
538*
539* apply rotations in 2nd set from both sides to diagonal
540* blocks
541*
542 CALL dlar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
543 $ ab( ka, j2+1 ), inca, work( n+j2 ),
544 $ work( j2 ), ka1 )
545*
546 END IF
547*
548* start applying rotations in 2nd set from the left
549*
550 DO 190 l = ka - 1, kb - k + 1, -1
551 nrt = ( n-j2+l ) / ka1
552 IF( nrt.GT.0 )
553 $ CALL dlartv( nrt, ab( l, j2+ka1-l ), inca,
554 $ ab( l+1, j2+ka1-l ), inca, work( n+j2 ),
555 $ work( j2 ), ka1 )
556 190 CONTINUE
557*
558 IF( wantx ) THEN
559*
560* post-multiply X by product of rotations in 2nd set
561*
562 DO 200 j = j2, j1, ka1
563 CALL drot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
564 $ work( n+j ), work( j ) )
565 200 CONTINUE
566 END IF
567 210 CONTINUE
568*
569 DO 230 k = 1, kb - 1
570 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
571*
572* finish applying rotations in 1st set from the left
573*
574 DO 220 l = kb - k, 1, -1
575 nrt = ( n-j2+l ) / ka1
576 IF( nrt.GT.0 )
577 $ CALL dlartv( nrt, ab( l, j2+ka1-l ), inca,
578 $ ab( l+1, j2+ka1-l ), inca,
579 $ work( n+j2-m ), work( j2-m ), ka1 )
580 220 CONTINUE
581 230 CONTINUE
582*
583 IF( kb.GT.1 ) THEN
584 DO 240 j = n - 1, i - kb + 2*ka + 1, -1
585 work( n+j-m ) = work( n+j-ka-m )
586 work( j-m ) = work( j-ka-m )
587 240 CONTINUE
588 END IF
589*
590 ELSE
591*
592* Transform A, working with the lower triangle
593*
594 IF( update ) THEN
595*
596* Form inv(S(i))**T * A * inv(S(i))
597*
598 bii = bb( 1, i )
599 DO 250 j = i, i1
600 ab( j-i+1, i ) = ab( j-i+1, i ) / bii
601 250 CONTINUE
602 DO 260 j = max( 1, i-ka ), i
603 ab( i-j+1, j ) = ab( i-j+1, j ) / bii
604 260 CONTINUE
605 DO 290 k = i - kbt, i - 1
606 DO 270 j = i - kbt, k
607 ab( k-j+1, j ) = ab( k-j+1, j ) -
608 $ bb( i-j+1, j )*ab( i-k+1, k ) -
609 $ bb( i-k+1, k )*ab( i-j+1, j ) +
610 $ ab( 1, i )*bb( i-j+1, j )*
611 $ bb( i-k+1, k )
612 270 CONTINUE
613 DO 280 j = max( 1, i-ka ), i - kbt - 1
614 ab( k-j+1, j ) = ab( k-j+1, j ) -
615 $ bb( i-k+1, k )*ab( i-j+1, j )
616 280 CONTINUE
617 290 CONTINUE
618 DO 310 j = i, i1
619 DO 300 k = max( j-ka, i-kbt ), i - 1
620 ab( j-k+1, k ) = ab( j-k+1, k ) -
621 $ bb( i-k+1, k )*ab( j-i+1, i )
622 300 CONTINUE
623 310 CONTINUE
624*
625 IF( wantx ) THEN
626*
627* post-multiply X by inv(S(i))
628*
629 CALL dscal( n-m, one / bii, x( m+1, i ), 1 )
630 IF( kbt.GT.0 )
631 $ CALL dger( n-m, kbt, -one, x( m+1, i ), 1,
632 $ bb( kbt+1, i-kbt ), ldbb-1,
633 $ x( m+1, i-kbt ), ldx )
634 END IF
635*
636* store a(i1,i) in RA1 for use in next loop over K
637*
638 ra1 = ab( i1-i+1, i )
639 END IF
640*
641* Generate and apply vectors of rotations to chase all the
642* existing bulges KA positions down toward the bottom of the
643* band
644*
645 DO 360 k = 1, kb - 1
646 IF( update ) THEN
647*
648* Determine the rotations which would annihilate the bulge
649* which has in theory just been created
650*
651 IF( i-k+ka.LT.n .AND. i-k.GT.1 ) THEN
652*
653* generate rotation to annihilate a(i-k+ka+1,i)
654*
655 CALL dlartg( ab( ka1-k, i ), ra1,
656 $ work( n+i-k+ka-m ),
657 $ work( i-k+ka-m ), ra )
658*
659* create nonzero element a(i-k+ka+1,i-k) outside the
660* band and store it in WORK(i-k)
661*
662 t = -bb( k+1, i-k )*ra1
663 work( i-k ) = work( n+i-k+ka-m )*t -
664 $ work( i-k+ka-m )*ab( ka1, i-k )
665 ab( ka1, i-k ) = work( i-k+ka-m )*t +
666 $ work( n+i-k+ka-m )*ab( ka1, i-k )
667 ra1 = ra
668 END IF
669 END IF
670 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
671 nr = ( n-j2+ka ) / ka1
672 j1 = j2 + ( nr-1 )*ka1
673 IF( update ) THEN
674 j2t = max( j2, i+2*ka-k+1 )
675 ELSE
676 j2t = j2
677 END IF
678 nrt = ( n-j2t+ka ) / ka1
679 DO 320 j = j2t, j1, ka1
680*
681* create nonzero element a(j+1,j-ka) outside the band
682* and store it in WORK(j-m)
683*
684 work( j-m ) = work( j-m )*ab( ka1, j-ka+1 )
685 ab( ka1, j-ka+1 ) = work( n+j-m )*ab( ka1, j-ka+1 )
686 320 CONTINUE
687*
688* generate rotations in 1st set to annihilate elements which
689* have been created outside the band
690*
691 IF( nrt.GT.0 )
692 $ CALL dlargv( nrt, ab( ka1, j2t-ka ), inca,
693 $ 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,
709 $ j2 ),
710 $ inca, work( n+j2-m ), work( j2-m ), ka1 )
711*
712 END IF
713*
714* start applying rotations in 1st set from the right
715*
716 DO 340 l = ka - 1, kb - k + 1, -1
717 nrt = ( n-j2+l ) / ka1
718 IF( nrt.GT.0 )
719 $ CALL dlartv( nrt, ab( ka1-l+1, j2 ), inca,
720 $ ab( ka1-l, j2+1 ), inca, work( n+j2-m ),
721 $ work( j2-m ), ka1 )
722 340 CONTINUE
723*
724 IF( wantx ) THEN
725*
726* post-multiply X by product of rotations in 1st set
727*
728 DO 350 j = j2, j1, ka1
729 CALL drot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
730 $ work( n+j-m ), work( j-m ) )
731 350 CONTINUE
732 END IF
733 360 CONTINUE
734*
735 IF( update ) THEN
736 IF( i2.LE.n .AND. kbt.GT.0 ) THEN
737*
738* create nonzero element a(i-kbt+ka+1,i-kbt) outside the
739* band and store it in WORK(i-kbt)
740*
741 work( i-kbt ) = -bb( kbt+1, i-kbt )*ra1
742 END IF
743 END IF
744*
745 DO 400 k = kb, 1, -1
746 IF( update ) THEN
747 j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
748 ELSE
749 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
750 END IF
751*
752* finish applying rotations in 2nd set from the right
753*
754 DO 370 l = kb - k, 1, -1
755 nrt = ( n-j2+ka+l ) / ka1
756 IF( nrt.GT.0 )
757 $ CALL dlartv( nrt, ab( ka1-l+1, j2-ka ), inca,
758 $ ab( ka1-l, j2-ka+1 ), inca,
759 $ work( n+j2-ka ), work( j2-ka ), ka1 )
760 370 CONTINUE
761 nr = ( n-j2+ka ) / ka1
762 j1 = j2 + ( nr-1 )*ka1
763 DO 380 j = j1, j2, -ka1
764 work( j ) = work( j-ka )
765 work( n+j ) = work( n+j-ka )
766 380 CONTINUE
767 DO 390 j = j2, j1, ka1
768*
769* create nonzero element a(j+1,j-ka) outside the band
770* and store it in WORK(j)
771*
772 work( j ) = work( j )*ab( ka1, j-ka+1 )
773 ab( ka1, j-ka+1 ) = work( n+j )*ab( ka1, j-ka+1 )
774 390 CONTINUE
775 IF( update ) THEN
776 IF( i-k.LT.n-ka .AND. k.LE.kbt )
777 $ work( i-k+ka ) = work( i-k )
778 END IF
779 400 CONTINUE
780*
781 DO 440 k = kb, 1, -1
782 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
783 nr = ( n-j2+ka ) / ka1
784 j1 = j2 + ( nr-1 )*ka1
785 IF( nr.GT.0 ) THEN
786*
787* generate rotations in 2nd set to annihilate elements
788* which have been created outside the band
789*
790 CALL dlargv( nr, ab( ka1, j2-ka ), inca, work( j2 ),
791 $ ka1,
792 $ work( n+j2 ), ka1 )
793*
794* apply rotations in 2nd set from the left
795*
796 DO 410 l = 1, ka - 1
797 CALL dlartv( nr, ab( l+1, j2-l ), inca,
798 $ ab( l+2, j2-l ), inca, work( n+j2 ),
799 $ work( j2 ), ka1 )
800 410 CONTINUE
801*
802* apply rotations in 2nd set from both sides to diagonal
803* blocks
804*
805 CALL dlar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2,
806 $ j2 ),
807 $ inca, work( n+j2 ), work( j2 ), ka1 )
808*
809 END IF
810*
811* start applying rotations in 2nd set from the right
812*
813 DO 420 l = ka - 1, kb - k + 1, -1
814 nrt = ( n-j2+l ) / ka1
815 IF( nrt.GT.0 )
816 $ CALL dlartv( nrt, ab( ka1-l+1, j2 ), inca,
817 $ ab( ka1-l, j2+1 ), inca, work( n+j2 ),
818 $ work( j2 ), ka1 )
819 420 CONTINUE
820*
821 IF( wantx ) THEN
822*
823* post-multiply X by product of rotations in 2nd set
824*
825 DO 430 j = j2, j1, ka1
826 CALL drot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
827 $ work( n+j ), work( j ) )
828 430 CONTINUE
829 END IF
830 440 CONTINUE
831*
832 DO 460 k = 1, kb - 1
833 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
834*
835* finish applying rotations in 1st set from the right
836*
837 DO 450 l = kb - k, 1, -1
838 nrt = ( n-j2+l ) / ka1
839 IF( nrt.GT.0 )
840 $ CALL dlartv( nrt, ab( ka1-l+1, j2 ), inca,
841 $ ab( ka1-l, j2+1 ), inca, work( n+j2-m ),
842 $ work( j2-m ), ka1 )
843 450 CONTINUE
844 460 CONTINUE
845*
846 IF( kb.GT.1 ) THEN
847 DO 470 j = n - 1, i - kb + 2*ka + 1, -1
848 work( n+j-m ) = work( n+j-ka-m )
849 work( j-m ) = work( j-ka-m )
850 470 CONTINUE
851 END IF
852*
853 END IF
854*
855 GO TO 10
856*
857 480 CONTINUE
858*
859* **************************** Phase 2 *****************************
860*
861* The logical structure of this phase is:
862*
863* UPDATE = .TRUE.
864* DO I = 1, M
865* use S(i) to update A and create a new bulge
866* apply rotations to push all bulges KA positions upward
867* END DO
868* UPDATE = .FALSE.
869* DO I = M - KA - 1, 2, -1
870* apply rotations to push all bulges KA positions upward
871* END DO
872*
873* To avoid duplicating code, the two loops are merged.
874*
875 update = .true.
876 i = 0
877 490 CONTINUE
878 IF( update ) THEN
879 i = i + 1
880 kbt = min( kb, m-i )
881 i0 = i + 1
882 i1 = max( 1, i-ka )
883 i2 = i + kbt - ka1
884 IF( i.GT.m ) THEN
885 update = .false.
886 i = i - 1
887 i0 = m + 1
888 IF( ka.EQ.0 )
889 $ RETURN
890 GO TO 490
891 END IF
892 ELSE
893 i = i - ka
894 IF( i.LT.2 )
895 $ RETURN
896 END IF
897*
898 IF( i.LT.m-kbt ) THEN
899 nx = m
900 ELSE
901 nx = n
902 END IF
903*
904 IF( upper ) THEN
905*
906* Transform A, working with the upper triangle
907*
908 IF( update ) THEN
909*
910* Form inv(S(i))**T * A * inv(S(i))
911*
912 bii = bb( kb1, i )
913 DO 500 j = i1, i
914 ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
915 500 CONTINUE
916 DO 510 j = i, min( n, i+ka )
917 ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
918 510 CONTINUE
919 DO 540 k = i + 1, i + kbt
920 DO 520 j = k, i + kbt
921 ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
922 $ bb( i-j+kb1, j )*ab( i-k+ka1, k ) -
923 $ bb( i-k+kb1, k )*ab( i-j+ka1, j ) +
924 $ ab( ka1, i )*bb( i-j+kb1, j )*
925 $ bb( i-k+kb1, k )
926 520 CONTINUE
927 DO 530 j = i + kbt + 1, min( n, i+ka )
928 ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
929 $ bb( i-k+kb1, k )*ab( i-j+ka1, j )
930 530 CONTINUE
931 540 CONTINUE
932 DO 560 j = i1, i
933 DO 550 k = i + 1, min( j+ka, i+kbt )
934 ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
935 $ bb( i-k+kb1, k )*ab( j-i+ka1, i )
936 550 CONTINUE
937 560 CONTINUE
938*
939 IF( wantx ) THEN
940*
941* post-multiply X by inv(S(i))
942*
943 CALL dscal( nx, one / bii, x( 1, i ), 1 )
944 IF( kbt.GT.0 )
945 $ CALL dger( nx, kbt, -one, x( 1, i ), 1, bb( kb,
946 $ i+1 ),
947 $ ldbb-1, x( 1, i+1 ), ldx )
948 END IF
949*
950* store a(i1,i) in RA1 for use in next loop over K
951*
952 ra1 = ab( i1-i+ka1, i )
953 END IF
954*
955* Generate and apply vectors of rotations to chase all the
956* existing bulges KA positions up toward the top of the band
957*
958 DO 610 k = 1, kb - 1
959 IF( update ) THEN
960*
961* Determine the rotations which would annihilate the bulge
962* which has in theory just been created
963*
964 IF( i+k-ka1.GT.0 .AND. i+k.LT.m ) THEN
965*
966* generate rotation to annihilate a(i+k-ka-1,i)
967*
968 CALL dlartg( ab( k+1, i ), ra1, work( n+i+k-ka ),
969 $ work( i+k-ka ), ra )
970*
971* create nonzero element a(i+k-ka-1,i+k) outside the
972* band and store it in WORK(m-kb+i+k)
973*
974 t = -bb( kb1-k, i+k )*ra1
975 work( m-kb+i+k ) = work( n+i+k-ka )*t -
976 $ work( i+k-ka )*ab( 1, i+k )
977 ab( 1, i+k ) = work( i+k-ka )*t +
978 $ work( n+i+k-ka )*ab( 1, i+k )
979 ra1 = ra
980 END IF
981 END IF
982 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
983 nr = ( j2+ka-1 ) / ka1
984 j1 = j2 - ( nr-1 )*ka1
985 IF( update ) THEN
986 j2t = min( j2, i-2*ka+k-1 )
987 ELSE
988 j2t = j2
989 END IF
990 nrt = ( j2t+ka-1 ) / ka1
991 DO 570 j = j1, j2t, ka1
992*
993* create nonzero element a(j-1,j+ka) outside the band
994* and store it in WORK(j)
995*
996 work( j ) = work( j )*ab( 1, j+ka-1 )
997 ab( 1, j+ka-1 ) = work( n+j )*ab( 1, j+ka-1 )
998 570 CONTINUE
999*
1000* generate rotations in 1st set to annihilate elements which
1001* have been created outside the band
1002*
1003 IF( nrt.GT.0 )
1004 $ CALL dlargv( nrt, ab( 1, j1+ka ), inca, work( j1 ),
1005 $ ka1,
1006 $ work( n+j1 ), ka1 )
1007 IF( nr.GT.0 ) THEN
1008*
1009* apply rotations in 1st set from the left
1010*
1011 DO 580 l = 1, ka - 1
1012 CALL dlartv( nr, ab( ka1-l, j1+l ), inca,
1013 $ ab( ka-l, j1+l ), inca, work( n+j1 ),
1014 $ work( j1 ), ka1 )
1015 580 CONTINUE
1016*
1017* apply rotations in 1st set from both sides to diagonal
1018* blocks
1019*
1020 CALL dlar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1021 $ ab( ka, j1 ), inca, work( n+j1 ),
1022 $ work( j1 ), ka1 )
1023*
1024 END IF
1025*
1026* start applying rotations in 1st set from the right
1027*
1028 DO 590 l = ka - 1, kb - k + 1, -1
1029 nrt = ( j2+l-1 ) / ka1
1030 j1t = j2 - ( nrt-1 )*ka1
1031 IF( nrt.GT.0 )
1032 $ CALL dlartv( nrt, ab( l, j1t ), inca,
1033 $ ab( l+1, j1t-1 ), inca, work( n+j1t ),
1034 $ work( j1t ), ka1 )
1035 590 CONTINUE
1036*
1037 IF( wantx ) THEN
1038*
1039* post-multiply X by product of rotations in 1st set
1040*
1041 DO 600 j = j1, j2, ka1
1042 CALL drot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1043 $ work( n+j ), work( j ) )
1044 600 CONTINUE
1045 END IF
1046 610 CONTINUE
1047*
1048 IF( update ) THEN
1049 IF( i2.GT.0 .AND. kbt.GT.0 ) THEN
1050*
1051* create nonzero element a(i+kbt-ka-1,i+kbt) outside the
1052* band and store it in WORK(m-kb+i+kbt)
1053*
1054 work( m-kb+i+kbt ) = -bb( kb1-kbt, i+kbt )*ra1
1055 END IF
1056 END IF
1057*
1058 DO 650 k = kb, 1, -1
1059 IF( update ) THEN
1060 j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1061 ELSE
1062 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1063 END IF
1064*
1065* finish applying rotations in 2nd set from the right
1066*
1067 DO 620 l = kb - k, 1, -1
1068 nrt = ( j2+ka+l-1 ) / ka1
1069 j1t = j2 - ( nrt-1 )*ka1
1070 IF( nrt.GT.0 )
1071 $ CALL dlartv( nrt, ab( l, j1t+ka ), inca,
1072 $ ab( l+1, j1t+ka-1 ), inca,
1073 $ work( n+m-kb+j1t+ka ),
1074 $ work( m-kb+j1t+ka ), ka1 )
1075 620 CONTINUE
1076 nr = ( j2+ka-1 ) / ka1
1077 j1 = j2 - ( nr-1 )*ka1
1078 DO 630 j = j1, j2, ka1
1079 work( m-kb+j ) = work( m-kb+j+ka )
1080 work( n+m-kb+j ) = work( n+m-kb+j+ka )
1081 630 CONTINUE
1082 DO 640 j = j1, j2, ka1
1083*
1084* create nonzero element a(j-1,j+ka) outside the band
1085* and store it in WORK(m-kb+j)
1086*
1087 work( m-kb+j ) = work( m-kb+j )*ab( 1, j+ka-1 )
1088 ab( 1, j+ka-1 ) = work( n+m-kb+j )*ab( 1, j+ka-1 )
1089 640 CONTINUE
1090 IF( update ) THEN
1091 IF( i+k.GT.ka1 .AND. k.LE.kbt )
1092 $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1093 END IF
1094 650 CONTINUE
1095*
1096 DO 690 k = kb, 1, -1
1097 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1098 nr = ( j2+ka-1 ) / ka1
1099 j1 = j2 - ( nr-1 )*ka1
1100 IF( nr.GT.0 ) THEN
1101*
1102* generate rotations in 2nd set to annihilate elements
1103* which have been created outside the band
1104*
1105 CALL dlargv( nr, ab( 1, j1+ka ), inca,
1106 $ work( m-kb+j1 ),
1107 $ ka1, work( n+m-kb+j1 ), ka1 )
1108*
1109* apply rotations in 2nd set from the left
1110*
1111 DO 660 l = 1, ka - 1
1112 CALL dlartv( nr, ab( ka1-l, j1+l ), inca,
1113 $ ab( ka-l, j1+l ), inca,
1114 $ work( n+m-kb+j1 ), work( m-kb+j1 ), ka1 )
1115 660 CONTINUE
1116*
1117* apply rotations in 2nd set from both sides to diagonal
1118* blocks
1119*
1120 CALL dlar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1121 $ ab( ka, j1 ), inca, work( n+m-kb+j1 ),
1122 $ work( m-kb+j1 ), ka1 )
1123*
1124 END IF
1125*
1126* start applying rotations in 2nd set from the right
1127*
1128 DO 670 l = ka - 1, kb - k + 1, -1
1129 nrt = ( j2+l-1 ) / ka1
1130 j1t = j2 - ( nrt-1 )*ka1
1131 IF( nrt.GT.0 )
1132 $ CALL dlartv( nrt, ab( l, j1t ), inca,
1133 $ ab( l+1, j1t-1 ), inca,
1134 $ work( n+m-kb+j1t ), work( m-kb+j1t ),
1135 $ ka1 )
1136 670 CONTINUE
1137*
1138 IF( wantx ) THEN
1139*
1140* post-multiply X by product of rotations in 2nd set
1141*
1142 DO 680 j = j1, j2, ka1
1143 CALL drot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1144 $ work( n+m-kb+j ), work( m-kb+j ) )
1145 680 CONTINUE
1146 END IF
1147 690 CONTINUE
1148*
1149 DO 710 k = 1, kb - 1
1150 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1151*
1152* finish applying rotations in 1st set from the right
1153*
1154 DO 700 l = kb - k, 1, -1
1155 nrt = ( j2+l-1 ) / ka1
1156 j1t = j2 - ( nrt-1 )*ka1
1157 IF( nrt.GT.0 )
1158 $ CALL dlartv( nrt, ab( l, j1t ), inca,
1159 $ ab( l+1, j1t-1 ), inca, work( n+j1t ),
1160 $ work( j1t ), ka1 )
1161 700 CONTINUE
1162 710 CONTINUE
1163*
1164 IF( kb.GT.1 ) THEN
1165 DO 720 j = 2, min( i+kb, m ) - 2*ka - 1
1166 work( n+j ) = work( n+j+ka )
1167 work( j ) = work( j+ka )
1168 720 CONTINUE
1169 END IF
1170*
1171 ELSE
1172*
1173* Transform A, working with the lower triangle
1174*
1175 IF( update ) THEN
1176*
1177* Form inv(S(i))**T * A * inv(S(i))
1178*
1179 bii = bb( 1, i )
1180 DO 730 j = i1, i
1181 ab( i-j+1, j ) = ab( i-j+1, j ) / bii
1182 730 CONTINUE
1183 DO 740 j = i, min( n, i+ka )
1184 ab( j-i+1, i ) = ab( j-i+1, i ) / bii
1185 740 CONTINUE
1186 DO 770 k = i + 1, i + kbt
1187 DO 750 j = k, i + kbt
1188 ab( j-k+1, k ) = ab( j-k+1, k ) -
1189 $ bb( j-i+1, i )*ab( k-i+1, i ) -
1190 $ bb( k-i+1, i )*ab( j-i+1, i ) +
1191 $ ab( 1, i )*bb( j-i+1, i )*
1192 $ bb( k-i+1, i )
1193 750 CONTINUE
1194 DO 760 j = i + kbt + 1, min( n, i+ka )
1195 ab( j-k+1, k ) = ab( j-k+1, k ) -
1196 $ bb( k-i+1, i )*ab( j-i+1, i )
1197 760 CONTINUE
1198 770 CONTINUE
1199 DO 790 j = i1, i
1200 DO 780 k = i + 1, min( j+ka, i+kbt )
1201 ab( k-j+1, j ) = ab( k-j+1, j ) -
1202 $ bb( k-i+1, i )*ab( i-j+1, j )
1203 780 CONTINUE
1204 790 CONTINUE
1205*
1206 IF( wantx ) THEN
1207*
1208* post-multiply X by inv(S(i))
1209*
1210 CALL dscal( nx, one / bii, x( 1, i ), 1 )
1211 IF( kbt.GT.0 )
1212 $ CALL dger( nx, kbt, -one, x( 1, i ), 1, bb( 2, i ),
1213 $ 1,
1214 $ x( 1, i+1 ), ldx )
1215 END IF
1216*
1217* store a(i,i1) in RA1 for use in next loop over K
1218*
1219 ra1 = ab( i-i1+1, i1 )
1220 END IF
1221*
1222* Generate and apply vectors of rotations to chase all the
1223* existing bulges KA positions up toward the top of the band
1224*
1225 DO 840 k = 1, kb - 1
1226 IF( update ) THEN
1227*
1228* Determine the rotations which would annihilate the bulge
1229* which has in theory just been created
1230*
1231 IF( i+k-ka1.GT.0 .AND. i+k.LT.m ) THEN
1232*
1233* generate rotation to annihilate a(i,i+k-ka-1)
1234*
1235 CALL dlartg( ab( ka1-k, i+k-ka ), ra1,
1236 $ work( n+i+k-ka ), work( i+k-ka ), ra )
1237*
1238* create nonzero element a(i+k,i+k-ka-1) outside the
1239* band and store it in WORK(m-kb+i+k)
1240*
1241 t = -bb( k+1, i )*ra1
1242 work( m-kb+i+k ) = work( n+i+k-ka )*t -
1243 $ work( i+k-ka )*ab( ka1, i+k-ka )
1244 ab( ka1, i+k-ka ) = work( i+k-ka )*t +
1245 $ work( n+i+k-ka )*ab( ka1, i+k-ka )
1246 ra1 = ra
1247 END IF
1248 END IF
1249 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1250 nr = ( j2+ka-1 ) / ka1
1251 j1 = j2 - ( nr-1 )*ka1
1252 IF( update ) THEN
1253 j2t = min( j2, i-2*ka+k-1 )
1254 ELSE
1255 j2t = j2
1256 END IF
1257 nrt = ( j2t+ka-1 ) / ka1
1258 DO 800 j = j1, j2t, ka1
1259*
1260* create nonzero element a(j+ka,j-1) outside the band
1261* and store it in WORK(j)
1262*
1263 work( j ) = work( j )*ab( ka1, j-1 )
1264 ab( ka1, j-1 ) = work( n+j )*ab( ka1, j-1 )
1265 800 CONTINUE
1266*
1267* generate rotations in 1st set to annihilate elements which
1268* have been created outside the band
1269*
1270 IF( nrt.GT.0 )
1271 $ CALL dlargv( nrt, ab( ka1, j1 ), inca, work( j1 ),
1272 $ ka1,
1273 $ work( n+j1 ), ka1 )
1274 IF( nr.GT.0 ) THEN
1275*
1276* apply rotations in 1st set from the right
1277*
1278 DO 810 l = 1, ka - 1
1279 CALL dlartv( nr, ab( l+1, j1 ), inca, ab( l+2,
1280 $ j1-1 ),
1281 $ inca, work( n+j1 ), work( j1 ), ka1 )
1282 810 CONTINUE
1283*
1284* apply rotations in 1st set from both sides to diagonal
1285* blocks
1286*
1287 CALL dlar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1288 $ ab( 2, j1-1 ), inca, work( n+j1 ),
1289 $ work( j1 ), ka1 )
1290*
1291 END IF
1292*
1293* start applying rotations in 1st set from the left
1294*
1295 DO 820 l = ka - 1, kb - k + 1, -1
1296 nrt = ( j2+l-1 ) / ka1
1297 j1t = j2 - ( nrt-1 )*ka1
1298 IF( nrt.GT.0 )
1299 $ CALL dlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1300 $ ab( ka1-l, j1t-ka1+l ), inca,
1301 $ work( n+j1t ), work( j1t ), ka1 )
1302 820 CONTINUE
1303*
1304 IF( wantx ) THEN
1305*
1306* post-multiply X by product of rotations in 1st set
1307*
1308 DO 830 j = j1, j2, ka1
1309 CALL drot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1310 $ work( n+j ), work( j ) )
1311 830 CONTINUE
1312 END IF
1313 840 CONTINUE
1314*
1315 IF( update ) THEN
1316 IF( i2.GT.0 .AND. kbt.GT.0 ) THEN
1317*
1318* create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1319* band and store it in WORK(m-kb+i+kbt)
1320*
1321 work( m-kb+i+kbt ) = -bb( kbt+1, i )*ra1
1322 END IF
1323 END IF
1324*
1325 DO 880 k = kb, 1, -1
1326 IF( update ) THEN
1327 j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1328 ELSE
1329 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1330 END IF
1331*
1332* finish applying rotations in 2nd set from the left
1333*
1334 DO 850 l = kb - k, 1, -1
1335 nrt = ( j2+ka+l-1 ) / ka1
1336 j1t = j2 - ( nrt-1 )*ka1
1337 IF( nrt.GT.0 )
1338 $ CALL dlartv( nrt, ab( ka1-l+1, j1t+l-1 ), inca,
1339 $ ab( ka1-l, j1t+l-1 ), inca,
1340 $ work( n+m-kb+j1t+ka ),
1341 $ work( m-kb+j1t+ka ), ka1 )
1342 850 CONTINUE
1343 nr = ( j2+ka-1 ) / ka1
1344 j1 = j2 - ( nr-1 )*ka1
1345 DO 860 j = j1, j2, ka1
1346 work( m-kb+j ) = work( m-kb+j+ka )
1347 work( n+m-kb+j ) = work( n+m-kb+j+ka )
1348 860 CONTINUE
1349 DO 870 j = j1, j2, ka1
1350*
1351* create nonzero element a(j+ka,j-1) outside the band
1352* and store it in WORK(m-kb+j)
1353*
1354 work( m-kb+j ) = work( m-kb+j )*ab( ka1, j-1 )
1355 ab( ka1, j-1 ) = work( n+m-kb+j )*ab( ka1, j-1 )
1356 870 CONTINUE
1357 IF( update ) THEN
1358 IF( i+k.GT.ka1 .AND. k.LE.kbt )
1359 $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1360 END IF
1361 880 CONTINUE
1362*
1363 DO 920 k = kb, 1, -1
1364 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1365 nr = ( j2+ka-1 ) / ka1
1366 j1 = j2 - ( nr-1 )*ka1
1367 IF( nr.GT.0 ) THEN
1368*
1369* generate rotations in 2nd set to annihilate elements
1370* which have been created outside the band
1371*
1372 CALL dlargv( nr, ab( ka1, j1 ), inca, work( m-kb+j1 ),
1373 $ ka1, work( n+m-kb+j1 ), ka1 )
1374*
1375* apply rotations in 2nd set from the right
1376*
1377 DO 890 l = 1, ka - 1
1378 CALL dlartv( nr, ab( l+1, j1 ), inca, ab( l+2,
1379 $ j1-1 ),
1380 $ inca, work( n+m-kb+j1 ), work( m-kb+j1 ),
1381 $ ka1 )
1382 890 CONTINUE
1383*
1384* apply rotations in 2nd set from both sides to diagonal
1385* blocks
1386*
1387 CALL dlar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1388 $ ab( 2, j1-1 ), inca, work( n+m-kb+j1 ),
1389 $ work( m-kb+j1 ), ka1 )
1390*
1391 END IF
1392*
1393* start applying rotations in 2nd set from the left
1394*
1395 DO 900 l = ka - 1, kb - k + 1, -1
1396 nrt = ( j2+l-1 ) / ka1
1397 j1t = j2 - ( nrt-1 )*ka1
1398 IF( nrt.GT.0 )
1399 $ CALL dlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1400 $ ab( ka1-l, j1t-ka1+l ), inca,
1401 $ work( n+m-kb+j1t ), work( m-kb+j1t ),
1402 $ ka1 )
1403 900 CONTINUE
1404*
1405 IF( wantx ) THEN
1406*
1407* post-multiply X by product of rotations in 2nd set
1408*
1409 DO 910 j = j1, j2, ka1
1410 CALL drot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1411 $ work( n+m-kb+j ), work( m-kb+j ) )
1412 910 CONTINUE
1413 END IF
1414 920 CONTINUE
1415*
1416 DO 940 k = 1, kb - 1
1417 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1418*
1419* finish applying rotations in 1st set from the left
1420*
1421 DO 930 l = kb - k, 1, -1
1422 nrt = ( j2+l-1 ) / ka1
1423 j1t = j2 - ( nrt-1 )*ka1
1424 IF( nrt.GT.0 )
1425 $ CALL dlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1426 $ ab( ka1-l, j1t-ka1+l ), inca,
1427 $ work( n+j1t ), work( j1t ), ka1 )
1428 930 CONTINUE
1429 940 CONTINUE
1430*
1431 IF( kb.GT.1 ) THEN
1432 DO 950 j = 2, min( i+kb, m ) - 2*ka - 1
1433 work( n+j ) = work( n+j+ka )
1434 work( j ) = work( j+ka )
1435 950 CONTINUE
1436 END IF
1437*
1438 END IF
1439*
1440 GO TO 490
1441*
1442* End of DSBGST
1443*
1444 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
dsbgst
subroutine dsbgst(vect, uplo, n, ka, kb, ab, ldab, bb, ldbb, x, ldx, work, info)
DSBGST
Definition dsbgst.f:158
dlar2v
subroutine dlar2v(n, x, y, z, incx, c, s, incc)
DLAR2V applies a vector of plane rotations with real cosines and real sines from both sides to a sequ...
Definition dlar2v.f:108
dlargv
subroutine dlargv(n, x, incx, y, incy, c, incc)
DLARGV generates a vector of plane rotations with real cosines and real sines.
Definition dlargv.f:102
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlartv
subroutine dlartv(n, x, incx, y, incy, c, s, incc)
DLARTV applies a vector of plane rotations with real cosines and real sines to the elements of a pair...
Definition dlartv.f:106
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:108
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79