LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssbtrd.f
Go to the documentation of this file.
1*> \brief \b SSBTRD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download SSBTRD + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssbtrd.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssbtrd.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssbtrd.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
20* WORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER UPLO, VECT
24* INTEGER INFO, KD, LDAB, LDQ, N
25* ..
26* .. Array Arguments ..
27* REAL AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
28* $ WORK( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SSBTRD reduces a real symmetric band matrix A to symmetric
38*> tridiagonal form T by an orthogonal similarity transformation:
39*> Q**T * A * Q = T.
40*> \endverbatim
41*
42* Arguments:
43* ==========
44*
45*> \param[in] VECT
46*> \verbatim
47*> VECT is CHARACTER*1
48*> = 'N': do not form Q;
49*> = 'V': form Q;
50*> = 'U': update a matrix X, by forming X*Q.
51*> \endverbatim
52*>
53*> \param[in] UPLO
54*> \verbatim
55*> UPLO is CHARACTER*1
56*> = 'U': Upper triangle of A is stored;
57*> = 'L': Lower triangle of A is stored.
58*> \endverbatim
59*>
60*> \param[in] N
61*> \verbatim
62*> N is INTEGER
63*> The order of the matrix A. N >= 0.
64*> \endverbatim
65*>
66*> \param[in] KD
67*> \verbatim
68*> KD is INTEGER
69*> The number of superdiagonals of the matrix A if UPLO = 'U',
70*> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
71*> \endverbatim
72*>
73*> \param[in,out] AB
74*> \verbatim
75*> AB is REAL array, dimension (LDAB,N)
76*> On entry, the upper or lower triangle of the symmetric band
77*> matrix A, stored in the first KD+1 rows of the array. The
78*> j-th column of A is stored in the j-th column of the array AB
79*> as follows:
80*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
81*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
82*> On exit, the diagonal elements of AB are overwritten by the
83*> diagonal elements of the tridiagonal matrix T; if KD > 0, the
84*> elements on the first superdiagonal (if UPLO = 'U') or the
85*> first subdiagonal (if UPLO = 'L') are overwritten by the
86*> off-diagonal elements of T; the rest of AB is overwritten by
87*> values generated during the reduction.
88*> \endverbatim
89*>
90*> \param[in] LDAB
91*> \verbatim
92*> LDAB is INTEGER
93*> The leading dimension of the array AB. LDAB >= KD+1.
94*> \endverbatim
95*>
96*> \param[out] D
97*> \verbatim
98*> D is REAL array, dimension (N)
99*> The diagonal elements of the tridiagonal matrix T.
100*> \endverbatim
101*>
102*> \param[out] E
103*> \verbatim
104*> E is REAL array, dimension (N-1)
105*> The off-diagonal elements of the tridiagonal matrix T:
106*> E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
107*> \endverbatim
108*>
109*> \param[in,out] Q
110*> \verbatim
111*> Q is REAL array, dimension (LDQ,N)
112*> On entry, if VECT = 'U', then Q must contain an N-by-N
113*> matrix X; if VECT = 'N' or 'V', then Q need not be set.
114*>
115*> On exit:
116*> if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
117*> if VECT = 'U', Q contains the product X*Q;
118*> if VECT = 'N', the array Q is not referenced.
119*> \endverbatim
120*>
121*> \param[in] LDQ
122*> \verbatim
123*> LDQ is INTEGER
124*> The leading dimension of the array Q.
125*> LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
126*> \endverbatim
127*>
128*> \param[out] WORK
129*> \verbatim
130*> WORK is REAL array, dimension (N)
131*> \endverbatim
132*>
133*> \param[out] INFO
134*> \verbatim
135*> INFO is INTEGER
136*> = 0: successful exit
137*> < 0: if INFO = -i, the i-th argument had an illegal value
138*> \endverbatim
139*
140* Authors:
141* ========
142*
143*> \author Univ. of Tennessee
144*> \author Univ. of California Berkeley
145*> \author Univ. of Colorado Denver
146*> \author NAG Ltd.
147*
148*> \ingroup hbtrd
149*
150*> \par Further Details:
151* =====================
152*>
153*> \verbatim
154*>
155*> Modified by Linda Kaufman, Bell Labs.
156*> \endverbatim
157*>
158* =====================================================================
159 SUBROUTINE ssbtrd( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
160 $ WORK, INFO )
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 REAL AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
172 $ work( * )
173* ..
174*
175* =====================================================================
176*
177* .. Parameters ..
178 REAL ZERO, ONE
179 parameter( zero = 0.0e+0, one = 1.0e+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 REAL TEMP
187* ..
188* .. External Subroutines ..
189 EXTERNAL slar2v, slargv, slartg, slartv, slaset,
190 $ srot,
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( 'SSBTRD', -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 slaset( '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 slargv( 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* SLARTV or SROT is used
281*
282 IF( nr.GE.2*kd-1 ) THEN
283 DO 10 l = 1, kd - 1
284 CALL slartv( 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 srot( 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 slartg( 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 srot( 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 slar2v( 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* SLARTV or SROT 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 slartv( 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 srot( 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 srot( 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 srot( 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 srot( 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 slargv( 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* SLARTV or SROT is used
477*
478 IF( nr.GT.2*kd-1 ) THEN
479 DO 130 l = 1, kd - 1
480 CALL slartv( 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 srot( 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 slartg( 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 srot( 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 slar2v( 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* SLARTV or SROT 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 slartv( 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 srot( 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 srot( 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 srot( 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 srot( 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 SSBTRD
645*
646 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
ssbtrd
subroutine ssbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
SSBTRD
Definition ssbtrd.f:161
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:108
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:102
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
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:106
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:108
srot
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92