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