LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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

Download DSBTRD + dependencies [TGZ] [ZIP] [TXT]

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
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Modified by Linda Kaufman, Bell Labs.
!>

Definition at line 159 of file dsbtrd.f.

161*
162* -- LAPACK computational routine --
163* -- LAPACK is a software package provided by Univ. of Tennessee, --
164* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165*
166* .. Scalar Arguments ..
167 CHARACTER UPLO, VECT
168 INTEGER INFO, KD, LDAB, LDQ, N
169* ..
170* .. Array Arguments ..
171 DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
172 $ WORK( * )
173* ..
174*
175* =====================================================================
176*
177* .. Parameters ..
178 DOUBLE PRECISION ZERO, ONE
179 parameter( zero = 0.0d+0, one = 1.0d+0 )
180* ..
181* .. Local Scalars ..
182 LOGICAL INITQ, UPPER, WANTQ
183 INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
184 $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
185 $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
186 DOUBLE PRECISION TEMP
187* ..
188* .. External Subroutines ..
189 EXTERNAL dlar2v, dlargv, dlartg, dlartv, dlaset,
190 $ drot,
191 $ xerbla
192* ..
193* .. Intrinsic Functions ..
194 INTRINSIC max, min
195* ..
196* .. External Functions ..
197 LOGICAL LSAME
198 EXTERNAL lsame
199* ..
200* .. Executable Statements ..
201*
202* Test the input parameters
203*
204 initq = lsame( vect, 'V' )
205 wantq = initq .OR. lsame( vect, 'U' )
206 upper = lsame( uplo, 'U' )
207 kd1 = kd + 1
208 kdm1 = kd - 1
209 incx = ldab - 1
210 iqend = 1
211*
212 info = 0
213 IF( .NOT.wantq .AND. .NOT.lsame( vect, 'N' ) ) THEN
214 info = -1
215 ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
216 info = -2
217 ELSE IF( n.LT.0 ) THEN
218 info = -3
219 ELSE IF( kd.LT.0 ) THEN
220 info = -4
221 ELSE IF( ldab.LT.kd1 ) THEN
222 info = -6
223 ELSE IF( ldq.LT.max( 1, n ) .AND. wantq ) THEN
224 info = -10
225 END IF
226 IF( info.NE.0 ) THEN
227 CALL xerbla( 'DSBTRD', -info )
228 RETURN
229 END IF
230*
231* Quick return if possible
232*
233 IF( n.EQ.0 )
234 $ RETURN
235*
236* Initialize Q to the unit matrix, if needed
237*
238 IF( initq )
239 $ CALL dlaset( 'Full', n, n, zero, one, q, ldq )
240*
241* Wherever possible, plane rotations are generated and applied in
242* vector operations of length NR over the index set J1:J2:KD1.
243*
244* The cosines and sines of the plane rotations are stored in the
245* arrays D and WORK.
246*
247 inca = kd1*ldab
248 kdn = min( n-1, kd )
249 IF( upper ) THEN
250*
251 IF( kd.GT.1 ) THEN
252*
253* Reduce to tridiagonal form, working with upper triangle
254*
255 nr = 0
256 j1 = kdn + 2
257 j2 = 1
258*
259 DO 90 i = 1, n - 2
260*
261* Reduce i-th row of matrix to tridiagonal form
262*
263 DO 80 k = kdn + 1, 2, -1
264 j1 = j1 + kdn
265 j2 = j2 + kdn
266*
267 IF( nr.GT.0 ) THEN
268*
269* generate plane rotations to annihilate nonzero
270* elements which have been created outside the band
271*
272 CALL dlargv( nr, ab( 1, j1-1 ), inca,
273 $ 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 ),
345 $ inca,
346 $ ab( kd-l+1, j1+l ), inca,
347 $ d( j1 ), work( j1 ), kd1 )
348 30 CONTINUE
349 ELSE
350 j1end = j1 + kd1*( nr-2 )
351 IF( j1end.GE.j1 ) THEN
352 DO 40 jin = j1, j1end, kd1
353 CALL drot( kd-1, ab( kd-1, jin+1 ),
354 $ incx,
355 $ ab( kd, jin+1 ), incx,
356 $ d( jin ), work( jin ) )
357 40 CONTINUE
358 END IF
359 lend = min( kdm1, n-j2 )
360 last = j1end + kd1
361 IF( lend.GT.0 )
362 $ CALL drot( lend, ab( kd-1, last+1 ), incx,
363 $ ab( kd, last+1 ), incx, d( last ),
364 $ work( last ) )
365 END IF
366 END IF
367*
368 IF( wantq ) THEN
369*
370* accumulate product of plane rotations in Q
371*
372 IF( initq ) THEN
373*
374* take advantage of the fact that Q was
375* initially the Identity matrix
376*
377 iqend = max( iqend, j2 )
378 i2 = max( 0, k-3 )
379 iqaend = 1 + i*kd
380 IF( k.EQ.2 )
381 $ iqaend = iqaend + kd
382 iqaend = min( iqaend, iqend )
383 DO 50 j = j1, j2, kd1
384 ibl = i - i2 / kdm1
385 i2 = i2 + 1
386 iqb = max( 1, j-ibl )
387 nq = 1 + iqaend - iqb
388 iqaend = min( iqaend+kd, iqend )
389 CALL drot( nq, q( iqb, j-1 ), 1, q( iqb,
390 $ 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 ),
397 $ 1,
398 $ d( j ), work( j ) )
399 60 CONTINUE
400 END IF
401*
402 END IF
403*
404 IF( j2+kdn.GT.n ) THEN
405*
406* adjust J2 to keep within the bounds of the matrix
407*
408 nr = nr - 1
409 j2 = j2 - kdn - 1
410 END IF
411*
412 DO 70 j = j1, j2, kd1
413*
414* create nonzero element a(j-1,j+kd) outside the band
415* and store it in WORK
416*
417 work( j+kd ) = work( j )*ab( 1, j+kd )
418 ab( 1, j+kd ) = d( j )*ab( 1, j+kd )
419 70 CONTINUE
420 80 CONTINUE
421 90 CONTINUE
422 END IF
423*
424 IF( kd.GT.0 ) THEN
425*
426* copy off-diagonal elements to E
427*
428 DO 100 i = 1, n - 1
429 e( i ) = ab( kd, i+1 )
430 100 CONTINUE
431 ELSE
432*
433* set E to zero if original matrix was diagonal
434*
435 DO 110 i = 1, n - 1
436 e( i ) = zero
437 110 CONTINUE
438 END IF
439*
440* copy diagonal elements to D
441*
442 DO 120 i = 1, n
443 d( i ) = ab( kd1, i )
444 120 CONTINUE
445*
446 ELSE
447*
448 IF( kd.GT.1 ) THEN
449*
450* Reduce to tridiagonal form, working with lower triangle
451*
452 nr = 0
453 j1 = kdn + 2
454 j2 = 1
455*
456 DO 210 i = 1, n - 2
457*
458* Reduce i-th column of matrix to tridiagonal form
459*
460 DO 200 k = kdn + 1, 2, -1
461 j1 = j1 + kdn
462 j2 = j2 + kdn
463*
464 IF( nr.GT.0 ) THEN
465*
466* generate plane rotations to annihilate nonzero
467* elements which have been created outside the band
468*
469 CALL dlargv( nr, ab( kd1, j1-kd1 ), inca,
470 $ work( j1 ), kd1, d( j1 ), kd1 )
471*
472* apply plane rotations from one side
473*
474*
475* Dependent on the the number of diagonals either
476* DLARTV or DROT is used
477*
478 IF( nr.GT.2*kd-1 ) THEN
479 DO 130 l = 1, kd - 1
480 CALL dlartv( nr, ab( kd1-l, j1-kd1+l ),
481 $ inca,
482 $ ab( kd1-l+1, j1-kd1+l ), inca,
483 $ d( j1 ), work( j1 ), kd1 )
484 130 CONTINUE
485 ELSE
486 jend = j1 + kd1*( nr-1 )
487 DO 140 jinc = j1, jend, kd1
488 CALL drot( kdm1, ab( kd, jinc-kd ), incx,
489 $ ab( kd1, jinc-kd ), incx,
490 $ d( jinc ), work( jinc ) )
491 140 CONTINUE
492 END IF
493*
494 END IF
495*
496 IF( k.GT.2 ) THEN
497 IF( k.LE.n-i+1 ) THEN
498*
499* generate plane rotation to annihilate a(i+k-1,i)
500* within the band
501*
502 CALL dlartg( ab( k-1, i ), ab( k, i ),
503 $ d( i+k-1 ), work( i+k-1 ), temp )
504 ab( k-1, i ) = temp
505*
506* apply rotation from the left
507*
508 CALL drot( k-3, ab( k-2, i+1 ), ldab-1,
509 $ ab( k-1, i+1 ), ldab-1, d( i+k-1 ),
510 $ work( i+k-1 ) )
511 END IF
512 nr = nr + 1
513 j1 = j1 - kdn - 1
514 END IF
515*
516* apply plane rotations from both sides to diagonal
517* blocks
518*
519 IF( nr.GT.0 )
520 $ CALL dlar2v( nr, ab( 1, j1-1 ), ab( 1, j1 ),
521 $ ab( 2, j1-1 ), inca, d( j1 ),
522 $ work( j1 ), kd1 )
523*
524* apply plane rotations from the right
525*
526*
527* Dependent on the the number of diagonals either
528* DLARTV or DROT is used
529*
530 IF( nr.GT.0 ) THEN
531 IF( nr.GT.2*kd-1 ) THEN
532 DO 150 l = 1, kd - 1
533 IF( j2+l.GT.n ) THEN
534 nrt = nr - 1
535 ELSE
536 nrt = nr
537 END IF
538 IF( nrt.GT.0 )
539 $ CALL dlartv( nrt, ab( l+2, j1-1 ),
540 $ inca,
541 $ ab( l+1, j1 ), inca, d( j1 ),
542 $ work( j1 ), kd1 )
543 150 CONTINUE
544 ELSE
545 j1end = j1 + kd1*( nr-2 )
546 IF( j1end.GE.j1 ) THEN
547 DO 160 j1inc = j1, j1end, kd1
548 CALL drot( kdm1, ab( 3, j1inc-1 ), 1,
549 $ ab( 2, j1inc ), 1, d( j1inc ),
550 $ work( j1inc ) )
551 160 CONTINUE
552 END IF
553 lend = min( kdm1, n-j2 )
554 last = j1end + kd1
555 IF( lend.GT.0 )
556 $ CALL drot( lend, ab( 3, last-1 ), 1,
557 $ ab( 2, last ), 1, d( last ),
558 $ work( last ) )
559 END IF
560 END IF
561*
562*
563*
564 IF( wantq ) THEN
565*
566* accumulate product of plane rotations in Q
567*
568 IF( initq ) THEN
569*
570* take advantage of the fact that Q was
571* initially the Identity matrix
572*
573 iqend = max( iqend, j2 )
574 i2 = max( 0, k-3 )
575 iqaend = 1 + i*kd
576 IF( k.EQ.2 )
577 $ iqaend = iqaend + kd
578 iqaend = min( iqaend, iqend )
579 DO 170 j = j1, j2, kd1
580 ibl = i - i2 / kdm1
581 i2 = i2 + 1
582 iqb = max( 1, j-ibl )
583 nq = 1 + iqaend - iqb
584 iqaend = min( iqaend+kd, iqend )
585 CALL drot( nq, q( iqb, j-1 ), 1, q( iqb,
586 $ j ),
587 $ 1, d( j ), work( j ) )
588 170 CONTINUE
589 ELSE
590*
591 DO 180 j = j1, j2, kd1
592 CALL drot( n, q( 1, j-1 ), 1, q( 1, j ),
593 $ 1,
594 $ d( j ), work( j ) )
595 180 CONTINUE
596 END IF
597 END IF
598*
599 IF( j2+kdn.GT.n ) THEN
600*
601* adjust J2 to keep within the bounds of the matrix
602*
603 nr = nr - 1
604 j2 = j2 - kdn - 1
605 END IF
606*
607 DO 190 j = j1, j2, kd1
608*
609* create nonzero element a(j+kd,j-1) outside the
610* band and store it in WORK
611*
612 work( j+kd ) = work( j )*ab( kd1, j )
613 ab( kd1, j ) = d( j )*ab( kd1, j )
614 190 CONTINUE
615 200 CONTINUE
616 210 CONTINUE
617 END IF
618*
619 IF( kd.GT.0 ) THEN
620*
621* copy off-diagonal elements to E
622*
623 DO 220 i = 1, n - 1
624 e( i ) = ab( 2, i )
625 220 CONTINUE
626 ELSE
627*
628* set E to zero if original matrix was diagonal
629*
630 DO 230 i = 1, n - 1
631 e( i ) = zero
632 230 CONTINUE
633 END IF
634*
635* copy diagonal elements to D
636*
637 DO 240 i = 1, n
638 d( i ) = ab( 1, i )
639 240 CONTINUE
640 END IF
641*
642 RETURN
643*
644* End of DSBTRD
645*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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:108
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:102
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:106
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:108
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: