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