LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

◆ dsbtrd()

 subroutine dsbtrd ( character VECT, character UPLO, integer N, integer KD, double precision, dimension( ldab, * ) AB, integer LDAB, double precision, dimension( * ) D, double precision, dimension( * ) E, double precision, dimension( ldq, * ) Q, integer LDQ, double precision, dimension( * ) WORK, integer INFO )

DSBTRD

Purpose:
``` DSBTRD reduces a real symmetric band matrix A to symmetric
tridiagonal form T by an orthogonal similarity transformation:
Q**T * A * Q = T.```
Parameters
 [in] VECT ``` VECT is CHARACTER*1 = 'N': do not form Q; = 'V': form Q; = 'U': update a matrix X, by forming X*Q.``` [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] KD ``` KD is INTEGER The number of superdiagonals of the matrix A if UPLO = 'U', or the number of subdiagonals if UPLO = 'L'. KD >= 0.``` [in,out] AB ``` AB is DOUBLE PRECISION array, dimension (LDAB,N) On entry, the upper or lower triangle of the symmetric band matrix A, stored in the first KD+1 rows of the array. The j-th column of A is stored in the j-th column of the array AB as follows: if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). On exit, the diagonal elements of AB are overwritten by the diagonal elements of the tridiagonal matrix T; if KD > 0, the elements on the first superdiagonal (if UPLO = 'U') or the first subdiagonal (if UPLO = 'L') are overwritten by the off-diagonal elements of T; the rest of AB is overwritten by values generated during the reduction.``` [in] LDAB ``` LDAB is INTEGER The leading dimension of the array AB. LDAB >= KD+1.``` [out] D ``` D is DOUBLE PRECISION array, dimension (N) The diagonal elements of the tridiagonal matrix T.``` [out] E ``` E is DOUBLE PRECISION array, dimension (N-1) The off-diagonal elements of the tridiagonal matrix T: E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.``` [in,out] Q ``` Q is DOUBLE PRECISION array, dimension (LDQ,N) On entry, if VECT = 'U', then Q must contain an N-by-N matrix X; if VECT = 'N' or 'V', then Q need not be set. On exit: if VECT = 'V', Q contains the N-by-N orthogonal matrix Q; if VECT = 'U', Q contains the product X*Q; if VECT = 'N', the array Q is not referenced.``` [in] LDQ ``` LDQ is INTEGER The leading dimension of the array Q. LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.``` [out] WORK ` WORK is DOUBLE PRECISION array, dimension (N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
Further Details:
`  Modified by Linda Kaufman, Bell Labs.`

Definition at line 161 of file dsbtrd.f.

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 *
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
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
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
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
Here is the call graph for this function:
Here is the caller graph for this function: