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