001:       SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
002:      $                   WORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          UPLO, VECT
010:       INTEGER            INFO, KD, LDAB, LDQ, N
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   D( * ), E( * )
014:       COMPLEX*16         AB( LDAB, * ), Q( LDQ, * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  ZHBTRD reduces a complex Hermitian band matrix A to real symmetric
021: *  tridiagonal form T by a unitary similarity transformation:
022: *  Q**H * A * Q = T.
023: *
024: *  Arguments
025: *  =========
026: *
027: *  VECT    (input) CHARACTER*1
028: *          = 'N':  do not form Q;
029: *          = 'V':  form Q;
030: *          = 'U':  update a matrix X, by forming X*Q.
031: *
032: *  UPLO    (input) CHARACTER*1
033: *          = 'U':  Upper triangle of A is stored;
034: *          = 'L':  Lower triangle of A is stored.
035: *
036: *  N       (input) INTEGER
037: *          The order of the matrix A.  N >= 0.
038: *
039: *  KD      (input) INTEGER
040: *          The number of superdiagonals of the matrix A if UPLO = 'U',
041: *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
042: *
043: *  AB      (input/output) COMPLEX*16 array, dimension (LDAB,N)
044: *          On entry, the upper or lower triangle of the Hermitian band
045: *          matrix A, stored in the first KD+1 rows of the array.  The
046: *          j-th column of A is stored in the j-th column of the array AB
047: *          as follows:
048: *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
049: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
050: *          On exit, the diagonal elements of AB are overwritten by the
051: *          diagonal elements of the tridiagonal matrix T; if KD > 0, the
052: *          elements on the first superdiagonal (if UPLO = 'U') or the
053: *          first subdiagonal (if UPLO = 'L') are overwritten by the
054: *          off-diagonal elements of T; the rest of AB is overwritten by
055: *          values generated during the reduction.
056: *
057: *  LDAB    (input) INTEGER
058: *          The leading dimension of the array AB.  LDAB >= KD+1.
059: *
060: *  D       (output) DOUBLE PRECISION array, dimension (N)
061: *          The diagonal elements of the tridiagonal matrix T.
062: *
063: *  E       (output) DOUBLE PRECISION array, dimension (N-1)
064: *          The off-diagonal elements of the tridiagonal matrix T:
065: *          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
066: *
067: *  Q       (input/output) COMPLEX*16 array, dimension (LDQ,N)
068: *          On entry, if VECT = 'U', then Q must contain an N-by-N
069: *          matrix X; if VECT = 'N' or 'V', then Q need not be set.
070: *
071: *          On exit:
072: *          if VECT = 'V', Q contains the N-by-N unitary matrix Q;
073: *          if VECT = 'U', Q contains the product X*Q;
074: *          if VECT = 'N', the array Q is not referenced.
075: *
076: *  LDQ     (input) INTEGER
077: *          The leading dimension of the array Q.
078: *          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
079: *
080: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
081: *
082: *  INFO    (output) INTEGER
083: *          = 0:  successful exit
084: *          < 0:  if INFO = -i, the i-th argument had an illegal value
085: *
086: *  Further Details
087: *  ===============
088: *
089: *  Modified by Linda Kaufman, Bell Labs.
090: *
091: *  =====================================================================
092: *
093: *     .. Parameters ..
094:       DOUBLE PRECISION   ZERO
095:       PARAMETER          ( ZERO = 0.0D+0 )
096:       COMPLEX*16         CZERO, CONE
097:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
098:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
099: *     ..
100: *     .. Local Scalars ..
101:       LOGICAL            INITQ, UPPER, WANTQ
102:       INTEGER            I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
103:      $                   J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
104:      $                   KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
105:       DOUBLE PRECISION   ABST
106:       COMPLEX*16         T, TEMP
107: *     ..
108: *     .. External Subroutines ..
109:       EXTERNAL           XERBLA, ZLACGV, ZLAR2V, ZLARGV, ZLARTG, ZLARTV,
110:      $                   ZLASET, ZROT, ZSCAL
111: *     ..
112: *     .. Intrinsic Functions ..
113:       INTRINSIC          ABS, DBLE, DCONJG, MAX, MIN
114: *     ..
115: *     .. External Functions ..
116:       LOGICAL            LSAME
117:       EXTERNAL           LSAME
118: *     ..
119: *     .. Executable Statements ..
120: *
121: *     Test the input parameters
122: *
123:       INITQ = LSAME( VECT, 'V' )
124:       WANTQ = INITQ .OR. LSAME( VECT, 'U' )
125:       UPPER = LSAME( UPLO, 'U' )
126:       KD1 = KD + 1
127:       KDM1 = KD - 1
128:       INCX = LDAB - 1
129:       IQEND = 1
130: *
131:       INFO = 0
132:       IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
133:          INFO = -1
134:       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
135:          INFO = -2
136:       ELSE IF( N.LT.0 ) THEN
137:          INFO = -3
138:       ELSE IF( KD.LT.0 ) THEN
139:          INFO = -4
140:       ELSE IF( LDAB.LT.KD1 ) THEN
141:          INFO = -6
142:       ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
143:          INFO = -10
144:       END IF
145:       IF( INFO.NE.0 ) THEN
146:          CALL XERBLA( 'ZHBTRD', -INFO )
147:          RETURN
148:       END IF
149: *
150: *     Quick return if possible
151: *
152:       IF( N.EQ.0 )
153:      $   RETURN
154: *
155: *     Initialize Q to the unit matrix, if needed
156: *
157:       IF( INITQ )
158:      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
159: *
160: *     Wherever possible, plane rotations are generated and applied in
161: *     vector operations of length NR over the index set J1:J2:KD1.
162: *
163: *     The real cosines and complex sines of the plane rotations are
164: *     stored in the arrays D and WORK.
165: *
166:       INCA = KD1*LDAB
167:       KDN = MIN( N-1, KD )
168:       IF( UPPER ) THEN
169: *
170:          IF( KD.GT.1 ) THEN
171: *
172: *           Reduce to complex Hermitian tridiagonal form, working with
173: *           the upper triangle
174: *
175:             NR = 0
176:             J1 = KDN + 2
177:             J2 = 1
178: *
179:             AB( KD1, 1 ) = DBLE( AB( KD1, 1 ) )
180:             DO 90 I = 1, N - 2
181: *
182: *              Reduce i-th row of matrix to tridiagonal form
183: *
184:                DO 80 K = KDN + 1, 2, -1
185:                   J1 = J1 + KDN
186:                   J2 = J2 + KDN
187: *
188:                   IF( NR.GT.0 ) THEN
189: *
190: *                    generate plane rotations to annihilate nonzero
191: *                    elements which have been created outside the band
192: *
193:                      CALL ZLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
194:      $                            KD1, D( J1 ), KD1 )
195: *
196: *                    apply rotations from the right
197: *
198: *
199: *                    Dependent on the the number of diagonals either
200: *                    ZLARTV or ZROT is used
201: *
202:                      IF( NR.GE.2*KD-1 ) THEN
203:                         DO 10 L = 1, KD - 1
204:                            CALL ZLARTV( NR, AB( L+1, J1-1 ), INCA,
205:      $                                  AB( L, J1 ), INCA, D( J1 ),
206:      $                                  WORK( J1 ), KD1 )
207:    10                   CONTINUE
208: *
209:                      ELSE
210:                         JEND = J1 + ( NR-1 )*KD1
211:                         DO 20 JINC = J1, JEND, KD1
212:                            CALL ZROT( KDM1, AB( 2, JINC-1 ), 1,
213:      $                                AB( 1, JINC ), 1, D( JINC ),
214:      $                                WORK( JINC ) )
215:    20                   CONTINUE
216:                      END IF
217:                   END IF
218: *
219: *
220:                   IF( K.GT.2 ) THEN
221:                      IF( K.LE.N-I+1 ) THEN
222: *
223: *                       generate plane rotation to annihilate a(i,i+k-1)
224: *                       within the band
225: *
226:                         CALL ZLARTG( AB( KD-K+3, I+K-2 ),
227:      $                               AB( KD-K+2, I+K-1 ), D( I+K-1 ),
228:      $                               WORK( I+K-1 ), TEMP )
229:                         AB( KD-K+3, I+K-2 ) = TEMP
230: *
231: *                       apply rotation from the right
232: *
233:                         CALL ZROT( K-3, AB( KD-K+4, I+K-2 ), 1,
234:      $                             AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
235:      $                             WORK( I+K-1 ) )
236:                      END IF
237:                      NR = NR + 1
238:                      J1 = J1 - KDN - 1
239:                   END IF
240: *
241: *                 apply plane rotations from both sides to diagonal
242: *                 blocks
243: *
244:                   IF( NR.GT.0 )
245:      $               CALL ZLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
246:      $                            AB( KD, J1 ), INCA, D( J1 ),
247:      $                            WORK( J1 ), KD1 )
248: *
249: *                 apply plane rotations from the left
250: *
251:                   IF( NR.GT.0 ) THEN
252:                      CALL ZLACGV( NR, WORK( J1 ), KD1 )
253:                      IF( 2*KD-1.LT.NR ) THEN
254: *
255: *                    Dependent on the the number of diagonals either
256: *                    ZLARTV or ZROT is used
257: *
258:                         DO 30 L = 1, KD - 1
259:                            IF( J2+L.GT.N ) THEN
260:                               NRT = NR - 1
261:                            ELSE
262:                               NRT = NR
263:                            END IF
264:                            IF( NRT.GT.0 )
265:      $                        CALL ZLARTV( NRT, AB( KD-L, J1+L ), INCA,
266:      $                                     AB( KD-L+1, J1+L ), INCA,
267:      $                                     D( J1 ), WORK( J1 ), KD1 )
268:    30                   CONTINUE
269:                      ELSE
270:                         J1END = J1 + KD1*( NR-2 )
271:                         IF( J1END.GE.J1 ) THEN
272:                            DO 40 JIN = J1, J1END, KD1
273:                               CALL ZROT( KD-1, AB( KD-1, JIN+1 ), INCX,
274:      $                                   AB( KD, JIN+1 ), INCX,
275:      $                                   D( JIN ), WORK( JIN ) )
276:    40                      CONTINUE
277:                         END IF
278:                         LEND = MIN( KDM1, N-J2 )
279:                         LAST = J1END + KD1
280:                         IF( LEND.GT.0 )
281:      $                     CALL ZROT( LEND, AB( KD-1, LAST+1 ), INCX,
282:      $                                AB( KD, LAST+1 ), INCX, D( LAST ),
283:      $                                WORK( LAST ) )
284:                      END IF
285:                   END IF
286: *
287:                   IF( WANTQ ) THEN
288: *
289: *                    accumulate product of plane rotations in Q
290: *
291:                      IF( INITQ ) THEN
292: *
293: *                 take advantage of the fact that Q was
294: *                 initially the Identity matrix
295: *
296:                         IQEND = MAX( IQEND, J2 )
297:                         I2 = MAX( 0, K-3 )
298:                         IQAEND = 1 + I*KD
299:                         IF( K.EQ.2 )
300:      $                     IQAEND = IQAEND + KD
301:                         IQAEND = MIN( IQAEND, IQEND )
302:                         DO 50 J = J1, J2, KD1
303:                            IBL = I - I2 / KDM1
304:                            I2 = I2 + 1
305:                            IQB = MAX( 1, J-IBL )
306:                            NQ = 1 + IQAEND - IQB
307:                            IQAEND = MIN( IQAEND+KD, IQEND )
308:                            CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
309:      $                                1, D( J ), DCONJG( WORK( J ) ) )
310:    50                   CONTINUE
311:                      ELSE
312: *
313:                         DO 60 J = J1, J2, KD1
314:                            CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
315:      $                                D( J ), DCONJG( WORK( J ) ) )
316:    60                   CONTINUE
317:                      END IF
318: *
319:                   END IF
320: *
321:                   IF( J2+KDN.GT.N ) THEN
322: *
323: *                    adjust J2 to keep within the bounds of the matrix
324: *
325:                      NR = NR - 1
326:                      J2 = J2 - KDN - 1
327:                   END IF
328: *
329:                   DO 70 J = J1, J2, KD1
330: *
331: *                    create nonzero element a(j-1,j+kd) outside the band
332: *                    and store it in WORK
333: *
334:                      WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
335:                      AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
336:    70             CONTINUE
337:    80          CONTINUE
338:    90       CONTINUE
339:          END IF
340: *
341:          IF( KD.GT.0 ) THEN
342: *
343: *           make off-diagonal elements real and copy them to E
344: *
345:             DO 100 I = 1, N - 1
346:                T = AB( KD, I+1 )
347:                ABST = ABS( T )
348:                AB( KD, I+1 ) = ABST
349:                E( I ) = ABST
350:                IF( ABST.NE.ZERO ) THEN
351:                   T = T / ABST
352:                ELSE
353:                   T = CONE
354:                END IF
355:                IF( I.LT.N-1 )
356:      $            AB( KD, I+2 ) = AB( KD, I+2 )*T
357:                IF( WANTQ ) THEN
358:                   CALL ZSCAL( N, DCONJG( T ), Q( 1, I+1 ), 1 )
359:                END IF
360:   100       CONTINUE
361:          ELSE
362: *
363: *           set E to zero if original matrix was diagonal
364: *
365:             DO 110 I = 1, N - 1
366:                E( I ) = ZERO
367:   110       CONTINUE
368:          END IF
369: *
370: *        copy diagonal elements to D
371: *
372:          DO 120 I = 1, N
373:             D( I ) = AB( KD1, I )
374:   120    CONTINUE
375: *
376:       ELSE
377: *
378:          IF( KD.GT.1 ) THEN
379: *
380: *           Reduce to complex Hermitian tridiagonal form, working with
381: *           the lower triangle
382: *
383:             NR = 0
384:             J1 = KDN + 2
385:             J2 = 1
386: *
387:             AB( 1, 1 ) = DBLE( AB( 1, 1 ) )
388:             DO 210 I = 1, N - 2
389: *
390: *              Reduce i-th column of matrix to tridiagonal form
391: *
392:                DO 200 K = KDN + 1, 2, -1
393:                   J1 = J1 + KDN
394:                   J2 = J2 + KDN
395: *
396:                   IF( NR.GT.0 ) THEN
397: *
398: *                    generate plane rotations to annihilate nonzero
399: *                    elements which have been created outside the band
400: *
401:                      CALL ZLARGV( NR, AB( KD1, J1-KD1 ), INCA,
402:      $                            WORK( J1 ), KD1, D( J1 ), KD1 )
403: *
404: *                    apply plane rotations from one side
405: *
406: *
407: *                    Dependent on the the number of diagonals either
408: *                    ZLARTV or ZROT is used
409: *
410:                      IF( NR.GT.2*KD-1 ) THEN
411:                         DO 130 L = 1, KD - 1
412:                            CALL ZLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
413:      $                                  AB( KD1-L+1, J1-KD1+L ), INCA,
414:      $                                  D( J1 ), WORK( J1 ), KD1 )
415:   130                   CONTINUE
416:                      ELSE
417:                         JEND = J1 + KD1*( NR-1 )
418:                         DO 140 JINC = J1, JEND, KD1
419:                            CALL ZROT( KDM1, AB( KD, JINC-KD ), INCX,
420:      $                                AB( KD1, JINC-KD ), INCX,
421:      $                                D( JINC ), WORK( JINC ) )
422:   140                   CONTINUE
423:                      END IF
424: *
425:                   END IF
426: *
427:                   IF( K.GT.2 ) THEN
428:                      IF( K.LE.N-I+1 ) THEN
429: *
430: *                       generate plane rotation to annihilate a(i+k-1,i)
431: *                       within the band
432: *
433:                         CALL ZLARTG( AB( K-1, I ), AB( K, I ),
434:      $                               D( I+K-1 ), WORK( I+K-1 ), TEMP )
435:                         AB( K-1, I ) = TEMP
436: *
437: *                       apply rotation from the left
438: *
439:                         CALL ZROT( K-3, AB( K-2, I+1 ), LDAB-1,
440:      $                             AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
441:      $                             WORK( I+K-1 ) )
442:                      END IF
443:                      NR = NR + 1
444:                      J1 = J1 - KDN - 1
445:                   END IF
446: *
447: *                 apply plane rotations from both sides to diagonal
448: *                 blocks
449: *
450:                   IF( NR.GT.0 )
451:      $               CALL ZLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
452:      $                            AB( 2, J1-1 ), INCA, D( J1 ),
453:      $                            WORK( J1 ), KD1 )
454: *
455: *                 apply plane rotations from the right
456: *
457: *
458: *                    Dependent on the the number of diagonals either
459: *                    ZLARTV or ZROT is used
460: *
461:                   IF( NR.GT.0 ) THEN
462:                      CALL ZLACGV( NR, WORK( J1 ), KD1 )
463:                      IF( NR.GT.2*KD-1 ) THEN
464:                         DO 150 L = 1, KD - 1
465:                            IF( J2+L.GT.N ) THEN
466:                               NRT = NR - 1
467:                            ELSE
468:                               NRT = NR
469:                            END IF
470:                            IF( NRT.GT.0 )
471:      $                        CALL ZLARTV( NRT, AB( L+2, J1-1 ), INCA,
472:      $                                     AB( L+1, J1 ), INCA, D( J1 ),
473:      $                                     WORK( J1 ), KD1 )
474:   150                   CONTINUE
475:                      ELSE
476:                         J1END = J1 + KD1*( NR-2 )
477:                         IF( J1END.GE.J1 ) THEN
478:                            DO 160 J1INC = J1, J1END, KD1
479:                               CALL ZROT( KDM1, AB( 3, J1INC-1 ), 1,
480:      $                                   AB( 2, J1INC ), 1, D( J1INC ),
481:      $                                   WORK( J1INC ) )
482:   160                      CONTINUE
483:                         END IF
484:                         LEND = MIN( KDM1, N-J2 )
485:                         LAST = J1END + KD1
486:                         IF( LEND.GT.0 )
487:      $                     CALL ZROT( LEND, AB( 3, LAST-1 ), 1,
488:      $                                AB( 2, LAST ), 1, D( LAST ),
489:      $                                WORK( LAST ) )
490:                      END IF
491:                   END IF
492: *
493: *
494: *
495:                   IF( WANTQ ) THEN
496: *
497: *                    accumulate product of plane rotations in Q
498: *
499:                      IF( INITQ ) THEN
500: *
501: *                 take advantage of the fact that Q was
502: *                 initially the Identity matrix
503: *
504:                         IQEND = MAX( IQEND, J2 )
505:                         I2 = MAX( 0, K-3 )
506:                         IQAEND = 1 + I*KD
507:                         IF( K.EQ.2 )
508:      $                     IQAEND = IQAEND + KD
509:                         IQAEND = MIN( IQAEND, IQEND )
510:                         DO 170 J = J1, J2, KD1
511:                            IBL = I - I2 / KDM1
512:                            I2 = I2 + 1
513:                            IQB = MAX( 1, J-IBL )
514:                            NQ = 1 + IQAEND - IQB
515:                            IQAEND = MIN( IQAEND+KD, IQEND )
516:                            CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
517:      $                                1, D( J ), WORK( J ) )
518:   170                   CONTINUE
519:                      ELSE
520: *
521:                         DO 180 J = J1, J2, KD1
522:                            CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
523:      $                                D( J ), WORK( J ) )
524:   180                   CONTINUE
525:                      END IF
526:                   END IF
527: *
528:                   IF( J2+KDN.GT.N ) THEN
529: *
530: *                    adjust J2 to keep within the bounds of the matrix
531: *
532:                      NR = NR - 1
533:                      J2 = J2 - KDN - 1
534:                   END IF
535: *
536:                   DO 190 J = J1, J2, KD1
537: *
538: *                    create nonzero element a(j+kd,j-1) outside the
539: *                    band and store it in WORK
540: *
541:                      WORK( J+KD ) = WORK( J )*AB( KD1, J )
542:                      AB( KD1, J ) = D( J )*AB( KD1, J )
543:   190             CONTINUE
544:   200          CONTINUE
545:   210       CONTINUE
546:          END IF
547: *
548:          IF( KD.GT.0 ) THEN
549: *
550: *           make off-diagonal elements real and copy them to E
551: *
552:             DO 220 I = 1, N - 1
553:                T = AB( 2, I )
554:                ABST = ABS( T )
555:                AB( 2, I ) = ABST
556:                E( I ) = ABST
557:                IF( ABST.NE.ZERO ) THEN
558:                   T = T / ABST
559:                ELSE
560:                   T = CONE
561:                END IF
562:                IF( I.LT.N-1 )
563:      $            AB( 2, I+1 ) = AB( 2, I+1 )*T
564:                IF( WANTQ ) THEN
565:                   CALL ZSCAL( N, T, Q( 1, I+1 ), 1 )
566:                END IF
567:   220       CONTINUE
568:          ELSE
569: *
570: *           set E to zero if original matrix was diagonal
571: *
572:             DO 230 I = 1, N - 1
573:                E( I ) = ZERO
574:   230       CONTINUE
575:          END IF
576: *
577: *        copy diagonal elements to D
578: *
579:          DO 240 I = 1, N
580:             D( I ) = AB( 1, I )
581:   240    CONTINUE
582:       END IF
583: *
584:       RETURN
585: *
586: *     End of ZHBTRD
587: *
588:       END
589: