1: *> \brief \b DSBTRD
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSBTRD + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbtrd.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbtrd.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbtrd.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
22: * WORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER UPLO, VECT
26: * INTEGER INFO, KD, LDAB, LDQ, N
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
30: * $ WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DSBTRD reduces a real symmetric band matrix A to symmetric
40: *> tridiagonal form T by an orthogonal similarity transformation:
41: *> Q**T * A * Q = T.
42: *> \endverbatim
43: *
44: * Arguments:
45: * ==========
46: *
47: *> \param[in] VECT
48: *> \verbatim
49: *> VECT is CHARACTER*1
50: *> = 'N': do not form Q;
51: *> = 'V': form Q;
52: *> = 'U': update a matrix X, by forming X*Q.
53: *> \endverbatim
54: *>
55: *> \param[in] UPLO
56: *> \verbatim
57: *> UPLO is CHARACTER*1
58: *> = 'U': Upper triangle of A is stored;
59: *> = 'L': Lower triangle of A is stored.
60: *> \endverbatim
61: *>
62: *> \param[in] N
63: *> \verbatim
64: *> N is INTEGER
65: *> The order of the matrix A. N >= 0.
66: *> \endverbatim
67: *>
68: *> \param[in] KD
69: *> \verbatim
70: *> KD is INTEGER
71: *> The number of superdiagonals of the matrix A if UPLO = 'U',
72: *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
73: *> \endverbatim
74: *>
75: *> \param[in,out] AB
76: *> \verbatim
77: *> AB is DOUBLE PRECISION array, dimension (LDAB,N)
78: *> On entry, the upper or lower triangle of the symmetric band
79: *> matrix A, stored in the first KD+1 rows of the array. The
80: *> j-th column of A is stored in the j-th column of the array AB
81: *> as follows:
82: *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
83: *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
84: *> On exit, the diagonal elements of AB are overwritten by the
85: *> diagonal elements of the tridiagonal matrix T; if KD > 0, the
86: *> elements on the first superdiagonal (if UPLO = 'U') or the
87: *> first subdiagonal (if UPLO = 'L') are overwritten by the
88: *> off-diagonal elements of T; the rest of AB is overwritten by
89: *> values generated during the reduction.
90: *> \endverbatim
91: *>
92: *> \param[in] LDAB
93: *> \verbatim
94: *> LDAB is INTEGER
95: *> The leading dimension of the array AB. LDAB >= KD+1.
96: *> \endverbatim
97: *>
98: *> \param[out] D
99: *> \verbatim
100: *> D is DOUBLE PRECISION array, dimension (N)
101: *> The diagonal elements of the tridiagonal matrix T.
102: *> \endverbatim
103: *>
104: *> \param[out] E
105: *> \verbatim
106: *> E is DOUBLE PRECISION array, dimension (N-1)
107: *> The off-diagonal elements of the tridiagonal matrix T:
108: *> E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
109: *> \endverbatim
110: *>
111: *> \param[in,out] Q
112: *> \verbatim
113: *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
114: *> On entry, if VECT = 'U', then Q must contain an N-by-N
115: *> matrix X; if VECT = 'N' or 'V', then Q need not be set.
116: *>
117: *> On exit:
118: *> if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
119: *> if VECT = 'U', Q contains the product X*Q;
120: *> if VECT = 'N', the array Q is not referenced.
121: *> \endverbatim
122: *>
123: *> \param[in] LDQ
124: *> \verbatim
125: *> LDQ is INTEGER
126: *> The leading dimension of the array Q.
127: *> LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
128: *> \endverbatim
129: *>
130: *> \param[out] WORK
131: *> \verbatim
132: *> WORK is DOUBLE PRECISION array, dimension (N)
133: *> \endverbatim
134: *>
135: *> \param[out] INFO
136: *> \verbatim
137: *> INFO is INTEGER
138: *> = 0: successful exit
139: *> < 0: if INFO = -i, the i-th argument had an illegal value
140: *> \endverbatim
141: *
142: * Authors:
143: * ========
144: *
145: *> \author Univ. of Tennessee
146: *> \author Univ. of California Berkeley
147: *> \author Univ. of Colorado Denver
148: *> \author NAG Ltd.
149: *
150: *> \ingroup doubleOTHERcomputational
151: *
152: *> \par Further Details:
153: * =====================
154: *>
155: *> \verbatim
156: *>
157: *> Modified by Linda Kaufman, Bell Labs.
158: *> \endverbatim
159: *>
160: * =====================================================================
161: SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
162: $ WORK, INFO )
163: *
164: * -- LAPACK computational routine --
165: * -- LAPACK is a software package provided by Univ. of Tennessee, --
166: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167: *
168: * .. Scalar Arguments ..
169: CHARACTER UPLO, VECT
170: INTEGER INFO, KD, LDAB, LDQ, N
171: * ..
172: * .. Array Arguments ..
173: DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
174: $ WORK( * )
175: * ..
176: *
177: * =====================================================================
178: *
179: * .. Parameters ..
180: DOUBLE PRECISION ZERO, ONE
181: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
182: * ..
183: * .. Local Scalars ..
184: LOGICAL INITQ, UPPER, WANTQ
185: INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
186: $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
187: $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
188: DOUBLE PRECISION TEMP
189: * ..
190: * .. External Subroutines ..
191: EXTERNAL DLAR2V, DLARGV, DLARTG, DLARTV, DLASET, DROT,
192: $ XERBLA
193: * ..
194: * .. Intrinsic Functions ..
195: INTRINSIC MAX, MIN
196: * ..
197: * .. External Functions ..
198: LOGICAL LSAME
199: EXTERNAL LSAME
200: * ..
201: * .. Executable Statements ..
202: *
203: * Test the input parameters
204: *
205: INITQ = LSAME( VECT, 'V' )
206: WANTQ = INITQ .OR. LSAME( VECT, 'U' )
207: UPPER = LSAME( UPLO, 'U' )
208: KD1 = KD + 1
209: KDM1 = KD - 1
210: INCX = LDAB - 1
211: IQEND = 1
212: *
213: INFO = 0
214: IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
215: INFO = -1
216: ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
217: INFO = -2
218: ELSE IF( N.LT.0 ) THEN
219: INFO = -3
220: ELSE IF( KD.LT.0 ) THEN
221: INFO = -4
222: ELSE IF( LDAB.LT.KD1 ) THEN
223: INFO = -6
224: ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
225: INFO = -10
226: END IF
227: IF( INFO.NE.0 ) THEN
228: CALL XERBLA( 'DSBTRD', -INFO )
229: RETURN
230: END IF
231: *
232: * Quick return if possible
233: *
234: IF( N.EQ.0 )
235: $ RETURN
236: *
237: * Initialize Q to the unit matrix, if needed
238: *
239: IF( INITQ )
240: $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
241: *
242: * Wherever possible, plane rotations are generated and applied in
243: * vector operations of length NR over the index set J1:J2:KD1.
244: *
245: * The cosines and sines of the plane rotations are stored in the
246: * arrays D and WORK.
247: *
248: INCA = KD1*LDAB
249: KDN = MIN( N-1, KD )
250: IF( UPPER ) THEN
251: *
252: IF( KD.GT.1 ) THEN
253: *
254: * Reduce to tridiagonal form, working with upper triangle
255: *
256: NR = 0
257: J1 = KDN + 2
258: J2 = 1
259: *
260: DO 90 I = 1, N - 2
261: *
262: * Reduce i-th row of matrix to tridiagonal form
263: *
264: DO 80 K = KDN + 1, 2, -1
265: J1 = J1 + KDN
266: J2 = J2 + KDN
267: *
268: IF( NR.GT.0 ) THEN
269: *
270: * generate plane rotations to annihilate nonzero
271: * elements which have been created outside the band
272: *
273: CALL DLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
274: $ KD1, D( J1 ), KD1 )
275: *
276: * apply rotations from the right
277: *
278: *
279: * Dependent on the the number of diagonals either
280: * DLARTV or DROT is used
281: *
282: IF( NR.GE.2*KD-1 ) THEN
283: DO 10 L = 1, KD - 1
284: CALL DLARTV( NR, AB( L+1, J1-1 ), INCA,
285: $ AB( L, J1 ), INCA, D( J1 ),
286: $ WORK( J1 ), KD1 )
287: 10 CONTINUE
288: *
289: ELSE
290: JEND = J1 + ( NR-1 )*KD1
291: DO 20 JINC = J1, JEND, KD1
292: CALL DROT( KDM1, AB( 2, JINC-1 ), 1,
293: $ AB( 1, JINC ), 1, D( JINC ),
294: $ WORK( JINC ) )
295: 20 CONTINUE
296: END IF
297: END IF
298: *
299: *
300: IF( K.GT.2 ) THEN
301: IF( K.LE.N-I+1 ) THEN
302: *
303: * generate plane rotation to annihilate a(i,i+k-1)
304: * within the band
305: *
306: CALL DLARTG( AB( KD-K+3, I+K-2 ),
307: $ AB( KD-K+2, I+K-1 ), D( I+K-1 ),
308: $ WORK( I+K-1 ), TEMP )
309: AB( KD-K+3, I+K-2 ) = TEMP
310: *
311: * apply rotation from the right
312: *
313: CALL DROT( K-3, AB( KD-K+4, I+K-2 ), 1,
314: $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
315: $ WORK( I+K-1 ) )
316: END IF
317: NR = NR + 1
318: J1 = J1 - KDN - 1
319: END IF
320: *
321: * apply plane rotations from both sides to diagonal
322: * blocks
323: *
324: IF( NR.GT.0 )
325: $ CALL DLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
326: $ AB( KD, J1 ), INCA, D( J1 ),
327: $ WORK( J1 ), KD1 )
328: *
329: * apply plane rotations from the left
330: *
331: IF( NR.GT.0 ) THEN
332: IF( 2*KD-1.LT.NR ) THEN
333: *
334: * Dependent on the the number of diagonals either
335: * DLARTV or DROT is used
336: *
337: DO 30 L = 1, KD - 1
338: IF( J2+L.GT.N ) THEN
339: NRT = NR - 1
340: ELSE
341: NRT = NR
342: END IF
343: IF( NRT.GT.0 )
344: $ CALL DLARTV( NRT, AB( KD-L, J1+L ), INCA,
345: $ AB( KD-L+1, J1+L ), INCA,
346: $ D( J1 ), WORK( J1 ), KD1 )
347: 30 CONTINUE
348: ELSE
349: J1END = J1 + KD1*( NR-2 )
350: IF( J1END.GE.J1 ) THEN
351: DO 40 JIN = J1, J1END, KD1
352: CALL DROT( KD-1, AB( KD-1, JIN+1 ), INCX,
353: $ AB( KD, JIN+1 ), INCX,
354: $ D( JIN ), WORK( JIN ) )
355: 40 CONTINUE
356: END IF
357: LEND = MIN( KDM1, N-J2 )
358: LAST = J1END + KD1
359: IF( LEND.GT.0 )
360: $ CALL DROT( LEND, AB( KD-1, LAST+1 ), INCX,
361: $ AB( KD, LAST+1 ), INCX, D( LAST ),
362: $ WORK( LAST ) )
363: END IF
364: END IF
365: *
366: IF( WANTQ ) THEN
367: *
368: * accumulate product of plane rotations in Q
369: *
370: IF( INITQ ) THEN
371: *
372: * take advantage of the fact that Q was
373: * initially the Identity matrix
374: *
375: IQEND = MAX( IQEND, J2 )
376: I2 = MAX( 0, K-3 )
377: IQAEND = 1 + I*KD
378: IF( K.EQ.2 )
379: $ IQAEND = IQAEND + KD
380: IQAEND = MIN( IQAEND, IQEND )
381: DO 50 J = J1, J2, KD1
382: IBL = I - I2 / KDM1
383: I2 = I2 + 1
384: IQB = MAX( 1, J-IBL )
385: NQ = 1 + IQAEND - IQB
386: IQAEND = MIN( IQAEND+KD, IQEND )
387: CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
388: $ 1, D( J ), WORK( J ) )
389: 50 CONTINUE
390: ELSE
391: *
392: DO 60 J = J1, J2, KD1
393: CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
394: $ D( J ), WORK( J ) )
395: 60 CONTINUE
396: END IF
397: *
398: END IF
399: *
400: IF( J2+KDN.GT.N ) THEN
401: *
402: * adjust J2 to keep within the bounds of the matrix
403: *
404: NR = NR - 1
405: J2 = J2 - KDN - 1
406: END IF
407: *
408: DO 70 J = J1, J2, KD1
409: *
410: * create nonzero element a(j-1,j+kd) outside the band
411: * and store it in WORK
412: *
413: WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
414: AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
415: 70 CONTINUE
416: 80 CONTINUE
417: 90 CONTINUE
418: END IF
419: *
420: IF( KD.GT.0 ) THEN
421: *
422: * copy off-diagonal elements to E
423: *
424: DO 100 I = 1, N - 1
425: E( I ) = AB( KD, I+1 )
426: 100 CONTINUE
427: ELSE
428: *
429: * set E to zero if original matrix was diagonal
430: *
431: DO 110 I = 1, N - 1
432: E( I ) = ZERO
433: 110 CONTINUE
434: END IF
435: *
436: * copy diagonal elements to D
437: *
438: DO 120 I = 1, N
439: D( I ) = AB( KD1, I )
440: 120 CONTINUE
441: *
442: ELSE
443: *
444: IF( KD.GT.1 ) THEN
445: *
446: * Reduce to tridiagonal form, working with lower triangle
447: *
448: NR = 0
449: J1 = KDN + 2
450: J2 = 1
451: *
452: DO 210 I = 1, N - 2
453: *
454: * Reduce i-th column of matrix to tridiagonal form
455: *
456: DO 200 K = KDN + 1, 2, -1
457: J1 = J1 + KDN
458: J2 = J2 + KDN
459: *
460: IF( NR.GT.0 ) THEN
461: *
462: * generate plane rotations to annihilate nonzero
463: * elements which have been created outside the band
464: *
465: CALL DLARGV( NR, AB( KD1, J1-KD1 ), INCA,
466: $ WORK( J1 ), KD1, D( J1 ), KD1 )
467: *
468: * apply plane rotations from one side
469: *
470: *
471: * Dependent on the the number of diagonals either
472: * DLARTV or DROT is used
473: *
474: IF( NR.GT.2*KD-1 ) THEN
475: DO 130 L = 1, KD - 1
476: CALL DLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
477: $ AB( KD1-L+1, J1-KD1+L ), INCA,
478: $ D( J1 ), WORK( J1 ), KD1 )
479: 130 CONTINUE
480: ELSE
481: JEND = J1 + KD1*( NR-1 )
482: DO 140 JINC = J1, JEND, KD1
483: CALL DROT( KDM1, AB( KD, JINC-KD ), INCX,
484: $ AB( KD1, JINC-KD ), INCX,
485: $ D( JINC ), WORK( JINC ) )
486: 140 CONTINUE
487: END IF
488: *
489: END IF
490: *
491: IF( K.GT.2 ) THEN
492: IF( K.LE.N-I+1 ) THEN
493: *
494: * generate plane rotation to annihilate a(i+k-1,i)
495: * within the band
496: *
497: CALL DLARTG( AB( K-1, I ), AB( K, I ),
498: $ D( I+K-1 ), WORK( I+K-1 ), TEMP )
499: AB( K-1, I ) = TEMP
500: *
501: * apply rotation from the left
502: *
503: CALL DROT( K-3, AB( K-2, I+1 ), LDAB-1,
504: $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
505: $ WORK( I+K-1 ) )
506: END IF
507: NR = NR + 1
508: J1 = J1 - KDN - 1
509: END IF
510: *
511: * apply plane rotations from both sides to diagonal
512: * blocks
513: *
514: IF( NR.GT.0 )
515: $ CALL DLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
516: $ AB( 2, J1-1 ), INCA, D( J1 ),
517: $ WORK( J1 ), KD1 )
518: *
519: * apply plane rotations from the right
520: *
521: *
522: * Dependent on the the number of diagonals either
523: * DLARTV or DROT is used
524: *
525: IF( NR.GT.0 ) THEN
526: IF( NR.GT.2*KD-1 ) THEN
527: DO 150 L = 1, KD - 1
528: IF( J2+L.GT.N ) THEN
529: NRT = NR - 1
530: ELSE
531: NRT = NR
532: END IF
533: IF( NRT.GT.0 )
534: $ CALL DLARTV( NRT, AB( L+2, J1-1 ), INCA,
535: $ AB( L+1, J1 ), INCA, D( J1 ),
536: $ WORK( J1 ), KD1 )
537: 150 CONTINUE
538: ELSE
539: J1END = J1 + KD1*( NR-2 )
540: IF( J1END.GE.J1 ) THEN
541: DO 160 J1INC = J1, J1END, KD1
542: CALL DROT( KDM1, AB( 3, J1INC-1 ), 1,
543: $ AB( 2, J1INC ), 1, D( J1INC ),
544: $ WORK( J1INC ) )
545: 160 CONTINUE
546: END IF
547: LEND = MIN( KDM1, N-J2 )
548: LAST = J1END + KD1
549: IF( LEND.GT.0 )
550: $ CALL DROT( LEND, AB( 3, LAST-1 ), 1,
551: $ AB( 2, LAST ), 1, D( LAST ),
552: $ WORK( LAST ) )
553: END IF
554: END IF
555: *
556: *
557: *
558: IF( WANTQ ) THEN
559: *
560: * accumulate product of plane rotations in Q
561: *
562: IF( INITQ ) THEN
563: *
564: * take advantage of the fact that Q was
565: * initially the Identity matrix
566: *
567: IQEND = MAX( IQEND, J2 )
568: I2 = MAX( 0, K-3 )
569: IQAEND = 1 + I*KD
570: IF( K.EQ.2 )
571: $ IQAEND = IQAEND + KD
572: IQAEND = MIN( IQAEND, IQEND )
573: DO 170 J = J1, J2, KD1
574: IBL = I - I2 / KDM1
575: I2 = I2 + 1
576: IQB = MAX( 1, J-IBL )
577: NQ = 1 + IQAEND - IQB
578: IQAEND = MIN( IQAEND+KD, IQEND )
579: CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
580: $ 1, D( J ), WORK( J ) )
581: 170 CONTINUE
582: ELSE
583: *
584: DO 180 J = J1, J2, KD1
585: CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
586: $ D( J ), WORK( J ) )
587: 180 CONTINUE
588: END IF
589: END IF
590: *
591: IF( J2+KDN.GT.N ) THEN
592: *
593: * adjust J2 to keep within the bounds of the matrix
594: *
595: NR = NR - 1
596: J2 = J2 - KDN - 1
597: END IF
598: *
599: DO 190 J = J1, J2, KD1
600: *
601: * create nonzero element a(j+kd,j-1) outside the
602: * band and store it in WORK
603: *
604: WORK( J+KD ) = WORK( J )*AB( KD1, J )
605: AB( KD1, J ) = D( J )*AB( KD1, J )
606: 190 CONTINUE
607: 200 CONTINUE
608: 210 CONTINUE
609: END IF
610: *
611: IF( KD.GT.0 ) THEN
612: *
613: * copy off-diagonal elements to E
614: *
615: DO 220 I = 1, N - 1
616: E( I ) = AB( 2, I )
617: 220 CONTINUE
618: ELSE
619: *
620: * set E to zero if original matrix was diagonal
621: *
622: DO 230 I = 1, N - 1
623: E( I ) = ZERO
624: 230 CONTINUE
625: END IF
626: *
627: * copy diagonal elements to D
628: *
629: DO 240 I = 1, N
630: D( I ) = AB( 1, I )
631: 240 CONTINUE
632: END IF
633: *
634: RETURN
635: *
636: * End of DSBTRD
637: *
638: 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