1: *> \brief \b ZHBGST
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZHBGST + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbgst.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbgst.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbgst.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHBGST( 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: * DOUBLE PRECISION RWORK( * )
30: * COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
31: * $ X( LDX, * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> ZHBGST 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 ZPBSTF, 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*16 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*16 array, dimension (LDBB,N)
110: *> The banded factor S from the split Cholesky factorization of
111: *> B, as returned by ZPBSTF, 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*16 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*16 array, dimension (N)
138: *> \endverbatim
139: *>
140: *> \param[out] RWORK
141: *> \verbatim
142: *> RWORK is DOUBLE PRECISION 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 complex16OTHERcomputational
161: *
162: * =====================================================================
163: SUBROUTINE ZHBGST( 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: DOUBLE PRECISION RWORK( * )
176: COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
177: $ X( LDX, * )
178: * ..
179: *
180: * =====================================================================
181: *
182: * .. Parameters ..
183: COMPLEX*16 CZERO, CONE
184: DOUBLE PRECISION ONE
185: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
186: $ CONE = ( 1.0D+0, 0.0D+0 ), ONE = 1.0D+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: DOUBLE PRECISION BII
193: COMPLEX*16 RA, RA1, T
194: * ..
195: * .. External Functions ..
196: LOGICAL LSAME
197: EXTERNAL LSAME
198: * ..
199: * .. External Subroutines ..
200: EXTERNAL XERBLA, ZDSCAL, ZGERC, ZGERU, ZLACGV, ZLAR2V,
201: $ ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT
202: * ..
203: * .. Intrinsic Functions ..
204: INTRINSIC DBLE, DCONJG, MAX, MIN
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( 'ZHBGST', -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 ZLASET( '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 ZPBSTF. 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 = DBLE( BB( KB1, I ) )
346: AB( KA1, I ) = ( DBLE( 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: $ DCONJG( AB( K-I+KA1, I ) ) -
358: $ DCONJG( BB( K-I+KB1, I ) )*
359: $ AB( J-I+KA1, I ) +
360: $ DBLE( AB( KA1, I ) )*
361: $ BB( J-I+KB1, I )*
362: $ DCONJG( 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: $ DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
382: IF( KBT.GT.0 )
383: $ CALL ZGERC( 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 ZLARTG( 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: $ DCONJG( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( 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 ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
482: $ RWORK( J-M ), DCONJG( 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 ZLARTV( 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 ZLARGV( 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 ZLARTV( 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 ZLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
557: $ AB( KA, J2+1 ), INCA, RWORK( J2 ),
558: $ WORK( J2 ), KA1 )
559: *
560: CALL ZLACGV( 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 ZLARTV( 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 ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
579: $ RWORK( J ), DCONJG( 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 ZLARTV( 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 = DBLE( BB( 1, I ) )
614: AB( 1, I ) = ( DBLE( 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 )*DCONJG( AB( I-K+1,
625: $ K ) ) - DCONJG( BB( I-K+1, K ) )*
626: $ AB( I-J+1, J ) + DBLE( AB( 1, I ) )*
627: $ BB( I-J+1, J )*DCONJG( 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: $ DCONJG( 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 ZDSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
648: IF( KBT.GT.0 )
649: $ CALL ZGERU( 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 ZLARTG( 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: $ DCONJG( WORK( I-K+KA-M ) )*
682: $ AB( KA1, I-K )
683: AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
684: $ RWORK( I-K+KA-M )*AB( KA1, I-K )
685: RA1 = RA
686: END IF
687: END IF
688: J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
689: NR = ( N-J2+KA ) / KA1
690: J1 = J2 + ( NR-1 )*KA1
691: IF( UPDATE ) THEN
692: J2T = MAX( J2, I+2*KA-K+1 )
693: ELSE
694: J2T = J2
695: END IF
696: NRT = ( N-J2T+KA ) / KA1
697: DO 320 J = J2T, J1, KA1
698: *
699: * create nonzero element a(j+1,j-ka) outside the band
700: * and store it in WORK(j-m)
701: *
702: WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
703: AB( KA1, J-KA+1 ) = RWORK( J-M )*AB( KA1, J-KA+1 )
704: 320 CONTINUE
705: *
706: * generate rotations in 1st set to annihilate elements which
707: * have been created outside the band
708: *
709: IF( NRT.GT.0 )
710: $ CALL ZLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
711: $ KA1, RWORK( J2T-M ), KA1 )
712: IF( NR.GT.0 ) THEN
713: *
714: * apply rotations in 1st set from the left
715: *
716: DO 330 L = 1, KA - 1
717: CALL ZLARTV( NR, AB( L+1, J2-L ), INCA,
718: $ AB( L+2, J2-L ), INCA, RWORK( J2-M ),
719: $ WORK( J2-M ), KA1 )
720: 330 CONTINUE
721: *
722: * apply rotations in 1st set from both sides to diagonal
723: * blocks
724: *
725: CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
726: $ INCA, RWORK( J2-M ), WORK( J2-M ), KA1 )
727: *
728: CALL ZLACGV( NR, WORK( J2-M ), KA1 )
729: END IF
730: *
731: * start applying rotations in 1st set from the right
732: *
733: DO 340 L = KA - 1, KB - K + 1, -1
734: NRT = ( N-J2+L ) / KA1
735: IF( NRT.GT.0 )
736: $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
737: $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
738: $ WORK( J2-M ), KA1 )
739: 340 CONTINUE
740: *
741: IF( WANTX ) THEN
742: *
743: * post-multiply X by product of rotations in 1st set
744: *
745: DO 350 J = J2, J1, KA1
746: CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
747: $ RWORK( J-M ), WORK( J-M ) )
748: 350 CONTINUE
749: END IF
750: 360 CONTINUE
751: *
752: IF( UPDATE ) THEN
753: IF( I2.LE.N .AND. KBT.GT.0 ) THEN
754: *
755: * create nonzero element a(i-kbt+ka+1,i-kbt) outside the
756: * band and store it in WORK(i-kbt)
757: *
758: WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
759: END IF
760: END IF
761: *
762: DO 400 K = KB, 1, -1
763: IF( UPDATE ) THEN
764: J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
765: ELSE
766: J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
767: END IF
768: *
769: * finish applying rotations in 2nd set from the right
770: *
771: DO 370 L = KB - K, 1, -1
772: NRT = ( N-J2+KA+L ) / KA1
773: IF( NRT.GT.0 )
774: $ CALL ZLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
775: $ AB( KA1-L, J2-KA+1 ), INCA,
776: $ RWORK( J2-KA ), WORK( J2-KA ), KA1 )
777: 370 CONTINUE
778: NR = ( N-J2+KA ) / KA1
779: J1 = J2 + ( NR-1 )*KA1
780: DO 380 J = J1, J2, -KA1
781: WORK( J ) = WORK( J-KA )
782: RWORK( J ) = RWORK( J-KA )
783: 380 CONTINUE
784: DO 390 J = J2, J1, KA1
785: *
786: * create nonzero element a(j+1,j-ka) outside the band
787: * and store it in WORK(j)
788: *
789: WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
790: AB( KA1, J-KA+1 ) = RWORK( J )*AB( KA1, J-KA+1 )
791: 390 CONTINUE
792: IF( UPDATE ) THEN
793: IF( I-K.LT.N-KA .AND. K.LE.KBT )
794: $ WORK( I-K+KA ) = WORK( I-K )
795: END IF
796: 400 CONTINUE
797: *
798: DO 440 K = KB, 1, -1
799: J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
800: NR = ( N-J2+KA ) / KA1
801: J1 = J2 + ( NR-1 )*KA1
802: IF( NR.GT.0 ) THEN
803: *
804: * generate rotations in 2nd set to annihilate elements
805: * which have been created outside the band
806: *
807: CALL ZLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
808: $ RWORK( J2 ), KA1 )
809: *
810: * apply rotations in 2nd set from the left
811: *
812: DO 410 L = 1, KA - 1
813: CALL ZLARTV( NR, AB( L+1, J2-L ), INCA,
814: $ AB( L+2, J2-L ), INCA, RWORK( J2 ),
815: $ WORK( J2 ), KA1 )
816: 410 CONTINUE
817: *
818: * apply rotations in 2nd set from both sides to diagonal
819: * blocks
820: *
821: CALL ZLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
822: $ INCA, RWORK( J2 ), WORK( J2 ), KA1 )
823: *
824: CALL ZLACGV( NR, WORK( J2 ), KA1 )
825: END IF
826: *
827: * start applying rotations in 2nd set from the right
828: *
829: DO 420 L = KA - 1, KB - K + 1, -1
830: NRT = ( N-J2+L ) / KA1
831: IF( NRT.GT.0 )
832: $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
833: $ AB( KA1-L, J2+1 ), INCA, RWORK( J2 ),
834: $ WORK( J2 ), KA1 )
835: 420 CONTINUE
836: *
837: IF( WANTX ) THEN
838: *
839: * post-multiply X by product of rotations in 2nd set
840: *
841: DO 430 J = J2, J1, KA1
842: CALL ZROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
843: $ RWORK( J ), WORK( J ) )
844: 430 CONTINUE
845: END IF
846: 440 CONTINUE
847: *
848: DO 460 K = 1, KB - 1
849: J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
850: *
851: * finish applying rotations in 1st set from the right
852: *
853: DO 450 L = KB - K, 1, -1
854: NRT = ( N-J2+L ) / KA1
855: IF( NRT.GT.0 )
856: $ CALL ZLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
857: $ AB( KA1-L, J2+1 ), INCA, RWORK( J2-M ),
858: $ WORK( J2-M ), KA1 )
859: 450 CONTINUE
860: 460 CONTINUE
861: *
862: IF( KB.GT.1 ) THEN
863: DO 470 J = N - 1, J2 + KA, -1
864: RWORK( J-M ) = RWORK( J-KA-M )
865: WORK( J-M ) = WORK( J-KA-M )
866: 470 CONTINUE
867: END IF
868: *
869: END IF
870: *
871: GO TO 10
872: *
873: 480 CONTINUE
874: *
875: * **************************** Phase 2 *****************************
876: *
877: * The logical structure of this phase is:
878: *
879: * UPDATE = .TRUE.
880: * DO I = 1, M
881: * use S(i) to update A and create a new bulge
882: * apply rotations to push all bulges KA positions upward
883: * END DO
884: * UPDATE = .FALSE.
885: * DO I = M - KA - 1, 2, -1
886: * apply rotations to push all bulges KA positions upward
887: * END DO
888: *
889: * To avoid duplicating code, the two loops are merged.
890: *
891: UPDATE = .TRUE.
892: I = 0
893: 490 CONTINUE
894: IF( UPDATE ) THEN
895: I = I + 1
896: KBT = MIN( KB, M-I )
897: I0 = I + 1
898: I1 = MAX( 1, I-KA )
899: I2 = I + KBT - KA1
900: IF( I.GT.M ) THEN
901: UPDATE = .FALSE.
902: I = I - 1
903: I0 = M + 1
904: IF( KA.EQ.0 )
905: $ RETURN
906: GO TO 490
907: END IF
908: ELSE
909: I = I - KA
910: IF( I.LT.2 )
911: $ RETURN
912: END IF
913: *
914: IF( I.LT.M-KBT ) THEN
915: NX = M
916: ELSE
917: NX = N
918: END IF
919: *
920: IF( UPPER ) THEN
921: *
922: * Transform A, working with the upper triangle
923: *
924: IF( UPDATE ) THEN
925: *
926: * Form inv(S(i))**H * A * inv(S(i))
927: *
928: BII = DBLE( BB( KB1, I ) )
929: AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII
930: DO 500 J = I1, I - 1
931: AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
932: 500 CONTINUE
933: DO 510 J = I + 1, MIN( N, I+KA )
934: AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
935: 510 CONTINUE
936: DO 540 K = I + 1, I + KBT
937: DO 520 J = K, I + KBT
938: AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
939: $ BB( I-J+KB1, J )*
940: $ DCONJG( AB( I-K+KA1, K ) ) -
941: $ DCONJG( BB( I-K+KB1, K ) )*
942: $ AB( I-J+KA1, J ) +
943: $ DBLE( AB( KA1, I ) )*
944: $ BB( I-J+KB1, J )*
945: $ DCONJG( BB( I-K+KB1, K ) )
946: 520 CONTINUE
947: DO 530 J = I + KBT + 1, MIN( N, I+KA )
948: AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
949: $ DCONJG( BB( I-K+KB1, K ) )*
950: $ AB( I-J+KA1, J )
951: 530 CONTINUE
952: 540 CONTINUE
953: DO 560 J = I1, I
954: DO 550 K = I + 1, MIN( J+KA, I+KBT )
955: AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
956: $ BB( I-K+KB1, K )*AB( J-I+KA1, I )
957: 550 CONTINUE
958: 560 CONTINUE
959: *
960: IF( WANTX ) THEN
961: *
962: * post-multiply X by inv(S(i))
963: *
964: CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
965: IF( KBT.GT.0 )
966: $ CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1,
967: $ BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX )
968: END IF
969: *
970: * store a(i1,i) in RA1 for use in next loop over K
971: *
972: RA1 = AB( I1-I+KA1, I )
973: END IF
974: *
975: * Generate and apply vectors of rotations to chase all the
976: * existing bulges KA positions up toward the top of the band
977: *
978: DO 610 K = 1, KB - 1
979: IF( UPDATE ) THEN
980: *
981: * Determine the rotations which would annihilate the bulge
982: * which has in theory just been created
983: *
984: IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
985: *
986: * generate rotation to annihilate a(i+k-ka-1,i)
987: *
988: CALL ZLARTG( AB( K+1, I ), RA1, RWORK( I+K-KA ),
989: $ WORK( I+K-KA ), RA )
990: *
991: * create nonzero element a(i+k-ka-1,i+k) outside the
992: * band and store it in WORK(m-kb+i+k)
993: *
994: T = -BB( KB1-K, I+K )*RA1
995: WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
996: $ DCONJG( WORK( I+K-KA ) )*
997: $ AB( 1, I+K )
998: AB( 1, I+K ) = WORK( I+K-KA )*T +
999: $ RWORK( I+K-KA )*AB( 1, I+K )
1000: RA1 = RA
1001: END IF
1002: END IF
1003: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1004: NR = ( J2+KA-1 ) / KA1
1005: J1 = J2 - ( NR-1 )*KA1
1006: IF( UPDATE ) THEN
1007: J2T = MIN( J2, I-2*KA+K-1 )
1008: ELSE
1009: J2T = J2
1010: END IF
1011: NRT = ( J2T+KA-1 ) / KA1
1012: DO 570 J = J1, J2T, KA1
1013: *
1014: * create nonzero element a(j-1,j+ka) outside the band
1015: * and store it in WORK(j)
1016: *
1017: WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
1018: AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 )
1019: 570 CONTINUE
1020: *
1021: * generate rotations in 1st set to annihilate elements which
1022: * have been created outside the band
1023: *
1024: IF( NRT.GT.0 )
1025: $ CALL ZLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
1026: $ RWORK( J1 ), KA1 )
1027: IF( NR.GT.0 ) THEN
1028: *
1029: * apply rotations in 1st set from the left
1030: *
1031: DO 580 L = 1, KA - 1
1032: CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA,
1033: $ AB( KA-L, J1+L ), INCA, RWORK( J1 ),
1034: $ WORK( J1 ), KA1 )
1035: 580 CONTINUE
1036: *
1037: * apply rotations in 1st set from both sides to diagonal
1038: * blocks
1039: *
1040: CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1041: $ AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ),
1042: $ KA1 )
1043: *
1044: CALL ZLACGV( NR, WORK( J1 ), KA1 )
1045: END IF
1046: *
1047: * start applying rotations in 1st set from the right
1048: *
1049: DO 590 L = KA - 1, KB - K + 1, -1
1050: NRT = ( J2+L-1 ) / KA1
1051: J1T = J2 - ( NRT-1 )*KA1
1052: IF( NRT.GT.0 )
1053: $ CALL ZLARTV( NRT, AB( L, J1T ), INCA,
1054: $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
1055: $ WORK( J1T ), KA1 )
1056: 590 CONTINUE
1057: *
1058: IF( WANTX ) THEN
1059: *
1060: * post-multiply X by product of rotations in 1st set
1061: *
1062: DO 600 J = J1, J2, KA1
1063: CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1064: $ RWORK( J ), WORK( J ) )
1065: 600 CONTINUE
1066: END IF
1067: 610 CONTINUE
1068: *
1069: IF( UPDATE ) THEN
1070: IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1071: *
1072: * create nonzero element a(i+kbt-ka-1,i+kbt) outside the
1073: * band and store it in WORK(m-kb+i+kbt)
1074: *
1075: WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
1076: END IF
1077: END IF
1078: *
1079: DO 650 K = KB, 1, -1
1080: IF( UPDATE ) THEN
1081: J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1082: ELSE
1083: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1084: END IF
1085: *
1086: * finish applying rotations in 2nd set from the right
1087: *
1088: DO 620 L = KB - K, 1, -1
1089: NRT = ( J2+KA+L-1 ) / KA1
1090: J1T = J2 - ( NRT-1 )*KA1
1091: IF( NRT.GT.0 )
1092: $ CALL ZLARTV( NRT, AB( L, J1T+KA ), INCA,
1093: $ AB( L+1, J1T+KA-1 ), INCA,
1094: $ RWORK( M-KB+J1T+KA ),
1095: $ WORK( M-KB+J1T+KA ), KA1 )
1096: 620 CONTINUE
1097: NR = ( J2+KA-1 ) / KA1
1098: J1 = J2 - ( NR-1 )*KA1
1099: DO 630 J = J1, J2, KA1
1100: WORK( M-KB+J ) = WORK( M-KB+J+KA )
1101: RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
1102: 630 CONTINUE
1103: DO 640 J = J1, J2, KA1
1104: *
1105: * create nonzero element a(j-1,j+ka) outside the band
1106: * and store it in WORK(m-kb+j)
1107: *
1108: WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
1109: AB( 1, J+KA-1 ) = RWORK( M-KB+J )*AB( 1, J+KA-1 )
1110: 640 CONTINUE
1111: IF( UPDATE ) THEN
1112: IF( I+K.GT.KA1 .AND. K.LE.KBT )
1113: $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1114: END IF
1115: 650 CONTINUE
1116: *
1117: DO 690 K = KB, 1, -1
1118: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1119: NR = ( J2+KA-1 ) / KA1
1120: J1 = J2 - ( NR-1 )*KA1
1121: IF( NR.GT.0 ) THEN
1122: *
1123: * generate rotations in 2nd set to annihilate elements
1124: * which have been created outside the band
1125: *
1126: CALL ZLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
1127: $ KA1, RWORK( M-KB+J1 ), KA1 )
1128: *
1129: * apply rotations in 2nd set from the left
1130: *
1131: DO 660 L = 1, KA - 1
1132: CALL ZLARTV( NR, AB( KA1-L, J1+L ), INCA,
1133: $ AB( KA-L, J1+L ), INCA, RWORK( M-KB+J1 ),
1134: $ WORK( M-KB+J1 ), KA1 )
1135: 660 CONTINUE
1136: *
1137: * apply rotations in 2nd set from both sides to diagonal
1138: * blocks
1139: *
1140: CALL ZLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
1141: $ AB( KA, J1 ), INCA, RWORK( M-KB+J1 ),
1142: $ WORK( M-KB+J1 ), KA1 )
1143: *
1144: CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 )
1145: END IF
1146: *
1147: * start applying rotations in 2nd set from the right
1148: *
1149: DO 670 L = KA - 1, KB - K + 1, -1
1150: NRT = ( J2+L-1 ) / KA1
1151: J1T = J2 - ( NRT-1 )*KA1
1152: IF( NRT.GT.0 )
1153: $ CALL ZLARTV( NRT, AB( L, J1T ), INCA,
1154: $ AB( L+1, J1T-1 ), INCA,
1155: $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
1156: $ KA1 )
1157: 670 CONTINUE
1158: *
1159: IF( WANTX ) THEN
1160: *
1161: * post-multiply X by product of rotations in 2nd set
1162: *
1163: DO 680 J = J1, J2, KA1
1164: CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1165: $ RWORK( M-KB+J ), WORK( M-KB+J ) )
1166: 680 CONTINUE
1167: END IF
1168: 690 CONTINUE
1169: *
1170: DO 710 K = 1, KB - 1
1171: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1172: *
1173: * finish applying rotations in 1st set from the right
1174: *
1175: DO 700 L = KB - K, 1, -1
1176: NRT = ( J2+L-1 ) / KA1
1177: J1T = J2 - ( NRT-1 )*KA1
1178: IF( NRT.GT.0 )
1179: $ CALL ZLARTV( NRT, AB( L, J1T ), INCA,
1180: $ AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
1181: $ WORK( J1T ), KA1 )
1182: 700 CONTINUE
1183: 710 CONTINUE
1184: *
1185: IF( KB.GT.1 ) THEN
1186: DO 720 J = 2, I2 - KA
1187: RWORK( J ) = RWORK( J+KA )
1188: WORK( J ) = WORK( J+KA )
1189: 720 CONTINUE
1190: END IF
1191: *
1192: ELSE
1193: *
1194: * Transform A, working with the lower triangle
1195: *
1196: IF( UPDATE ) THEN
1197: *
1198: * Form inv(S(i))**H * A * inv(S(i))
1199: *
1200: BII = DBLE( BB( 1, I ) )
1201: AB( 1, I ) = ( DBLE( AB( 1, I ) ) / BII ) / BII
1202: DO 730 J = I1, I - 1
1203: AB( I-J+1, J ) = AB( I-J+1, J ) / BII
1204: 730 CONTINUE
1205: DO 740 J = I + 1, MIN( N, I+KA )
1206: AB( J-I+1, I ) = AB( J-I+1, I ) / BII
1207: 740 CONTINUE
1208: DO 770 K = I + 1, I + KBT
1209: DO 750 J = K, I + KBT
1210: AB( J-K+1, K ) = AB( J-K+1, K ) -
1211: $ BB( J-I+1, I )*DCONJG( AB( K-I+1,
1212: $ I ) ) - DCONJG( BB( K-I+1, I ) )*
1213: $ AB( J-I+1, I ) + DBLE( AB( 1, I ) )*
1214: $ BB( J-I+1, I )*DCONJG( BB( K-I+1,
1215: $ I ) )
1216: 750 CONTINUE
1217: DO 760 J = I + KBT + 1, MIN( N, I+KA )
1218: AB( J-K+1, K ) = AB( J-K+1, K ) -
1219: $ DCONJG( BB( K-I+1, I ) )*
1220: $ AB( J-I+1, I )
1221: 760 CONTINUE
1222: 770 CONTINUE
1223: DO 790 J = I1, I
1224: DO 780 K = I + 1, MIN( J+KA, I+KBT )
1225: AB( K-J+1, J ) = AB( K-J+1, J ) -
1226: $ BB( K-I+1, I )*AB( I-J+1, J )
1227: 780 CONTINUE
1228: 790 CONTINUE
1229: *
1230: IF( WANTX ) THEN
1231: *
1232: * post-multiply X by inv(S(i))
1233: *
1234: CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
1235: IF( KBT.GT.0 )
1236: $ CALL ZGERC( NX, KBT, -CONE, X( 1, I ), 1, BB( 2, I ),
1237: $ 1, X( 1, I+1 ), LDX )
1238: END IF
1239: *
1240: * store a(i,i1) in RA1 for use in next loop over K
1241: *
1242: RA1 = AB( I-I1+1, I1 )
1243: END IF
1244: *
1245: * Generate and apply vectors of rotations to chase all the
1246: * existing bulges KA positions up toward the top of the band
1247: *
1248: DO 840 K = 1, KB - 1
1249: IF( UPDATE ) THEN
1250: *
1251: * Determine the rotations which would annihilate the bulge
1252: * which has in theory just been created
1253: *
1254: IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
1255: *
1256: * generate rotation to annihilate a(i,i+k-ka-1)
1257: *
1258: CALL ZLARTG( AB( KA1-K, I+K-KA ), RA1,
1259: $ RWORK( I+K-KA ), WORK( I+K-KA ), RA )
1260: *
1261: * create nonzero element a(i+k,i+k-ka-1) outside the
1262: * band and store it in WORK(m-kb+i+k)
1263: *
1264: T = -BB( K+1, I )*RA1
1265: WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
1266: $ DCONJG( WORK( I+K-KA ) )*
1267: $ AB( KA1, I+K-KA )
1268: AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
1269: $ RWORK( I+K-KA )*AB( KA1, I+K-KA )
1270: RA1 = RA
1271: END IF
1272: END IF
1273: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1274: NR = ( J2+KA-1 ) / KA1
1275: J1 = J2 - ( NR-1 )*KA1
1276: IF( UPDATE ) THEN
1277: J2T = MIN( J2, I-2*KA+K-1 )
1278: ELSE
1279: J2T = J2
1280: END IF
1281: NRT = ( J2T+KA-1 ) / KA1
1282: DO 800 J = J1, J2T, KA1
1283: *
1284: * create nonzero element a(j+ka,j-1) outside the band
1285: * and store it in WORK(j)
1286: *
1287: WORK( J ) = WORK( J )*AB( KA1, J-1 )
1288: AB( KA1, J-1 ) = RWORK( J )*AB( KA1, J-1 )
1289: 800 CONTINUE
1290: *
1291: * generate rotations in 1st set to annihilate elements which
1292: * have been created outside the band
1293: *
1294: IF( NRT.GT.0 )
1295: $ CALL ZLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
1296: $ RWORK( J1 ), KA1 )
1297: IF( NR.GT.0 ) THEN
1298: *
1299: * apply rotations in 1st set from the right
1300: *
1301: DO 810 L = 1, KA - 1
1302: CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1303: $ INCA, RWORK( J1 ), WORK( J1 ), KA1 )
1304: 810 CONTINUE
1305: *
1306: * apply rotations in 1st set from both sides to diagonal
1307: * blocks
1308: *
1309: CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1310: $ AB( 2, J1-1 ), INCA, RWORK( J1 ),
1311: $ WORK( J1 ), KA1 )
1312: *
1313: CALL ZLACGV( NR, WORK( J1 ), KA1 )
1314: END IF
1315: *
1316: * start applying rotations in 1st set from the left
1317: *
1318: DO 820 L = KA - 1, KB - K + 1, -1
1319: NRT = ( J2+L-1 ) / KA1
1320: J1T = J2 - ( NRT-1 )*KA1
1321: IF( NRT.GT.0 )
1322: $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1323: $ AB( KA1-L, J1T-KA1+L ), INCA,
1324: $ RWORK( J1T ), WORK( J1T ), KA1 )
1325: 820 CONTINUE
1326: *
1327: IF( WANTX ) THEN
1328: *
1329: * post-multiply X by product of rotations in 1st set
1330: *
1331: DO 830 J = J1, J2, KA1
1332: CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1333: $ RWORK( J ), DCONJG( WORK( J ) ) )
1334: 830 CONTINUE
1335: END IF
1336: 840 CONTINUE
1337: *
1338: IF( UPDATE ) THEN
1339: IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
1340: *
1341: * create nonzero element a(i+kbt,i+kbt-ka-1) outside the
1342: * band and store it in WORK(m-kb+i+kbt)
1343: *
1344: WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
1345: END IF
1346: END IF
1347: *
1348: DO 880 K = KB, 1, -1
1349: IF( UPDATE ) THEN
1350: J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
1351: ELSE
1352: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1353: END IF
1354: *
1355: * finish applying rotations in 2nd set from the left
1356: *
1357: DO 850 L = KB - K, 1, -1
1358: NRT = ( J2+KA+L-1 ) / KA1
1359: J1T = J2 - ( NRT-1 )*KA1
1360: IF( NRT.GT.0 )
1361: $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
1362: $ AB( KA1-L, J1T+L-1 ), INCA,
1363: $ RWORK( M-KB+J1T+KA ),
1364: $ WORK( M-KB+J1T+KA ), KA1 )
1365: 850 CONTINUE
1366: NR = ( J2+KA-1 ) / KA1
1367: J1 = J2 - ( NR-1 )*KA1
1368: DO 860 J = J1, J2, KA1
1369: WORK( M-KB+J ) = WORK( M-KB+J+KA )
1370: RWORK( M-KB+J ) = RWORK( M-KB+J+KA )
1371: 860 CONTINUE
1372: DO 870 J = J1, J2, KA1
1373: *
1374: * create nonzero element a(j+ka,j-1) outside the band
1375: * and store it in WORK(m-kb+j)
1376: *
1377: WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
1378: AB( KA1, J-1 ) = RWORK( M-KB+J )*AB( KA1, J-1 )
1379: 870 CONTINUE
1380: IF( UPDATE ) THEN
1381: IF( I+K.GT.KA1 .AND. K.LE.KBT )
1382: $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
1383: END IF
1384: 880 CONTINUE
1385: *
1386: DO 920 K = KB, 1, -1
1387: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
1388: NR = ( J2+KA-1 ) / KA1
1389: J1 = J2 - ( NR-1 )*KA1
1390: IF( NR.GT.0 ) THEN
1391: *
1392: * generate rotations in 2nd set to annihilate elements
1393: * which have been created outside the band
1394: *
1395: CALL ZLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
1396: $ KA1, RWORK( M-KB+J1 ), KA1 )
1397: *
1398: * apply rotations in 2nd set from the right
1399: *
1400: DO 890 L = 1, KA - 1
1401: CALL ZLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
1402: $ INCA, RWORK( M-KB+J1 ), WORK( M-KB+J1 ),
1403: $ KA1 )
1404: 890 CONTINUE
1405: *
1406: * apply rotations in 2nd set from both sides to diagonal
1407: * blocks
1408: *
1409: CALL ZLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
1410: $ AB( 2, J1-1 ), INCA, RWORK( M-KB+J1 ),
1411: $ WORK( M-KB+J1 ), KA1 )
1412: *
1413: CALL ZLACGV( NR, WORK( M-KB+J1 ), KA1 )
1414: END IF
1415: *
1416: * start applying rotations in 2nd set from the left
1417: *
1418: DO 900 L = KA - 1, KB - K + 1, -1
1419: NRT = ( J2+L-1 ) / KA1
1420: J1T = J2 - ( NRT-1 )*KA1
1421: IF( NRT.GT.0 )
1422: $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1423: $ AB( KA1-L, J1T-KA1+L ), INCA,
1424: $ RWORK( M-KB+J1T ), WORK( M-KB+J1T ),
1425: $ KA1 )
1426: 900 CONTINUE
1427: *
1428: IF( WANTX ) THEN
1429: *
1430: * post-multiply X by product of rotations in 2nd set
1431: *
1432: DO 910 J = J1, J2, KA1
1433: CALL ZROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
1434: $ RWORK( M-KB+J ), DCONJG( WORK( M-KB+J ) ) )
1435: 910 CONTINUE
1436: END IF
1437: 920 CONTINUE
1438: *
1439: DO 940 K = 1, KB - 1
1440: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
1441: *
1442: * finish applying rotations in 1st set from the left
1443: *
1444: DO 930 L = KB - K, 1, -1
1445: NRT = ( J2+L-1 ) / KA1
1446: J1T = J2 - ( NRT-1 )*KA1
1447: IF( NRT.GT.0 )
1448: $ CALL ZLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
1449: $ AB( KA1-L, J1T-KA1+L ), INCA,
1450: $ RWORK( J1T ), WORK( J1T ), KA1 )
1451: 930 CONTINUE
1452: 940 CONTINUE
1453: *
1454: IF( KB.GT.1 ) THEN
1455: DO 950 J = 2, I2 - KA
1456: RWORK( J ) = RWORK( J+KA )
1457: WORK( J ) = WORK( J+KA )
1458: 950 CONTINUE
1459: END IF
1460: *
1461: END IF
1462: *
1463: GO TO 490
1464: *
1465: * End of ZHBGST
1466: *
1467: 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