LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ssbtrd ( character  VECT,
character  UPLO,
integer  N,
integer  KD,
real, dimension( ldab, * )  AB,
integer  LDAB,
real, dimension( * )  D,
real, dimension( * )  E,
real, dimension( ldq, * )  Q,
integer  LDQ,
real, dimension( * )  WORK,
integer  INFO 
)

SSBTRD

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

Purpose: 
 SSBTRD 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 REAL 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 REAL array, dimension (N)
          The diagonal elements of the tridiagonal matrix T.
[out]E
          E is REAL 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 REAL 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 REAL 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.
Date
November 2011
Further Details: 
  Modified by Linda Kaufman, Bell Labs.

Definition at line 165 of file ssbtrd.f.

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  REAL ab( ldab, * ), d( * ), e( * ), q( ldq, * ),
177  $ work( * )
178 * ..
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183  REAL zero, one
184  parameter ( zero = 0.0e+0, one = 1.0e+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  REAL temp
192 * ..
193 * .. External Subroutines ..
194  EXTERNAL slar2v, slargv, slartg, slartv, slaset, srot,
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( 'SSBTRD', -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 slaset( '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 slargv( 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 * SLARTV or SROT is used
284 *
285  IF( nr.GE.2*kd-1 ) THEN
286  DO 10 l = 1, kd - 1
287  CALL slartv( 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 srot( 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 slartg( 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 srot( 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 slar2v( 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 * SLARTV or SROT 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 slartv( 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 srot( 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 srot( 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 srot( 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 srot( 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 slargv( 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 * SLARTV or SROT is used
476 *
477  IF( nr.GT.2*kd-1 ) THEN
478  DO 130 l = 1, kd - 1
479  CALL slartv( 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 srot( 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 slartg( 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 srot( 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 slar2v( 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 * SLARTV or SROT 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 slartv( 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 srot( 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 srot( 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 srot( 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 srot( 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 SSBTRD
640 *
slar2v
subroutine slar2v(N, X, Y, Z, INCX, C, S, INCC)
SLAR2V applies a vector of plane rotations with real cosines and real sines from both sides to a sequ...
Definition: slar2v.f:112
slartv
subroutine slartv(N, X, INCX, Y, INCY, C, S, INCC)
SLARTV applies a vector of plane rotations with real cosines and real sines to the elements of a pair...
Definition: slartv.f:110
slargv
subroutine slargv(N, X, INCX, Y, INCY, C, INCC)
SLARGV generates a vector of plane rotations with real cosines and real sines.
Definition: slargv.f:106
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
srot
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:53
slartg
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
slaset
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: