001:       SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
002:      $                   WORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          UPLO, VECT
011:       INTEGER            INFO, KD, LDAB, LDQ, N
012: *     ..
013: *     .. Array Arguments ..
014:       DOUBLE PRECISION   D( * ), E( * )
015:       COMPLEX*16         AB( LDAB, * ), Q( LDQ, * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  ZHBTRD reduces a complex Hermitian band matrix A to real symmetric
022: *  tridiagonal form T by a unitary similarity transformation:
023: *  Q**H * A * Q = T.
024: *
025: *  Arguments
026: *  =========
027: *
028: *  VECT    (input) CHARACTER*1
029: *          = 'N':  do not form Q;
030: *          = 'V':  form Q;
031: *          = 'U':  update a matrix X, by forming X*Q.
032: *
033: *  UPLO    (input) CHARACTER*1
034: *          = 'U':  Upper triangle of A is stored;
035: *          = 'L':  Lower triangle of A is stored.
036: *
037: *  N       (input) INTEGER
038: *          The order of the matrix A.  N >= 0.
039: *
040: *  KD      (input) INTEGER
041: *          The number of superdiagonals of the matrix A if UPLO = 'U',
042: *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
043: *
044: *  AB      (input/output) COMPLEX*16 array, dimension (LDAB,N)
045: *          On entry, the upper or lower triangle of the Hermitian band
046: *          matrix A, stored in the first KD+1 rows of the array.  The
047: *          j-th column of A is stored in the j-th column of the array AB
048: *          as follows:
049: *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
050: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
051: *          On exit, the diagonal elements of AB are overwritten by the
052: *          diagonal elements of the tridiagonal matrix T; if KD > 0, the
053: *          elements on the first superdiagonal (if UPLO = 'U') or the
054: *          first subdiagonal (if UPLO = 'L') are overwritten by the
055: *          off-diagonal elements of T; the rest of AB is overwritten by
056: *          values generated during the reduction.
057: *
058: *  LDAB    (input) INTEGER
059: *          The leading dimension of the array AB.  LDAB >= KD+1.
060: *
061: *  D       (output) DOUBLE PRECISION array, dimension (N)
062: *          The diagonal elements of the tridiagonal matrix T.
063: *
064: *  E       (output) DOUBLE PRECISION array, dimension (N-1)
065: *          The off-diagonal elements of the tridiagonal matrix T:
066: *          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
067: *
068: *  Q       (input/output) COMPLEX*16 array, dimension (LDQ,N)
069: *          On entry, if VECT = 'U', then Q must contain an N-by-N
070: *          matrix X; if VECT = 'N' or 'V', then Q need not be set.
071: *
072: *          On exit:
073: *          if VECT = 'V', Q contains the N-by-N unitary matrix Q;
074: *          if VECT = 'U', Q contains the product X*Q;
075: *          if VECT = 'N', the array Q is not referenced.
076: *
077: *  LDQ     (input) INTEGER
078: *          The leading dimension of the array Q.
079: *          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
080: *
081: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
082: *
083: *  INFO    (output) INTEGER
084: *          = 0:  successful exit
085: *          < 0:  if INFO = -i, the i-th argument had an illegal value
086: *
087: *  Further Details
088: *  ===============
089: *
090: *  Modified by Linda Kaufman, Bell Labs.
091: *
092: *  =====================================================================
093: *
094: *     .. Parameters ..
095:       DOUBLE PRECISION   ZERO
096:       PARAMETER          ( ZERO = 0.0D+0 )
097:       COMPLEX*16         CZERO, CONE
098:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
099:      $                   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
590: