OpenRadioss 2025.1.11
OpenRadioss project
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*> \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*> \ingroup doubleOTHERcomputational
155*
156* =====================================================================
157 SUBROUTINE dsbgst( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
158 $ LDX, WORK, INFO )
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 DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
170 $ x( ldx, * )
171* ..
172*
173* =====================================================================
174*
175* .. Parameters ..
176 DOUBLE PRECISION ZERO, ONE
177 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION BII, RA, RA1, T
184* ..
185* .. External Functions ..
186 LOGICAL LSAME
187 EXTERNAL lsame
188* ..
189* .. External Subroutines ..
190 EXTERNAL dger, dlar2v, dlargv, dlartg, dlartv, 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 ), 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 dlartv( 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 dlar2v( 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 dlartv( 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 drot( 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 dlartv( 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 dlargv( 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 dlartv( 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 dlar2v( 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 dlartv( 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 drot( 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 dlartv( 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 dscal( n-m, one / bii, x( m+1, i ), 1 )
629 IF( kbt.GT.0 )
630 $ CALL dger( 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 dlartg( 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 dlargv( 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 dlartv( 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 dlar2v( 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 dlartv( 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 drot( 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 dlartv( 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 dlargv( 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 dlartv( 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 dlar2v( 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 dlartv( 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 drot( 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 dlartv( 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 dscal( nx, one / bii, x( 1, i ), 1 )
938 IF( kbt.GT.0 )
939 $ CALL dger( 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 dlartg( 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 dlargv( 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 dlartv( 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 dlar2v( 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 dlartv( 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 drot( 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 dlartv( 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 dlargv( 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 dlartv( 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 dlar2v( 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 dlartv( 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 drot( 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 dlartv( 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 dscal( nx, one / bii, x( 1, i ), 1 )
1202 IF( kbt.GT.0 )
1203 $ CALL dger( 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 dlartg( 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 dlargv( 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 dlartv( 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 dlar2v( 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 dlartv( 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 drot( 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 dlartv( 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 dlargv( 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 dlartv( 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 dlar2v( 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 dlartv( 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 drot( 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 dlartv( 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 DSBGST
1430*
1431 END
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:113
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:110
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
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:110
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: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:104
dsbgst
subroutine dsbgst(vect, uplo, n, ka, kb, ab, ldab, bb, ldbb, x, ldx, work, info)
DSBGST
Definition dsbgst.f:159
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21