001:       SUBROUTINE DSBTRD( 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   AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
014:      $                   WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DSBTRD reduces a real symmetric band matrix A to symmetric
021: *  tridiagonal form T by an orthogonal similarity transformation:
022: *  Q**T * 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) DOUBLE PRECISION array, dimension (LDAB,N)
044: *          On entry, the upper or lower triangle of the symmetric 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) DOUBLE PRECISION 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 orthogonal 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) DOUBLE PRECISION 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, ONE
095:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
096: *     ..
097: *     .. Local Scalars ..
098:       LOGICAL            INITQ, UPPER, WANTQ
099:       INTEGER            I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
100:      $                   J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
101:      $                   KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
102:       DOUBLE PRECISION   TEMP
103: *     ..
104: *     .. External Subroutines ..
105:       EXTERNAL           DLAR2V, DLARGV, DLARTG, DLARTV, DLASET, DROT,
106:      $                   XERBLA
107: *     ..
108: *     .. Intrinsic Functions ..
109:       INTRINSIC          MAX, MIN
110: *     ..
111: *     .. External Functions ..
112:       LOGICAL            LSAME
113:       EXTERNAL           LSAME
114: *     ..
115: *     .. Executable Statements ..
116: *
117: *     Test the input parameters
118: *
119:       INITQ = LSAME( VECT, 'V' )
120:       WANTQ = INITQ .OR. LSAME( VECT, 'U' )
121:       UPPER = LSAME( UPLO, 'U' )
122:       KD1 = KD + 1
123:       KDM1 = KD - 1
124:       INCX = LDAB - 1
125:       IQEND = 1
126: *
127:       INFO = 0
128:       IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
129:          INFO = -1
130:       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
131:          INFO = -2
132:       ELSE IF( N.LT.0 ) THEN
133:          INFO = -3
134:       ELSE IF( KD.LT.0 ) THEN
135:          INFO = -4
136:       ELSE IF( LDAB.LT.KD1 ) THEN
137:          INFO = -6
138:       ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
139:          INFO = -10
140:       END IF
141:       IF( INFO.NE.0 ) THEN
142:          CALL XERBLA( 'DSBTRD', -INFO )
143:          RETURN
144:       END IF
145: *
146: *     Quick return if possible
147: *
148:       IF( N.EQ.0 )
149:      $   RETURN
150: *
151: *     Initialize Q to the unit matrix, if needed
152: *
153:       IF( INITQ )
154:      $   CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
155: *
156: *     Wherever possible, plane rotations are generated and applied in
157: *     vector operations of length NR over the index set J1:J2:KD1.
158: *
159: *     The cosines and sines of the plane rotations are stored in the
160: *     arrays D and WORK.
161: *
162:       INCA = KD1*LDAB
163:       KDN = MIN( N-1, KD )
164:       IF( UPPER ) THEN
165: *
166:          IF( KD.GT.1 ) THEN
167: *
168: *           Reduce to tridiagonal form, working with upper triangle
169: *
170:             NR = 0
171:             J1 = KDN + 2
172:             J2 = 1
173: *
174:             DO 90 I = 1, N - 2
175: *
176: *              Reduce i-th row of matrix to tridiagonal form
177: *
178:                DO 80 K = KDN + 1, 2, -1
179:                   J1 = J1 + KDN
180:                   J2 = J2 + KDN
181: *
182:                   IF( NR.GT.0 ) THEN
183: *
184: *                    generate plane rotations to annihilate nonzero
185: *                    elements which have been created outside the band
186: *
187:                      CALL DLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
188:      $                            KD1, D( J1 ), KD1 )
189: *
190: *                    apply rotations from the right
191: *
192: *
193: *                    Dependent on the the number of diagonals either
194: *                    DLARTV or DROT is used
195: *
196:                      IF( NR.GE.2*KD-1 ) THEN
197:                         DO 10 L = 1, KD - 1
198:                            CALL DLARTV( NR, AB( L+1, J1-1 ), INCA,
199:      $                                  AB( L, J1 ), INCA, D( J1 ),
200:      $                                  WORK( J1 ), KD1 )
201:    10                   CONTINUE
202: *
203:                      ELSE
204:                         JEND = J1 + ( NR-1 )*KD1
205:                         DO 20 JINC = J1, JEND, KD1
206:                            CALL DROT( KDM1, AB( 2, JINC-1 ), 1,
207:      $                                AB( 1, JINC ), 1, D( JINC ),
208:      $                                WORK( JINC ) )
209:    20                   CONTINUE
210:                      END IF
211:                   END IF
212: *
213: *
214:                   IF( K.GT.2 ) THEN
215:                      IF( K.LE.N-I+1 ) THEN
216: *
217: *                       generate plane rotation to annihilate a(i,i+k-1)
218: *                       within the band
219: *
220:                         CALL DLARTG( AB( KD-K+3, I+K-2 ),
221:      $                               AB( KD-K+2, I+K-1 ), D( I+K-1 ),
222:      $                               WORK( I+K-1 ), TEMP )
223:                         AB( KD-K+3, I+K-2 ) = TEMP
224: *
225: *                       apply rotation from the right
226: *
227:                         CALL DROT( K-3, AB( KD-K+4, I+K-2 ), 1,
228:      $                             AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
229:      $                             WORK( I+K-1 ) )
230:                      END IF
231:                      NR = NR + 1
232:                      J1 = J1 - KDN - 1
233:                   END IF
234: *
235: *                 apply plane rotations from both sides to diagonal
236: *                 blocks
237: *
238:                   IF( NR.GT.0 )
239:      $               CALL DLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
240:      $                            AB( KD, J1 ), INCA, D( J1 ),
241:      $                            WORK( J1 ), KD1 )
242: *
243: *                 apply plane rotations from the left
244: *
245:                   IF( NR.GT.0 ) THEN
246:                      IF( 2*KD-1.LT.NR ) THEN
247: *
248: *                    Dependent on the the number of diagonals either
249: *                    DLARTV or DROT is used
250: *
251:                         DO 30 L = 1, KD - 1
252:                            IF( J2+L.GT.N ) THEN
253:                               NRT = NR - 1
254:                            ELSE
255:                               NRT = NR
256:                            END IF
257:                            IF( NRT.GT.0 )
258:      $                        CALL DLARTV( NRT, AB( KD-L, J1+L ), INCA,
259:      $                                     AB( KD-L+1, J1+L ), INCA,
260:      $                                     D( J1 ), WORK( J1 ), KD1 )
261:    30                   CONTINUE
262:                      ELSE
263:                         J1END = J1 + KD1*( NR-2 )
264:                         IF( J1END.GE.J1 ) THEN
265:                            DO 40 JIN = J1, J1END, KD1
266:                               CALL DROT( KD-1, AB( KD-1, JIN+1 ), INCX,
267:      $                                   AB( KD, JIN+1 ), INCX,
268:      $                                   D( JIN ), WORK( JIN ) )
269:    40                      CONTINUE
270:                         END IF
271:                         LEND = MIN( KDM1, N-J2 )
272:                         LAST = J1END + KD1
273:                         IF( LEND.GT.0 )
274:      $                     CALL DROT( LEND, AB( KD-1, LAST+1 ), INCX,
275:      $                                AB( KD, LAST+1 ), INCX, D( LAST ),
276:      $                                WORK( LAST ) )
277:                      END IF
278:                   END IF
279: *
280:                   IF( WANTQ ) THEN
281: *
282: *                    accumulate product of plane rotations in Q
283: *
284:                      IF( INITQ ) THEN
285: *
286: *                 take advantage of the fact that Q was
287: *                 initially the Identity matrix
288: *
289:                         IQEND = MAX( IQEND, J2 )
290:                         I2 = MAX( 0, K-3 )
291:                         IQAEND = 1 + I*KD
292:                         IF( K.EQ.2 )
293:      $                     IQAEND = IQAEND + KD
294:                         IQAEND = MIN( IQAEND, IQEND )
295:                         DO 50 J = J1, J2, KD1
296:                            IBL = I - I2 / KDM1
297:                            I2 = I2 + 1
298:                            IQB = MAX( 1, J-IBL )
299:                            NQ = 1 + IQAEND - IQB
300:                            IQAEND = MIN( IQAEND+KD, IQEND )
301:                            CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
302:      $                                1, D( J ), WORK( J ) )
303:    50                   CONTINUE
304:                      ELSE
305: *
306:                         DO 60 J = J1, J2, KD1
307:                            CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
308:      $                                D( J ), WORK( J ) )
309:    60                   CONTINUE
310:                      END IF
311: *
312:                   END IF
313: *
314:                   IF( J2+KDN.GT.N ) THEN
315: *
316: *                    adjust J2 to keep within the bounds of the matrix
317: *
318:                      NR = NR - 1
319:                      J2 = J2 - KDN - 1
320:                   END IF
321: *
322:                   DO 70 J = J1, J2, KD1
323: *
324: *                    create nonzero element a(j-1,j+kd) outside the band
325: *                    and store it in WORK
326: *
327:                      WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
328:                      AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
329:    70             CONTINUE
330:    80          CONTINUE
331:    90       CONTINUE
332:          END IF
333: *
334:          IF( KD.GT.0 ) THEN
335: *
336: *           copy off-diagonal elements to E
337: *
338:             DO 100 I = 1, N - 1
339:                E( I ) = AB( KD, I+1 )
340:   100       CONTINUE
341:          ELSE
342: *
343: *           set E to zero if original matrix was diagonal
344: *
345:             DO 110 I = 1, N - 1
346:                E( I ) = ZERO
347:   110       CONTINUE
348:          END IF
349: *
350: *        copy diagonal elements to D
351: *
352:          DO 120 I = 1, N
353:             D( I ) = AB( KD1, I )
354:   120    CONTINUE
355: *
356:       ELSE
357: *
358:          IF( KD.GT.1 ) THEN
359: *
360: *           Reduce to tridiagonal form, working with lower triangle
361: *
362:             NR = 0
363:             J1 = KDN + 2
364:             J2 = 1
365: *
366:             DO 210 I = 1, N - 2
367: *
368: *              Reduce i-th column of matrix to tridiagonal form
369: *
370:                DO 200 K = KDN + 1, 2, -1
371:                   J1 = J1 + KDN
372:                   J2 = J2 + KDN
373: *
374:                   IF( NR.GT.0 ) THEN
375: *
376: *                    generate plane rotations to annihilate nonzero
377: *                    elements which have been created outside the band
378: *
379:                      CALL DLARGV( NR, AB( KD1, J1-KD1 ), INCA,
380:      $                            WORK( J1 ), KD1, D( J1 ), KD1 )
381: *
382: *                    apply plane rotations from one side
383: *
384: *
385: *                    Dependent on the the number of diagonals either
386: *                    DLARTV or DROT is used
387: *
388:                      IF( NR.GT.2*KD-1 ) THEN
389:                         DO 130 L = 1, KD - 1
390:                            CALL DLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
391:      $                                  AB( KD1-L+1, J1-KD1+L ), INCA,
392:      $                                  D( J1 ), WORK( J1 ), KD1 )
393:   130                   CONTINUE
394:                      ELSE
395:                         JEND = J1 + KD1*( NR-1 )
396:                         DO 140 JINC = J1, JEND, KD1
397:                            CALL DROT( KDM1, AB( KD, JINC-KD ), INCX,
398:      $                                AB( KD1, JINC-KD ), INCX,
399:      $                                D( JINC ), WORK( JINC ) )
400:   140                   CONTINUE
401:                      END IF
402: *
403:                   END IF
404: *
405:                   IF( K.GT.2 ) THEN
406:                      IF( K.LE.N-I+1 ) THEN
407: *
408: *                       generate plane rotation to annihilate a(i+k-1,i)
409: *                       within the band
410: *
411:                         CALL DLARTG( AB( K-1, I ), AB( K, I ),
412:      $                               D( I+K-1 ), WORK( I+K-1 ), TEMP )
413:                         AB( K-1, I ) = TEMP
414: *
415: *                       apply rotation from the left
416: *
417:                         CALL DROT( K-3, AB( K-2, I+1 ), LDAB-1,
418:      $                             AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
419:      $                             WORK( I+K-1 ) )
420:                      END IF
421:                      NR = NR + 1
422:                      J1 = J1 - KDN - 1
423:                   END IF
424: *
425: *                 apply plane rotations from both sides to diagonal
426: *                 blocks
427: *
428:                   IF( NR.GT.0 )
429:      $               CALL DLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
430:      $                            AB( 2, J1-1 ), INCA, D( J1 ),
431:      $                            WORK( J1 ), KD1 )
432: *
433: *                 apply plane rotations from the right
434: *
435: *
436: *                    Dependent on the the number of diagonals either
437: *                    DLARTV or DROT is used
438: *
439:                   IF( NR.GT.0 ) THEN
440:                      IF( NR.GT.2*KD-1 ) THEN
441:                         DO 150 L = 1, KD - 1
442:                            IF( J2+L.GT.N ) THEN
443:                               NRT = NR - 1
444:                            ELSE
445:                               NRT = NR
446:                            END IF
447:                            IF( NRT.GT.0 )
448:      $                        CALL DLARTV( NRT, AB( L+2, J1-1 ), INCA,
449:      $                                     AB( L+1, J1 ), INCA, D( J1 ),
450:      $                                     WORK( J1 ), KD1 )
451:   150                   CONTINUE
452:                      ELSE
453:                         J1END = J1 + KD1*( NR-2 )
454:                         IF( J1END.GE.J1 ) THEN
455:                            DO 160 J1INC = J1, J1END, KD1
456:                               CALL DROT( KDM1, AB( 3, J1INC-1 ), 1,
457:      $                                   AB( 2, J1INC ), 1, D( J1INC ),
458:      $                                   WORK( J1INC ) )
459:   160                      CONTINUE
460:                         END IF
461:                         LEND = MIN( KDM1, N-J2 )
462:                         LAST = J1END + KD1
463:                         IF( LEND.GT.0 )
464:      $                     CALL DROT( LEND, AB( 3, LAST-1 ), 1,
465:      $                                AB( 2, LAST ), 1, D( LAST ),
466:      $                                WORK( LAST ) )
467:                      END IF
468:                   END IF
469: *
470: *
471: *
472:                   IF( WANTQ ) THEN
473: *
474: *                    accumulate product of plane rotations in Q
475: *
476:                      IF( INITQ ) THEN
477: *
478: *                 take advantage of the fact that Q was
479: *                 initially the Identity matrix
480: *
481:                         IQEND = MAX( IQEND, J2 )
482:                         I2 = MAX( 0, K-3 )
483:                         IQAEND = 1 + I*KD
484:                         IF( K.EQ.2 )
485:      $                     IQAEND = IQAEND + KD
486:                         IQAEND = MIN( IQAEND, IQEND )
487:                         DO 170 J = J1, J2, KD1
488:                            IBL = I - I2 / KDM1
489:                            I2 = I2 + 1
490:                            IQB = MAX( 1, J-IBL )
491:                            NQ = 1 + IQAEND - IQB
492:                            IQAEND = MIN( IQAEND+KD, IQEND )
493:                            CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
494:      $                                1, D( J ), WORK( J ) )
495:   170                   CONTINUE
496:                      ELSE
497: *
498:                         DO 180 J = J1, J2, KD1
499:                            CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
500:      $                                D( J ), WORK( J ) )
501:   180                   CONTINUE
502:                      END IF
503:                   END IF
504: *
505:                   IF( J2+KDN.GT.N ) THEN
506: *
507: *                    adjust J2 to keep within the bounds of the matrix
508: *
509:                      NR = NR - 1
510:                      J2 = J2 - KDN - 1
511:                   END IF
512: *
513:                   DO 190 J = J1, J2, KD1
514: *
515: *                    create nonzero element a(j+kd,j-1) outside the
516: *                    band and store it in WORK
517: *
518:                      WORK( J+KD ) = WORK( J )*AB( KD1, J )
519:                      AB( KD1, J ) = D( J )*AB( KD1, J )
520:   190             CONTINUE
521:   200          CONTINUE
522:   210       CONTINUE
523:          END IF
524: *
525:          IF( KD.GT.0 ) THEN
526: *
527: *           copy off-diagonal elements to E
528: *
529:             DO 220 I = 1, N - 1
530:                E( I ) = AB( 2, I )
531:   220       CONTINUE
532:          ELSE
533: *
534: *           set E to zero if original matrix was diagonal
535: *
536:             DO 230 I = 1, N - 1
537:                E( I ) = ZERO
538:   230       CONTINUE
539:          END IF
540: *
541: *        copy diagonal elements to D
542: *
543:          DO 240 I = 1, N
544:             D( I ) = AB( 1, I )
545:   240    CONTINUE
546:       END IF
547: *
548:       RETURN
549: *
550: *     End of DSBTRD
551: *
552:       END
553: