LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
dsbtrd.f
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
dsbtrd
subroutine dsbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
DSBTRD
Definition: dsbtrd.f:163
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
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
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
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