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: *> \date November 2011
151: *
152: *> \ingroup doubleOTHERcomputational
153: *
154: *> \par Further Details:
155: * =====================
156: *>
157: *> \verbatim
158: *>
159: *> Modified by Linda Kaufman, Bell Labs.
160: *> \endverbatim
161: *>
162: * =====================================================================
163: SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
164: $ WORK, INFO )
165: *
166: * -- LAPACK computational routine (version 3.4.0) --
167: * -- LAPACK is a software package provided by Univ. of Tennessee, --
168: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169: * November 2011
170: *
171: * .. Scalar Arguments ..
172: CHARACTER UPLO, VECT
173: INTEGER INFO, KD, LDAB, LDQ, N
174: * ..
175: * .. Array Arguments ..
176: DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
177: $ WORK( * )
178: * ..
179: *
180: * =====================================================================
181: *
182: * .. Parameters ..
183: DOUBLE PRECISION ZERO, ONE
184: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
185: * ..
186: * .. Local Scalars ..
187: LOGICAL INITQ, UPPER, WANTQ
188: INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
189: $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
190: $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
191: DOUBLE PRECISION TEMP
192: * ..
193: * .. External Subroutines ..
194: EXTERNAL DLAR2V, DLARGV, DLARTG, DLARTV, DLASET, DROT,
195: $ XERBLA
196: * ..
197: * .. Intrinsic Functions ..
198: INTRINSIC MAX, MIN
199: * ..
200: * .. External Functions ..
201: LOGICAL LSAME
202: EXTERNAL LSAME
203: * ..
204: * .. Executable Statements ..
205: *
206: * Test the input parameters
207: *
208: INITQ = LSAME( VECT, 'V' )
209: WANTQ = INITQ .OR. LSAME( VECT, 'U' )
210: UPPER = LSAME( UPLO, 'U' )
211: KD1 = KD + 1
212: KDM1 = KD - 1
213: INCX = LDAB - 1
214: IQEND = 1
215: *
216: INFO = 0
217: IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
218: INFO = -1
219: ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
220: INFO = -2
221: ELSE IF( N.LT.0 ) THEN
222: INFO = -3
223: ELSE IF( KD.LT.0 ) THEN
224: INFO = -4
225: ELSE IF( LDAB.LT.KD1 ) THEN
226: INFO = -6
227: ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
228: INFO = -10
229: END IF
230: IF( INFO.NE.0 ) THEN
231: CALL XERBLA( 'DSBTRD', -INFO )
232: RETURN
233: END IF
234: *
235: * Quick return if possible
236: *
237: IF( N.EQ.0 )
238: $ RETURN
239: *
240: * Initialize Q to the unit matrix, if needed
241: *
242: IF( INITQ )
243: $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
244: *
245: * Wherever possible, plane rotations are generated and applied in
246: * vector operations of length NR over the index set J1:J2:KD1.
247: *
248: * The cosines and sines of the plane rotations are stored in the
249: * arrays D and WORK.
250: *
251: INCA = KD1*LDAB
252: KDN = MIN( N-1, KD )
253: IF( UPPER ) THEN
254: *
255: IF( KD.GT.1 ) THEN
256: *
257: * Reduce to tridiagonal form, working with upper triangle
258: *
259: NR = 0
260: J1 = KDN + 2
261: J2 = 1
262: *
263: DO 90 I = 1, N - 2
264: *
265: * Reduce i-th row of matrix to tridiagonal form
266: *
267: DO 80 K = KDN + 1, 2, -1
268: J1 = J1 + KDN
269: J2 = J2 + KDN
270: *
271: IF( NR.GT.0 ) THEN
272: *
273: * generate plane rotations to annihilate nonzero
274: * elements which have been created outside the band
275: *
276: CALL DLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
277: $ KD1, D( J1 ), KD1 )
278: *
279: * apply rotations from the right
280: *
281: *
282: * Dependent on the the number of diagonals either
283: * DLARTV or DROT is used
284: *
285: IF( NR.GE.2*KD-1 ) THEN
286: DO 10 L = 1, KD - 1
287: CALL DLARTV( NR, AB( L+1, J1-1 ), INCA,
288: $ AB( L, J1 ), INCA, D( J1 ),
289: $ WORK( J1 ), KD1 )
290: 10 CONTINUE
291: *
292: ELSE
293: JEND = J1 + ( NR-1 )*KD1
294: DO 20 JINC = J1, JEND, KD1
295: CALL DROT( KDM1, AB( 2, JINC-1 ), 1,
296: $ AB( 1, JINC ), 1, D( JINC ),
297: $ WORK( JINC ) )
298: 20 CONTINUE
299: END IF
300: END IF
301: *
302: *
303: IF( K.GT.2 ) THEN
304: IF( K.LE.N-I+1 ) THEN
305: *
306: * generate plane rotation to annihilate a(i,i+k-1)
307: * within the band
308: *
309: CALL DLARTG( AB( KD-K+3, I+K-2 ),
310: $ AB( KD-K+2, I+K-1 ), D( I+K-1 ),
311: $ WORK( I+K-1 ), TEMP )
312: AB( KD-K+3, I+K-2 ) = TEMP
313: *
314: * apply rotation from the right
315: *
316: CALL DROT( K-3, AB( KD-K+4, I+K-2 ), 1,
317: $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
318: $ WORK( I+K-1 ) )
319: END IF
320: NR = NR + 1
321: J1 = J1 - KDN - 1
322: END IF
323: *
324: * apply plane rotations from both sides to diagonal
325: * blocks
326: *
327: IF( NR.GT.0 )
328: $ CALL DLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
329: $ AB( KD, J1 ), INCA, D( J1 ),
330: $ WORK( J1 ), KD1 )
331: *
332: * apply plane rotations from the left
333: *
334: IF( NR.GT.0 ) THEN
335: IF( 2*KD-1.LT.NR ) THEN
336: *
337: * Dependent on the the number of diagonals either
338: * DLARTV or DROT is used
339: *
340: DO 30 L = 1, KD - 1
341: IF( J2+L.GT.N ) THEN
342: NRT = NR - 1
343: ELSE
344: NRT = NR
345: END IF
346: IF( NRT.GT.0 )
347: $ CALL DLARTV( NRT, AB( KD-L, J1+L ), INCA,
348: $ AB( KD-L+1, J1+L ), INCA,
349: $ D( J1 ), WORK( J1 ), KD1 )
350: 30 CONTINUE
351: ELSE
352: J1END = J1 + KD1*( NR-2 )
353: IF( J1END.GE.J1 ) THEN
354: DO 40 JIN = J1, J1END, KD1
355: CALL DROT( KD-1, AB( KD-1, JIN+1 ), INCX,
356: $ AB( KD, JIN+1 ), INCX,
357: $ D( JIN ), WORK( JIN ) )
358: 40 CONTINUE
359: END IF
360: LEND = MIN( KDM1, N-J2 )
361: LAST = J1END + KD1
362: IF( LEND.GT.0 )
363: $ CALL DROT( LEND, AB( KD-1, LAST+1 ), INCX,
364: $ AB( KD, LAST+1 ), INCX, D( LAST ),
365: $ WORK( LAST ) )
366: END IF
367: END IF
368: *
369: IF( WANTQ ) THEN
370: *
371: * accumulate product of plane rotations in Q
372: *
373: IF( INITQ ) THEN
374: *
375: * take advantage of the fact that Q was
376: * initially the Identity matrix
377: *
378: IQEND = MAX( IQEND, J2 )
379: I2 = MAX( 0, K-3 )
380: IQAEND = 1 + I*KD
381: IF( K.EQ.2 )
382: $ IQAEND = IQAEND + KD
383: IQAEND = MIN( IQAEND, IQEND )
384: DO 50 J = J1, J2, KD1
385: IBL = I - I2 / KDM1
386: I2 = I2 + 1
387: IQB = MAX( 1, J-IBL )
388: NQ = 1 + IQAEND - IQB
389: IQAEND = MIN( IQAEND+KD, IQEND )
390: CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
391: $ 1, D( J ), WORK( J ) )
392: 50 CONTINUE
393: ELSE
394: *
395: DO 60 J = J1, J2, KD1
396: CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
397: $ D( J ), WORK( J ) )
398: 60 CONTINUE
399: END IF
400: *
401: END IF
402: *
403: IF( J2+KDN.GT.N ) THEN
404: *
405: * adjust J2 to keep within the bounds of the matrix
406: *
407: NR = NR - 1
408: J2 = J2 - KDN - 1
409: END IF
410: *
411: DO 70 J = J1, J2, KD1
412: *
413: * create nonzero element a(j-1,j+kd) outside the band
414: * and store it in WORK
415: *
416: WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
417: AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
418: 70 CONTINUE
419: 80 CONTINUE
420: 90 CONTINUE
421: END IF
422: *
423: IF( KD.GT.0 ) THEN
424: *
425: * copy off-diagonal elements to E
426: *
427: DO 100 I = 1, N - 1
428: E( I ) = AB( KD, I+1 )
429: 100 CONTINUE
430: ELSE
431: *
432: * set E to zero if original matrix was diagonal
433: *
434: DO 110 I = 1, N - 1
435: E( I ) = ZERO
436: 110 CONTINUE
437: END IF
438: *
439: * copy diagonal elements to D
440: *
441: DO 120 I = 1, N
442: D( I ) = AB( KD1, I )
443: 120 CONTINUE
444: *
445: ELSE
446: *
447: IF( KD.GT.1 ) THEN
448: *
449: * Reduce to tridiagonal form, working with lower triangle
450: *
451: NR = 0
452: J1 = KDN + 2
453: J2 = 1
454: *
455: DO 210 I = 1, N - 2
456: *
457: * Reduce i-th column of matrix to tridiagonal form
458: *
459: DO 200 K = KDN + 1, 2, -1
460: J1 = J1 + KDN
461: J2 = J2 + KDN
462: *
463: IF( NR.GT.0 ) THEN
464: *
465: * generate plane rotations to annihilate nonzero
466: * elements which have been created outside the band
467: *
468: CALL DLARGV( NR, AB( KD1, J1-KD1 ), INCA,
469: $ WORK( J1 ), KD1, D( J1 ), KD1 )
470: *
471: * apply plane rotations from one side
472: *
473: *
474: * Dependent on the the number of diagonals either
475: * DLARTV or DROT is used
476: *
477: IF( NR.GT.2*KD-1 ) THEN
478: DO 130 L = 1, KD - 1
479: CALL DLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
480: $ AB( KD1-L+1, J1-KD1+L ), INCA,
481: $ D( J1 ), WORK( J1 ), KD1 )
482: 130 CONTINUE
483: ELSE
484: JEND = J1 + KD1*( NR-1 )
485: DO 140 JINC = J1, JEND, KD1
486: CALL DROT( KDM1, AB( KD, JINC-KD ), INCX,
487: $ AB( KD1, JINC-KD ), INCX,
488: $ D( JINC ), WORK( JINC ) )
489: 140 CONTINUE
490: END IF
491: *
492: END IF
493: *
494: IF( K.GT.2 ) THEN
495: IF( K.LE.N-I+1 ) THEN
496: *
497: * generate plane rotation to annihilate a(i+k-1,i)
498: * within the band
499: *
500: CALL DLARTG( AB( K-1, I ), AB( K, I ),
501: $ D( I+K-1 ), WORK( I+K-1 ), TEMP )
502: AB( K-1, I ) = TEMP
503: *
504: * apply rotation from the left
505: *
506: CALL DROT( K-3, AB( K-2, I+1 ), LDAB-1,
507: $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
508: $ WORK( I+K-1 ) )
509: END IF
510: NR = NR + 1
511: J1 = J1 - KDN - 1
512: END IF
513: *
514: * apply plane rotations from both sides to diagonal
515: * blocks
516: *
517: IF( NR.GT.0 )
518: $ CALL DLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
519: $ AB( 2, J1-1 ), INCA, D( J1 ),
520: $ WORK( J1 ), KD1 )
521: *
522: * apply plane rotations from the right
523: *
524: *
525: * Dependent on the the number of diagonals either
526: * DLARTV or DROT is used
527: *
528: IF( NR.GT.0 ) THEN
529: IF( NR.GT.2*KD-1 ) THEN
530: DO 150 L = 1, KD - 1
531: IF( J2+L.GT.N ) THEN
532: NRT = NR - 1
533: ELSE
534: NRT = NR
535: END IF
536: IF( NRT.GT.0 )
537: $ CALL DLARTV( NRT, AB( L+2, J1-1 ), INCA,
538: $ AB( L+1, J1 ), INCA, D( J1 ),
539: $ WORK( J1 ), KD1 )
540: 150 CONTINUE
541: ELSE
542: J1END = J1 + KD1*( NR-2 )
543: IF( J1END.GE.J1 ) THEN
544: DO 160 J1INC = J1, J1END, KD1
545: CALL DROT( KDM1, AB( 3, J1INC-1 ), 1,
546: $ AB( 2, J1INC ), 1, D( J1INC ),
547: $ WORK( J1INC ) )
548: 160 CONTINUE
549: END IF
550: LEND = MIN( KDM1, N-J2 )
551: LAST = J1END + KD1
552: IF( LEND.GT.0 )
553: $ CALL DROT( LEND, AB( 3, LAST-1 ), 1,
554: $ AB( 2, LAST ), 1, D( LAST ),
555: $ WORK( LAST ) )
556: END IF
557: END IF
558: *
559: *
560: *
561: IF( WANTQ ) THEN
562: *
563: * accumulate product of plane rotations in Q
564: *
565: IF( INITQ ) THEN
566: *
567: * take advantage of the fact that Q was
568: * initially the Identity matrix
569: *
570: IQEND = MAX( IQEND, J2 )
571: I2 = MAX( 0, K-3 )
572: IQAEND = 1 + I*KD
573: IF( K.EQ.2 )
574: $ IQAEND = IQAEND + KD
575: IQAEND = MIN( IQAEND, IQEND )
576: DO 170 J = J1, J2, KD1
577: IBL = I - I2 / KDM1
578: I2 = I2 + 1
579: IQB = MAX( 1, J-IBL )
580: NQ = 1 + IQAEND - IQB
581: IQAEND = MIN( IQAEND+KD, IQEND )
582: CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
583: $ 1, D( J ), WORK( J ) )
584: 170 CONTINUE
585: ELSE
586: *
587: DO 180 J = J1, J2, KD1
588: CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
589: $ D( J ), WORK( J ) )
590: 180 CONTINUE
591: END IF
592: END IF
593: *
594: IF( J2+KDN.GT.N ) THEN
595: *
596: * adjust J2 to keep within the bounds of the matrix
597: *
598: NR = NR - 1
599: J2 = J2 - KDN - 1
600: END IF
601: *
602: DO 190 J = J1, J2, KD1
603: *
604: * create nonzero element a(j+kd,j-1) outside the
605: * band and store it in WORK
606: *
607: WORK( J+KD ) = WORK( J )*AB( KD1, J )
608: AB( KD1, J ) = D( J )*AB( KD1, J )
609: 190 CONTINUE
610: 200 CONTINUE
611: 210 CONTINUE
612: END IF
613: *
614: IF( KD.GT.0 ) THEN
615: *
616: * copy off-diagonal elements to E
617: *
618: DO 220 I = 1, N - 1
619: E( I ) = AB( 2, I )
620: 220 CONTINUE
621: ELSE
622: *
623: * set E to zero if original matrix was diagonal
624: *
625: DO 230 I = 1, N - 1
626: E( I ) = ZERO
627: 230 CONTINUE
628: END IF
629: *
630: * copy diagonal elements to D
631: *
632: DO 240 I = 1, N
633: D( I ) = AB( 1, I )
634: 240 CONTINUE
635: END IF
636: *
637: RETURN
638: *
639: * End of DSBTRD
640: *
641: 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