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