![[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: *
150: *> \ingroup complex16OTHERcomputational
151: *
152: *> \par Further Details:
153: * =====================
154: *>
155: *> \verbatim
156: *>
157: *> Modified by Linda Kaufman, Bell Labs.
158: *> \endverbatim
159: *>
160: * =====================================================================
1.1 bertrand 161: SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
162: $ WORK, INFO )
163: *
1.17 ! bertrand 164: * -- LAPACK computational routine --
1.1 bertrand 165: * -- LAPACK is a software package provided by Univ. of Tennessee, --
166: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167: *
168: * .. Scalar Arguments ..
169: CHARACTER UPLO, VECT
170: INTEGER INFO, KD, LDAB, LDQ, N
171: * ..
172: * .. Array Arguments ..
173: DOUBLE PRECISION D( * ), E( * )
174: COMPLEX*16 AB( LDAB, * ), Q( LDQ, * ), WORK( * )
175: * ..
176: *
177: * =====================================================================
178: *
179: * .. Parameters ..
180: DOUBLE PRECISION ZERO
181: PARAMETER ( ZERO = 0.0D+0 )
182: COMPLEX*16 CZERO, CONE
183: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
184: $ CONE = ( 1.0D+0, 0.0D+0 ) )
185: * ..
186: * .. Local Scalars ..
187: LOGICAL INITQ, UPPER, WANTQ
188: INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
189: $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
190: $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
191: DOUBLE PRECISION ABST
192: COMPLEX*16 T, TEMP
193: * ..
194: * .. External Subroutines ..
195: EXTERNAL XERBLA, ZLACGV, ZLAR2V, ZLARGV, ZLARTG, ZLARTV,
196: $ ZLASET, ZROT, ZSCAL
197: * ..
198: * .. Intrinsic Functions ..
199: INTRINSIC ABS, DBLE, DCONJG, MAX, MIN
200: * ..
201: * .. External Functions ..
202: LOGICAL LSAME
203: EXTERNAL LSAME
204: * ..
205: * .. Executable Statements ..
206: *
207: * Test the input parameters
208: *
209: INITQ = LSAME( VECT, 'V' )
210: WANTQ = INITQ .OR. LSAME( VECT, 'U' )
211: UPPER = LSAME( UPLO, 'U' )
212: KD1 = KD + 1
213: KDM1 = KD - 1
214: INCX = LDAB - 1
215: IQEND = 1
216: *
217: INFO = 0
218: IF( .NOT.WANTQ .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( KD.LT.0 ) THEN
225: INFO = -4
226: ELSE IF( LDAB.LT.KD1 ) THEN
227: INFO = -6
228: ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
229: INFO = -10
230: END IF
231: IF( INFO.NE.0 ) THEN
232: CALL XERBLA( 'ZHBTRD', -INFO )
233: RETURN
234: END IF
235: *
236: * Quick return if possible
237: *
238: IF( N.EQ.0 )
239: $ RETURN
240: *
241: * Initialize Q to the unit matrix, if needed
242: *
243: IF( INITQ )
244: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
245: *
246: * Wherever possible, plane rotations are generated and applied in
247: * vector operations of length NR over the index set J1:J2:KD1.
248: *
249: * The real cosines and complex sines of the plane rotations are
250: * stored in the arrays D and WORK.
251: *
252: INCA = KD1*LDAB
253: KDN = MIN( N-1, KD )
254: IF( UPPER ) THEN
255: *
256: IF( KD.GT.1 ) THEN
257: *
258: * Reduce to complex Hermitian tridiagonal form, working with
259: * the upper triangle
260: *
261: NR = 0
262: J1 = KDN + 2
263: J2 = 1
264: *
265: AB( KD1, 1 ) = DBLE( AB( KD1, 1 ) )
266: DO 90 I = 1, N - 2
267: *
268: * Reduce i-th row of matrix to tridiagonal form
269: *
270: DO 80 K = KDN + 1, 2, -1
271: J1 = J1 + KDN
272: J2 = J2 + KDN
273: *
274: IF( NR.GT.0 ) THEN
275: *
276: * generate plane rotations to annihilate nonzero
277: * elements which have been created outside the band
278: *
279: CALL ZLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
280: $ KD1, D( J1 ), KD1 )
281: *
282: * apply rotations from the right
283: *
284: *
285: * Dependent on the the number of diagonals either
286: * ZLARTV or ZROT is used
287: *
288: IF( NR.GE.2*KD-1 ) THEN
289: DO 10 L = 1, KD - 1
290: CALL ZLARTV( NR, AB( L+1, J1-1 ), INCA,
291: $ AB( L, J1 ), INCA, D( J1 ),
292: $ WORK( J1 ), KD1 )
293: 10 CONTINUE
294: *
295: ELSE
296: JEND = J1 + ( NR-1 )*KD1
297: DO 20 JINC = J1, JEND, KD1
298: CALL ZROT( KDM1, AB( 2, JINC-1 ), 1,
299: $ AB( 1, JINC ), 1, D( JINC ),
300: $ WORK( JINC ) )
301: 20 CONTINUE
302: END IF
303: END IF
304: *
305: *
306: IF( K.GT.2 ) THEN
307: IF( K.LE.N-I+1 ) THEN
308: *
309: * generate plane rotation to annihilate a(i,i+k-1)
310: * within the band
311: *
312: CALL ZLARTG( AB( KD-K+3, I+K-2 ),
313: $ AB( KD-K+2, I+K-1 ), D( I+K-1 ),
314: $ WORK( I+K-1 ), TEMP )
315: AB( KD-K+3, I+K-2 ) = TEMP
316: *
317: * apply rotation from the right
318: *
319: CALL ZROT( K-3, AB( KD-K+4, I+K-2 ), 1,
320: $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
321: $ WORK( I+K-1 ) )
322: END IF
323: NR = NR + 1
324: J1 = J1 - KDN - 1
325: END IF
326: *
327: * apply plane rotations from both sides to diagonal
328: * blocks
329: *
330: IF( NR.GT.0 )
331: $ CALL ZLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
332: $ AB( KD, J1 ), INCA, D( J1 ),
333: $ WORK( J1 ), KD1 )
334: *
335: * apply plane rotations from the left
336: *
337: IF( NR.GT.0 ) THEN
338: CALL ZLACGV( NR, WORK( J1 ), KD1 )
339: IF( 2*KD-1.LT.NR ) THEN
340: *
341: * Dependent on the the number of diagonals either
342: * ZLARTV or ZROT is used
343: *
344: DO 30 L = 1, KD - 1
345: IF( J2+L.GT.N ) THEN
346: NRT = NR - 1
347: ELSE
348: NRT = NR
349: END IF
350: IF( NRT.GT.0 )
351: $ CALL ZLARTV( NRT, AB( KD-L, J1+L ), INCA,
352: $ AB( KD-L+1, J1+L ), INCA,
353: $ D( J1 ), WORK( J1 ), KD1 )
354: 30 CONTINUE
355: ELSE
356: J1END = J1 + KD1*( NR-2 )
357: IF( J1END.GE.J1 ) THEN
358: DO 40 JIN = J1, J1END, KD1
359: CALL ZROT( KD-1, AB( KD-1, JIN+1 ), INCX,
360: $ AB( KD, JIN+1 ), INCX,
361: $ D( JIN ), WORK( JIN ) )
362: 40 CONTINUE
363: END IF
364: LEND = MIN( KDM1, N-J2 )
365: LAST = J1END + KD1
366: IF( LEND.GT.0 )
367: $ CALL ZROT( LEND, AB( KD-1, LAST+1 ), INCX,
368: $ AB( KD, LAST+1 ), INCX, D( LAST ),
369: $ WORK( LAST ) )
370: END IF
371: END IF
372: *
373: IF( WANTQ ) THEN
374: *
375: * accumulate product of plane rotations in Q
376: *
377: IF( INITQ ) THEN
378: *
379: * take advantage of the fact that Q was
380: * initially the Identity matrix
381: *
382: IQEND = MAX( IQEND, J2 )
383: I2 = MAX( 0, K-3 )
384: IQAEND = 1 + I*KD
385: IF( K.EQ.2 )
386: $ IQAEND = IQAEND + KD
387: IQAEND = MIN( IQAEND, IQEND )
388: DO 50 J = J1, J2, KD1
389: IBL = I - I2 / KDM1
390: I2 = I2 + 1
391: IQB = MAX( 1, J-IBL )
392: NQ = 1 + IQAEND - IQB
393: IQAEND = MIN( IQAEND+KD, IQEND )
394: CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
395: $ 1, D( J ), DCONJG( WORK( J ) ) )
396: 50 CONTINUE
397: ELSE
398: *
399: DO 60 J = J1, J2, KD1
400: CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
401: $ D( J ), DCONJG( WORK( J ) ) )
402: 60 CONTINUE
403: END IF
404: *
405: END IF
406: *
407: IF( J2+KDN.GT.N ) THEN
408: *
409: * adjust J2 to keep within the bounds of the matrix
410: *
411: NR = NR - 1
412: J2 = J2 - KDN - 1
413: END IF
414: *
415: DO 70 J = J1, J2, KD1
416: *
417: * create nonzero element a(j-1,j+kd) outside the band
418: * and store it in WORK
419: *
420: WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
421: AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
422: 70 CONTINUE
423: 80 CONTINUE
424: 90 CONTINUE
425: END IF
426: *
427: IF( KD.GT.0 ) THEN
428: *
429: * make off-diagonal elements real and copy them to E
430: *
431: DO 100 I = 1, N - 1
432: T = AB( KD, I+1 )
433: ABST = ABS( T )
434: AB( KD, I+1 ) = ABST
435: E( I ) = ABST
436: IF( ABST.NE.ZERO ) THEN
437: T = T / ABST
438: ELSE
439: T = CONE
440: END IF
441: IF( I.LT.N-1 )
442: $ AB( KD, I+2 ) = AB( KD, I+2 )*T
443: IF( WANTQ ) THEN
444: CALL ZSCAL( N, DCONJG( T ), Q( 1, I+1 ), 1 )
445: END IF
446: 100 CONTINUE
447: ELSE
448: *
449: * set E to zero if original matrix was diagonal
450: *
451: DO 110 I = 1, N - 1
452: E( I ) = ZERO
453: 110 CONTINUE
454: END IF
455: *
456: * copy diagonal elements to D
457: *
458: DO 120 I = 1, N
1.17 ! bertrand 459: D( I ) = DBLE( AB( KD1, I ) )
1.1 bertrand 460: 120 CONTINUE
461: *
462: ELSE
463: *
464: IF( KD.GT.1 ) THEN
465: *
466: * Reduce to complex Hermitian tridiagonal form, working with
467: * the lower triangle
468: *
469: NR = 0
470: J1 = KDN + 2
471: J2 = 1
472: *
473: AB( 1, 1 ) = DBLE( AB( 1, 1 ) )
474: DO 210 I = 1, N - 2
475: *
476: * Reduce i-th column of matrix to tridiagonal form
477: *
478: DO 200 K = KDN + 1, 2, -1
479: J1 = J1 + KDN
480: J2 = J2 + KDN
481: *
482: IF( NR.GT.0 ) THEN
483: *
484: * generate plane rotations to annihilate nonzero
485: * elements which have been created outside the band
486: *
487: CALL ZLARGV( NR, AB( KD1, J1-KD1 ), INCA,
488: $ WORK( J1 ), KD1, D( J1 ), KD1 )
489: *
490: * apply plane rotations from one side
491: *
492: *
493: * Dependent on the the number of diagonals either
494: * ZLARTV or ZROT is used
495: *
496: IF( NR.GT.2*KD-1 ) THEN
497: DO 130 L = 1, KD - 1
498: CALL ZLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
499: $ AB( KD1-L+1, J1-KD1+L ), INCA,
500: $ D( J1 ), WORK( J1 ), KD1 )
501: 130 CONTINUE
502: ELSE
503: JEND = J1 + KD1*( NR-1 )
504: DO 140 JINC = J1, JEND, KD1
505: CALL ZROT( KDM1, AB( KD, JINC-KD ), INCX,
506: $ AB( KD1, JINC-KD ), INCX,
507: $ D( JINC ), WORK( JINC ) )
508: 140 CONTINUE
509: END IF
510: *
511: END IF
512: *
513: IF( K.GT.2 ) THEN
514: IF( K.LE.N-I+1 ) THEN
515: *
516: * generate plane rotation to annihilate a(i+k-1,i)
517: * within the band
518: *
519: CALL ZLARTG( AB( K-1, I ), AB( K, I ),
520: $ D( I+K-1 ), WORK( I+K-1 ), TEMP )
521: AB( K-1, I ) = TEMP
522: *
523: * apply rotation from the left
524: *
525: CALL ZROT( K-3, AB( K-2, I+1 ), LDAB-1,
526: $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
527: $ WORK( I+K-1 ) )
528: END IF
529: NR = NR + 1
530: J1 = J1 - KDN - 1
531: END IF
532: *
533: * apply plane rotations from both sides to diagonal
534: * blocks
535: *
536: IF( NR.GT.0 )
537: $ CALL ZLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
538: $ AB( 2, J1-1 ), INCA, D( J1 ),
539: $ WORK( J1 ), KD1 )
540: *
541: * apply plane rotations from the right
542: *
543: *
544: * Dependent on the the number of diagonals either
545: * ZLARTV or ZROT is used
546: *
547: IF( NR.GT.0 ) THEN
548: CALL ZLACGV( NR, WORK( J1 ), KD1 )
549: IF( NR.GT.2*KD-1 ) THEN
550: DO 150 L = 1, KD - 1
551: IF( J2+L.GT.N ) THEN
552: NRT = NR - 1
553: ELSE
554: NRT = NR
555: END IF
556: IF( NRT.GT.0 )
557: $ CALL ZLARTV( NRT, AB( L+2, J1-1 ), INCA,
558: $ AB( L+1, J1 ), INCA, D( J1 ),
559: $ WORK( J1 ), KD1 )
560: 150 CONTINUE
561: ELSE
562: J1END = J1 + KD1*( NR-2 )
563: IF( J1END.GE.J1 ) THEN
564: DO 160 J1INC = J1, J1END, KD1
565: CALL ZROT( KDM1, AB( 3, J1INC-1 ), 1,
566: $ AB( 2, J1INC ), 1, D( J1INC ),
567: $ WORK( J1INC ) )
568: 160 CONTINUE
569: END IF
570: LEND = MIN( KDM1, N-J2 )
571: LAST = J1END + KD1
572: IF( LEND.GT.0 )
573: $ CALL ZROT( LEND, AB( 3, LAST-1 ), 1,
574: $ AB( 2, LAST ), 1, D( LAST ),
575: $ WORK( LAST ) )
576: END IF
577: END IF
578: *
579: *
580: *
581: IF( WANTQ ) THEN
582: *
583: * accumulate product of plane rotations in Q
584: *
585: IF( INITQ ) THEN
586: *
587: * take advantage of the fact that Q was
588: * initially the Identity matrix
589: *
590: IQEND = MAX( IQEND, J2 )
591: I2 = MAX( 0, K-3 )
592: IQAEND = 1 + I*KD
593: IF( K.EQ.2 )
594: $ IQAEND = IQAEND + KD
595: IQAEND = MIN( IQAEND, IQEND )
596: DO 170 J = J1, J2, KD1
597: IBL = I - I2 / KDM1
598: I2 = I2 + 1
599: IQB = MAX( 1, J-IBL )
600: NQ = 1 + IQAEND - IQB
601: IQAEND = MIN( IQAEND+KD, IQEND )
602: CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
603: $ 1, D( J ), WORK( J ) )
604: 170 CONTINUE
605: ELSE
606: *
607: DO 180 J = J1, J2, KD1
608: CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
609: $ D( J ), WORK( J ) )
610: 180 CONTINUE
611: END IF
612: END IF
613: *
614: IF( J2+KDN.GT.N ) THEN
615: *
616: * adjust J2 to keep within the bounds of the matrix
617: *
618: NR = NR - 1
619: J2 = J2 - KDN - 1
620: END IF
621: *
622: DO 190 J = J1, J2, KD1
623: *
624: * create nonzero element a(j+kd,j-1) outside the
625: * band and store it in WORK
626: *
627: WORK( J+KD ) = WORK( J )*AB( KD1, J )
628: AB( KD1, J ) = D( J )*AB( KD1, J )
629: 190 CONTINUE
630: 200 CONTINUE
631: 210 CONTINUE
632: END IF
633: *
634: IF( KD.GT.0 ) THEN
635: *
636: * make off-diagonal elements real and copy them to E
637: *
638: DO 220 I = 1, N - 1
639: T = AB( 2, I )
640: ABST = ABS( T )
641: AB( 2, I ) = ABST
642: E( I ) = ABST
643: IF( ABST.NE.ZERO ) THEN
644: T = T / ABST
645: ELSE
646: T = CONE
647: END IF
648: IF( I.LT.N-1 )
649: $ AB( 2, I+1 ) = AB( 2, I+1 )*T
650: IF( WANTQ ) THEN
651: CALL ZSCAL( N, T, Q( 1, I+1 ), 1 )
652: END IF
653: 220 CONTINUE
654: ELSE
655: *
656: * set E to zero if original matrix was diagonal
657: *
658: DO 230 I = 1, N - 1
659: E( I ) = ZERO
660: 230 CONTINUE
661: END IF
662: *
663: * copy diagonal elements to D
664: *
665: DO 240 I = 1, N
1.17 ! bertrand 666: D( I ) = DBLE( AB( 1, I ) )
1.1 bertrand 667: 240 CONTINUE
668: END IF
669: *
670: RETURN
671: *
672: * End of ZHBTRD
673: *
674: END