OpenRadioss 2025.1.11
OpenRadioss project
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*> \htmlonly
9*> Download SSBTRD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssbtrd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssbtrd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssbtrd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SSBTRD( 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* REAL AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
30* $ WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SSBTRD 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 REAL 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 REAL array, dimension (N)
101*> The diagonal elements of the tridiagonal matrix T.
102*> \endverbatim
103*>
104*> \param[out] E
105*> \verbatim
106*> E is REAL 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 REAL 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 REAL 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*> \ingroup realOTHERcomputational
151*
152*> \par Further Details:
153* =====================
154*>
155*> \verbatim
156*>
157*> Modified by Linda Kaufman, Bell Labs.
158*> \endverbatim
159*>
160* =====================================================================
161 SUBROUTINE ssbtrd( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
162 $ WORK, INFO )
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 REAL AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
174 $ work( * )
175* ..
176*
177* =====================================================================
178*
179* .. Parameters ..
180 REAL ZERO, ONE
181 parameter( zero = 0.0e+0, one = 1.0e+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 REAL TEMP
189* ..
190* .. External Subroutines ..
191 EXTERNAL slar2v, slargv, slartg, slartv, slaset, srot,
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.NOT..AND..NOT. IF( WANTQ LSAME( VECT, 'n' ) ) THEN
215 INFO = -1
216.NOT..AND..NOT. ELSE IF( UPPER LSAME( UPLO, 'l' ) ) THEN
217 INFO = -2
218.LT. ELSE IF( N0 ) THEN
219 INFO = -3
220.LT. ELSE IF( KD0 ) THEN
221 INFO = -4
222.LT. ELSE IF( LDABKD1 ) THEN
223 INFO = -6
224.LT..AND. ELSE IF( LDQMAX( 1, N ) WANTQ ) THEN
225 INFO = -10
226 END IF
227.NE. IF( INFO0 ) THEN
228 CALL XERBLA( 'ssbtrd', -INFO )
229 RETURN
230 END IF
231*
232* Quick return if possible
233*
234.EQ. IF( N0 )
235 $ RETURN
236*
237* Initialize Q to the unit matrix, if needed
238*
239 IF( INITQ )
240 $ CALL SLASET( '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.GT. IF( KD1 ) 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.GT. IF( NR0 ) THEN
269*
270* generate plane rotations to annihilate nonzero
271* elements which have been created outside the band
272*
273 CALL SLARGV( 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* SLARTV or SROT is used
281*
282.GE. IF( NR2*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.GT. IF( K2 ) THEN
301.LE. IF( KN-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.GT. IF( NR0 )
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.GT. IF( NR0 ) THEN
332.LT. IF( 2*KD-1NR ) 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.GT. IF( J2+LN ) THEN
339 NRT = NR - 1
340 ELSE
341 NRT = NR
342 END IF
343.GT. IF( NRT0 )
344 $ CALL SLARTV( 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.GE. IF( J1ENDJ1 ) THEN
351 DO 40 JIN = J1, J1END, KD1
352 CALL SROT( 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.GT. IF( LEND0 )
360 $ CALL SROT( 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.EQ. IF( K2 )
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 SROT( 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 SROT( 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.GT. IF( J2+KDNN ) 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.GT. IF( KD0 ) 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.GT. IF( KD1 ) 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.GT. IF( NR0 ) THEN
461*
462* generate plane rotations to annihilate nonzero
463* elements which have been created outside the band
464*
465 CALL SLARGV( 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* SLARTV or SROT is used
473*
474.GT. IF( NR2*KD-1 ) THEN
475 DO 130 L = 1, KD - 1
476 CALL SLARTV( 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 SROT( 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.GT. IF( K2 ) THEN
492.LE. IF( KN-I+1 ) THEN
493*
494* generate plane rotation to annihilate a(i+k-1,i)
495* within the band
496*
497 CALL SLARTG( 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 SROT( 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.GT. IF( NR0 )
515 $ CALL SLAR2V( 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* SLARTV or SROT is used
524*
525.GT. IF( NR0 ) THEN
526.GT. IF( NR2*KD-1 ) THEN
527 DO 150 L = 1, KD - 1
528.GT. IF( J2+LN ) THEN
529 NRT = NR - 1
530 ELSE
531 NRT = NR
532 END IF
533.GT. IF( NRT0 )
534 $ CALL SLARTV( 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.GE. IF( J1ENDJ1 ) THEN
541 DO 160 J1INC = J1, J1END, KD1
542 CALL SROT( 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.GT. IF( LEND0 )
550 $ CALL SROT( 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.EQ. IF( K2 )
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 SROT( 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 SROT( 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.GT. IF( J2+KDNN ) 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.GT. IF( KD0 ) 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 SSBTRD
637*
638 END
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:110
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:113
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
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: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:104
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:108
ssbtrd
subroutine ssbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
SSBTRD
Definition ssbtrd.f:163
srot
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21