LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsbtrd.f
Go to the documentation of this file.
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 hbtrd
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
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dsbtrd
subroutine dsbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
DSBTRD
Definition dsbtrd.f:163
dlar2v
subroutine dlar2v(n, x, y, z, incx, c, s, incc)
DLAR2V applies a vector of plane rotations with real cosines and real sines from both sides to a sequ...
Definition dlar2v.f:110
dlargv
subroutine dlargv(n, x, incx, y, incy, c, incc)
DLARGV generates a vector of plane rotations with real cosines and real sines.
Definition dlargv.f:104
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlartv
subroutine dlartv(n, x, incx, y, incy, c, s, incc)
DLARTV applies a vector of plane rotations with real cosines and real sines to the elements of a pair...
Definition dlartv.f:108
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92