1: SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
2: $ WORK, 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, KD, LDAB, LDQ, N
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION D( * ), E( * )
15: COMPLEX*16 AB( LDAB, * ), Q( LDQ, * ), WORK( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * ZHBTRD reduces a complex Hermitian band matrix A to real symmetric
22: * tridiagonal form T by a unitary similarity transformation:
23: * Q**H * A * Q = T.
24: *
25: * Arguments
26: * =========
27: *
28: * VECT (input) CHARACTER*1
29: * = 'N': do not form Q;
30: * = 'V': form Q;
31: * = 'U': update a matrix X, by forming X*Q.
32: *
33: * UPLO (input) CHARACTER*1
34: * = 'U': Upper triangle of A is stored;
35: * = 'L': Lower triangle of A is stored.
36: *
37: * N (input) INTEGER
38: * The order of the matrix A. N >= 0.
39: *
40: * KD (input) INTEGER
41: * The number of superdiagonals of the matrix A if UPLO = 'U',
42: * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
43: *
44: * AB (input/output) COMPLEX*16 array, dimension (LDAB,N)
45: * On entry, the upper or lower triangle of the Hermitian band
46: * matrix A, stored in the first KD+1 rows of the array. The
47: * j-th column of A is stored in the j-th column of the array AB
48: * as follows:
49: * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
50: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
51: * On exit, the diagonal elements of AB are overwritten by the
52: * diagonal elements of the tridiagonal matrix T; if KD > 0, the
53: * elements on the first superdiagonal (if UPLO = 'U') or the
54: * first subdiagonal (if UPLO = 'L') are overwritten by the
55: * off-diagonal elements of T; the rest of AB is overwritten by
56: * values generated during the reduction.
57: *
58: * LDAB (input) INTEGER
59: * The leading dimension of the array AB. LDAB >= KD+1.
60: *
61: * D (output) DOUBLE PRECISION array, dimension (N)
62: * The diagonal elements of the tridiagonal matrix T.
63: *
64: * E (output) DOUBLE PRECISION array, dimension (N-1)
65: * The off-diagonal elements of the tridiagonal matrix T:
66: * E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
67: *
68: * Q (input/output) COMPLEX*16 array, dimension (LDQ,N)
69: * On entry, if VECT = 'U', then Q must contain an N-by-N
70: * matrix X; if VECT = 'N' or 'V', then Q need not be set.
71: *
72: * On exit:
73: * if VECT = 'V', Q contains the N-by-N unitary matrix Q;
74: * if VECT = 'U', Q contains the product X*Q;
75: * if VECT = 'N', the array Q is not referenced.
76: *
77: * LDQ (input) INTEGER
78: * The leading dimension of the array Q.
79: * LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
80: *
81: * WORK (workspace) COMPLEX*16 array, dimension (N)
82: *
83: * INFO (output) INTEGER
84: * = 0: successful exit
85: * < 0: if INFO = -i, the i-th argument had an illegal value
86: *
87: * Further Details
88: * ===============
89: *
90: * Modified by Linda Kaufman, Bell Labs.
91: *
92: * =====================================================================
93: *
94: * .. Parameters ..
95: DOUBLE PRECISION ZERO
96: PARAMETER ( ZERO = 0.0D+0 )
97: COMPLEX*16 CZERO, CONE
98: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
99: $ CONE = ( 1.0D+0, 0.0D+0 ) )
100: * ..
101: * .. Local Scalars ..
102: LOGICAL INITQ, UPPER, WANTQ
103: INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
104: $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
105: $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
106: DOUBLE PRECISION ABST
107: COMPLEX*16 T, TEMP
108: * ..
109: * .. External Subroutines ..
110: EXTERNAL XERBLA, ZLACGV, ZLAR2V, ZLARGV, ZLARTG, ZLARTV,
111: $ ZLASET, ZROT, ZSCAL
112: * ..
113: * .. Intrinsic Functions ..
114: INTRINSIC ABS, DBLE, DCONJG, MAX, MIN
115: * ..
116: * .. External Functions ..
117: LOGICAL LSAME
118: EXTERNAL LSAME
119: * ..
120: * .. Executable Statements ..
121: *
122: * Test the input parameters
123: *
124: INITQ = LSAME( VECT, 'V' )
125: WANTQ = INITQ .OR. LSAME( VECT, 'U' )
126: UPPER = LSAME( UPLO, 'U' )
127: KD1 = KD + 1
128: KDM1 = KD - 1
129: INCX = LDAB - 1
130: IQEND = 1
131: *
132: INFO = 0
133: IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
134: INFO = -1
135: ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
136: INFO = -2
137: ELSE IF( N.LT.0 ) THEN
138: INFO = -3
139: ELSE IF( KD.LT.0 ) THEN
140: INFO = -4
141: ELSE IF( LDAB.LT.KD1 ) THEN
142: INFO = -6
143: ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
144: INFO = -10
145: END IF
146: IF( INFO.NE.0 ) THEN
147: CALL XERBLA( 'ZHBTRD', -INFO )
148: RETURN
149: END IF
150: *
151: * Quick return if possible
152: *
153: IF( N.EQ.0 )
154: $ RETURN
155: *
156: * Initialize Q to the unit matrix, if needed
157: *
158: IF( INITQ )
159: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
160: *
161: * Wherever possible, plane rotations are generated and applied in
162: * vector operations of length NR over the index set J1:J2:KD1.
163: *
164: * The real cosines and complex sines of the plane rotations are
165: * stored in the arrays D and WORK.
166: *
167: INCA = KD1*LDAB
168: KDN = MIN( N-1, KD )
169: IF( UPPER ) THEN
170: *
171: IF( KD.GT.1 ) THEN
172: *
173: * Reduce to complex Hermitian tridiagonal form, working with
174: * the upper triangle
175: *
176: NR = 0
177: J1 = KDN + 2
178: J2 = 1
179: *
180: AB( KD1, 1 ) = DBLE( AB( KD1, 1 ) )
181: DO 90 I = 1, N - 2
182: *
183: * Reduce i-th row of matrix to tridiagonal form
184: *
185: DO 80 K = KDN + 1, 2, -1
186: J1 = J1 + KDN
187: J2 = J2 + KDN
188: *
189: IF( NR.GT.0 ) THEN
190: *
191: * generate plane rotations to annihilate nonzero
192: * elements which have been created outside the band
193: *
194: CALL ZLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
195: $ KD1, D( J1 ), KD1 )
196: *
197: * apply rotations from the right
198: *
199: *
200: * Dependent on the the number of diagonals either
201: * ZLARTV or ZROT is used
202: *
203: IF( NR.GE.2*KD-1 ) THEN
204: DO 10 L = 1, KD - 1
205: CALL ZLARTV( NR, AB( L+1, J1-1 ), INCA,
206: $ AB( L, J1 ), INCA, D( J1 ),
207: $ WORK( J1 ), KD1 )
208: 10 CONTINUE
209: *
210: ELSE
211: JEND = J1 + ( NR-1 )*KD1
212: DO 20 JINC = J1, JEND, KD1
213: CALL ZROT( KDM1, AB( 2, JINC-1 ), 1,
214: $ AB( 1, JINC ), 1, D( JINC ),
215: $ WORK( JINC ) )
216: 20 CONTINUE
217: END IF
218: END IF
219: *
220: *
221: IF( K.GT.2 ) THEN
222: IF( K.LE.N-I+1 ) THEN
223: *
224: * generate plane rotation to annihilate a(i,i+k-1)
225: * within the band
226: *
227: CALL ZLARTG( AB( KD-K+3, I+K-2 ),
228: $ AB( KD-K+2, I+K-1 ), D( I+K-1 ),
229: $ WORK( I+K-1 ), TEMP )
230: AB( KD-K+3, I+K-2 ) = TEMP
231: *
232: * apply rotation from the right
233: *
234: CALL ZROT( K-3, AB( KD-K+4, I+K-2 ), 1,
235: $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
236: $ WORK( I+K-1 ) )
237: END IF
238: NR = NR + 1
239: J1 = J1 - KDN - 1
240: END IF
241: *
242: * apply plane rotations from both sides to diagonal
243: * blocks
244: *
245: IF( NR.GT.0 )
246: $ CALL ZLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
247: $ AB( KD, J1 ), INCA, D( J1 ),
248: $ WORK( J1 ), KD1 )
249: *
250: * apply plane rotations from the left
251: *
252: IF( NR.GT.0 ) THEN
253: CALL ZLACGV( NR, WORK( J1 ), KD1 )
254: IF( 2*KD-1.LT.NR ) THEN
255: *
256: * Dependent on the the number of diagonals either
257: * ZLARTV or ZROT is used
258: *
259: DO 30 L = 1, KD - 1
260: IF( J2+L.GT.N ) THEN
261: NRT = NR - 1
262: ELSE
263: NRT = NR
264: END IF
265: IF( NRT.GT.0 )
266: $ CALL ZLARTV( NRT, AB( KD-L, J1+L ), INCA,
267: $ AB( KD-L+1, J1+L ), INCA,
268: $ D( J1 ), WORK( J1 ), KD1 )
269: 30 CONTINUE
270: ELSE
271: J1END = J1 + KD1*( NR-2 )
272: IF( J1END.GE.J1 ) THEN
273: DO 40 JIN = J1, J1END, KD1
274: CALL ZROT( KD-1, AB( KD-1, JIN+1 ), INCX,
275: $ AB( KD, JIN+1 ), INCX,
276: $ D( JIN ), WORK( JIN ) )
277: 40 CONTINUE
278: END IF
279: LEND = MIN( KDM1, N-J2 )
280: LAST = J1END + KD1
281: IF( LEND.GT.0 )
282: $ CALL ZROT( LEND, AB( KD-1, LAST+1 ), INCX,
283: $ AB( KD, LAST+1 ), INCX, D( LAST ),
284: $ WORK( LAST ) )
285: END IF
286: END IF
287: *
288: IF( WANTQ ) THEN
289: *
290: * accumulate product of plane rotations in Q
291: *
292: IF( INITQ ) THEN
293: *
294: * take advantage of the fact that Q was
295: * initially the Identity matrix
296: *
297: IQEND = MAX( IQEND, J2 )
298: I2 = MAX( 0, K-3 )
299: IQAEND = 1 + I*KD
300: IF( K.EQ.2 )
301: $ IQAEND = IQAEND + KD
302: IQAEND = MIN( IQAEND, IQEND )
303: DO 50 J = J1, J2, KD1
304: IBL = I - I2 / KDM1
305: I2 = I2 + 1
306: IQB = MAX( 1, J-IBL )
307: NQ = 1 + IQAEND - IQB
308: IQAEND = MIN( IQAEND+KD, IQEND )
309: CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
310: $ 1, D( J ), DCONJG( WORK( J ) ) )
311: 50 CONTINUE
312: ELSE
313: *
314: DO 60 J = J1, J2, KD1
315: CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
316: $ D( J ), DCONJG( WORK( J ) ) )
317: 60 CONTINUE
318: END IF
319: *
320: END IF
321: *
322: IF( J2+KDN.GT.N ) THEN
323: *
324: * adjust J2 to keep within the bounds of the matrix
325: *
326: NR = NR - 1
327: J2 = J2 - KDN - 1
328: END IF
329: *
330: DO 70 J = J1, J2, KD1
331: *
332: * create nonzero element a(j-1,j+kd) outside the band
333: * and store it in WORK
334: *
335: WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
336: AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
337: 70 CONTINUE
338: 80 CONTINUE
339: 90 CONTINUE
340: END IF
341: *
342: IF( KD.GT.0 ) THEN
343: *
344: * make off-diagonal elements real and copy them to E
345: *
346: DO 100 I = 1, N - 1
347: T = AB( KD, I+1 )
348: ABST = ABS( T )
349: AB( KD, I+1 ) = ABST
350: E( I ) = ABST
351: IF( ABST.NE.ZERO ) THEN
352: T = T / ABST
353: ELSE
354: T = CONE
355: END IF
356: IF( I.LT.N-1 )
357: $ AB( KD, I+2 ) = AB( KD, I+2 )*T
358: IF( WANTQ ) THEN
359: CALL ZSCAL( N, DCONJG( T ), Q( 1, I+1 ), 1 )
360: END IF
361: 100 CONTINUE
362: ELSE
363: *
364: * set E to zero if original matrix was diagonal
365: *
366: DO 110 I = 1, N - 1
367: E( I ) = ZERO
368: 110 CONTINUE
369: END IF
370: *
371: * copy diagonal elements to D
372: *
373: DO 120 I = 1, N
374: D( I ) = AB( KD1, I )
375: 120 CONTINUE
376: *
377: ELSE
378: *
379: IF( KD.GT.1 ) THEN
380: *
381: * Reduce to complex Hermitian tridiagonal form, working with
382: * the lower triangle
383: *
384: NR = 0
385: J1 = KDN + 2
386: J2 = 1
387: *
388: AB( 1, 1 ) = DBLE( AB( 1, 1 ) )
389: DO 210 I = 1, N - 2
390: *
391: * Reduce i-th column of matrix to tridiagonal form
392: *
393: DO 200 K = KDN + 1, 2, -1
394: J1 = J1 + KDN
395: J2 = J2 + KDN
396: *
397: IF( NR.GT.0 ) THEN
398: *
399: * generate plane rotations to annihilate nonzero
400: * elements which have been created outside the band
401: *
402: CALL ZLARGV( NR, AB( KD1, J1-KD1 ), INCA,
403: $ WORK( J1 ), KD1, D( J1 ), KD1 )
404: *
405: * apply plane rotations from one side
406: *
407: *
408: * Dependent on the the number of diagonals either
409: * ZLARTV or ZROT is used
410: *
411: IF( NR.GT.2*KD-1 ) THEN
412: DO 130 L = 1, KD - 1
413: CALL ZLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
414: $ AB( KD1-L+1, J1-KD1+L ), INCA,
415: $ D( J1 ), WORK( J1 ), KD1 )
416: 130 CONTINUE
417: ELSE
418: JEND = J1 + KD1*( NR-1 )
419: DO 140 JINC = J1, JEND, KD1
420: CALL ZROT( KDM1, AB( KD, JINC-KD ), INCX,
421: $ AB( KD1, JINC-KD ), INCX,
422: $ D( JINC ), WORK( JINC ) )
423: 140 CONTINUE
424: END IF
425: *
426: END IF
427: *
428: IF( K.GT.2 ) THEN
429: IF( K.LE.N-I+1 ) THEN
430: *
431: * generate plane rotation to annihilate a(i+k-1,i)
432: * within the band
433: *
434: CALL ZLARTG( AB( K-1, I ), AB( K, I ),
435: $ D( I+K-1 ), WORK( I+K-1 ), TEMP )
436: AB( K-1, I ) = TEMP
437: *
438: * apply rotation from the left
439: *
440: CALL ZROT( K-3, AB( K-2, I+1 ), LDAB-1,
441: $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
442: $ WORK( I+K-1 ) )
443: END IF
444: NR = NR + 1
445: J1 = J1 - KDN - 1
446: END IF
447: *
448: * apply plane rotations from both sides to diagonal
449: * blocks
450: *
451: IF( NR.GT.0 )
452: $ CALL ZLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
453: $ AB( 2, J1-1 ), INCA, D( J1 ),
454: $ WORK( J1 ), KD1 )
455: *
456: * apply plane rotations from the right
457: *
458: *
459: * Dependent on the the number of diagonals either
460: * ZLARTV or ZROT is used
461: *
462: IF( NR.GT.0 ) THEN
463: CALL ZLACGV( NR, WORK( J1 ), KD1 )
464: IF( NR.GT.2*KD-1 ) THEN
465: DO 150 L = 1, KD - 1
466: IF( J2+L.GT.N ) THEN
467: NRT = NR - 1
468: ELSE
469: NRT = NR
470: END IF
471: IF( NRT.GT.0 )
472: $ CALL ZLARTV( NRT, AB( L+2, J1-1 ), INCA,
473: $ AB( L+1, J1 ), INCA, D( J1 ),
474: $ WORK( J1 ), KD1 )
475: 150 CONTINUE
476: ELSE
477: J1END = J1 + KD1*( NR-2 )
478: IF( J1END.GE.J1 ) THEN
479: DO 160 J1INC = J1, J1END, KD1
480: CALL ZROT( KDM1, AB( 3, J1INC-1 ), 1,
481: $ AB( 2, J1INC ), 1, D( J1INC ),
482: $ WORK( J1INC ) )
483: 160 CONTINUE
484: END IF
485: LEND = MIN( KDM1, N-J2 )
486: LAST = J1END + KD1
487: IF( LEND.GT.0 )
488: $ CALL ZROT( LEND, AB( 3, LAST-1 ), 1,
489: $ AB( 2, LAST ), 1, D( LAST ),
490: $ WORK( LAST ) )
491: END IF
492: END IF
493: *
494: *
495: *
496: IF( WANTQ ) THEN
497: *
498: * accumulate product of plane rotations in Q
499: *
500: IF( INITQ ) THEN
501: *
502: * take advantage of the fact that Q was
503: * initially the Identity matrix
504: *
505: IQEND = MAX( IQEND, J2 )
506: I2 = MAX( 0, K-3 )
507: IQAEND = 1 + I*KD
508: IF( K.EQ.2 )
509: $ IQAEND = IQAEND + KD
510: IQAEND = MIN( IQAEND, IQEND )
511: DO 170 J = J1, J2, KD1
512: IBL = I - I2 / KDM1
513: I2 = I2 + 1
514: IQB = MAX( 1, J-IBL )
515: NQ = 1 + IQAEND - IQB
516: IQAEND = MIN( IQAEND+KD, IQEND )
517: CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
518: $ 1, D( J ), WORK( J ) )
519: 170 CONTINUE
520: ELSE
521: *
522: DO 180 J = J1, J2, KD1
523: CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
524: $ D( J ), WORK( J ) )
525: 180 CONTINUE
526: END IF
527: END IF
528: *
529: IF( J2+KDN.GT.N ) THEN
530: *
531: * adjust J2 to keep within the bounds of the matrix
532: *
533: NR = NR - 1
534: J2 = J2 - KDN - 1
535: END IF
536: *
537: DO 190 J = J1, J2, KD1
538: *
539: * create nonzero element a(j+kd,j-1) outside the
540: * band and store it in WORK
541: *
542: WORK( J+KD ) = WORK( J )*AB( KD1, J )
543: AB( KD1, J ) = D( J )*AB( KD1, J )
544: 190 CONTINUE
545: 200 CONTINUE
546: 210 CONTINUE
547: END IF
548: *
549: IF( KD.GT.0 ) THEN
550: *
551: * make off-diagonal elements real and copy them to E
552: *
553: DO 220 I = 1, N - 1
554: T = AB( 2, I )
555: ABST = ABS( T )
556: AB( 2, I ) = ABST
557: E( I ) = ABST
558: IF( ABST.NE.ZERO ) THEN
559: T = T / ABST
560: ELSE
561: T = CONE
562: END IF
563: IF( I.LT.N-1 )
564: $ AB( 2, I+1 ) = AB( 2, I+1 )*T
565: IF( WANTQ ) THEN
566: CALL ZSCAL( N, T, Q( 1, I+1 ), 1 )
567: END IF
568: 220 CONTINUE
569: ELSE
570: *
571: * set E to zero if original matrix was diagonal
572: *
573: DO 230 I = 1, N - 1
574: E( I ) = ZERO
575: 230 CONTINUE
576: END IF
577: *
578: * copy diagonal elements to D
579: *
580: DO 240 I = 1, N
581: D( I ) = AB( 1, I )
582: 240 CONTINUE
583: END IF
584: *
585: RETURN
586: *
587: * End of ZHBTRD
588: *
589: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhbtrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack