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