LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \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