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