1: SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
2: $ WORK, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER UPLO, VECT
11: INTEGER INFO, KD, LDAB, LDQ, N
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
15: $ WORK( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DSBTRD reduces a real symmetric band matrix A to symmetric
22: * tridiagonal form T by an orthogonal similarity transformation:
23: * Q**T * A * Q = T.
24: *
25: * Arguments
26: * =========
27: *
28: * VECT (input) CHARACTER*1
29: * = 'N': do not form Q;
30: * = 'V': form Q;
31: * = 'U': update a matrix X, by forming X*Q.
32: *
33: * UPLO (input) CHARACTER*1
34: * = 'U': Upper triangle of A is stored;
35: * = 'L': Lower triangle of A is stored.
36: *
37: * N (input) INTEGER
38: * The order of the matrix A. N >= 0.
39: *
40: * KD (input) INTEGER
41: * The number of superdiagonals of the matrix A if UPLO = 'U',
42: * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
43: *
44: * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
45: * On entry, the upper or lower triangle of the symmetric band
46: * matrix A, stored in the first KD+1 rows of the array. The
47: * j-th column of A is stored in the j-th column of the array AB
48: * as follows:
49: * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
50: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
51: * On exit, the diagonal elements of AB are overwritten by the
52: * diagonal elements of the tridiagonal matrix T; if KD > 0, the
53: * elements on the first superdiagonal (if UPLO = 'U') or the
54: * first subdiagonal (if UPLO = 'L') are overwritten by the
55: * off-diagonal elements of T; the rest of AB is overwritten by
56: * values generated during the reduction.
57: *
58: * LDAB (input) INTEGER
59: * The leading dimension of the array AB. LDAB >= KD+1.
60: *
61: * D (output) DOUBLE PRECISION array, dimension (N)
62: * The diagonal elements of the tridiagonal matrix T.
63: *
64: * E (output) DOUBLE PRECISION array, dimension (N-1)
65: * The off-diagonal elements of the tridiagonal matrix T:
66: * E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
67: *
68: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
69: * On entry, if VECT = 'U', then Q must contain an N-by-N
70: * matrix X; if VECT = 'N' or 'V', then Q need not be set.
71: *
72: * On exit:
73: * if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
74: * if VECT = 'U', Q contains the product X*Q;
75: * if VECT = 'N', the array Q is not referenced.
76: *
77: * LDQ (input) INTEGER
78: * The leading dimension of the array Q.
79: * LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
80: *
81: * WORK (workspace) DOUBLE PRECISION array, dimension (N)
82: *
83: * INFO (output) INTEGER
84: * = 0: successful exit
85: * < 0: if INFO = -i, the i-th argument had an illegal value
86: *
87: * Further Details
88: * ===============
89: *
90: * Modified by Linda Kaufman, Bell Labs.
91: *
92: * =====================================================================
93: *
94: * .. Parameters ..
95: DOUBLE PRECISION ZERO, ONE
96: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
97: * ..
98: * .. Local Scalars ..
99: 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: DOUBLE PRECISION TEMP
104: * ..
105: * .. External Subroutines ..
106: EXTERNAL DLAR2V, DLARGV, DLARTG, DLARTV, DLASET, DROT,
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( 'DSBTRD', -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 DLASET( '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 DLARGV( 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: * DLARTV or DROT is used
196: *
197: IF( NR.GE.2*KD-1 ) THEN
198: DO 10 L = 1, KD - 1
199: CALL DLARTV( 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 DROT( 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 DLARTG( 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 DROT( 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 DLAR2V( 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: * DLARTV or DROT 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 DLARTV( 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 DROT( 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 DROT( 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 DROT( 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 DROT( 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 DLARGV( 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: * DLARTV or DROT is used
388: *
389: IF( NR.GT.2*KD-1 ) THEN
390: DO 130 L = 1, KD - 1
391: CALL DLARTV( 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 DROT( 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 DLARTG( 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 DROT( 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 DLAR2V( 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: * DLARTV or DROT 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 DLARTV( 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 DROT( 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 DROT( 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 DROT( 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 DROT( 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 DSBTRD
552: *
553: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsbtrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack