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
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