1 SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
2 $ WORK, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER UPLO, VECT
11 INTEGER INFO, KD, LDAB, LDQ, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
15 $ WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSBTRD reduces a real symmetric band matrix A to symmetric
22 * tridiagonal form T by an orthogonal similarity transformation:
23 * Q**T * A * Q = T.
24 *
25 * Arguments
26 * =========
27 *
28 * VECT (input) CHARACTER*1
29 * = 'N': do not form Q;
30 * = 'V': form Q;
31 * = 'U': update a matrix X, by forming X*Q.
32 *
33 * UPLO (input) CHARACTER*1
34 * = 'U': Upper triangle of A is stored;
35 * = 'L': Lower triangle of A is stored.
36 *
37 * N (input) INTEGER
38 * The order of the matrix A. N >= 0.
39 *
40 * KD (input) INTEGER
41 * The number of superdiagonals of the matrix A if UPLO = 'U',
42 * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
43 *
44 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
45 * On entry, the upper or lower triangle of the symmetric band
46 * matrix A, stored in the first KD+1 rows of the array. The
47 * j-th column of A is stored in the j-th column of the array AB
48 * as follows:
49 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
50 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
51 * On exit, the diagonal elements of AB are overwritten by the
52 * diagonal elements of the tridiagonal matrix T; if KD > 0, the
53 * elements on the first superdiagonal (if UPLO = 'U') or the
54 * first subdiagonal (if UPLO = 'L') are overwritten by the
55 * off-diagonal elements of T; the rest of AB is overwritten by
56 * values generated during the reduction.
57 *
58 * LDAB (input) INTEGER
59 * The leading dimension of the array AB. LDAB >= KD+1.
60 *
61 * D (output) DOUBLE PRECISION array, dimension (N)
62 * The diagonal elements of the tridiagonal matrix T.
63 *
64 * E (output) DOUBLE PRECISION array, dimension (N-1)
65 * The off-diagonal elements of the tridiagonal matrix T:
66 * E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
67 *
68 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
69 * On entry, if VECT = 'U', then Q must contain an N-by-N
70 * matrix X; if VECT = 'N' or 'V', then Q need not be set.
71 *
72 * On exit:
73 * if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
74 * if VECT = 'U', Q contains the product X*Q;
75 * if VECT = 'N', the array Q is not referenced.
76 *
77 * LDQ (input) INTEGER
78 * The leading dimension of the array Q.
79 * LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
80 *
81 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
82 *
83 * INFO (output) INTEGER
84 * = 0: successful exit
85 * < 0: if INFO = -i, the i-th argument had an illegal value
86 *
87 * Further Details
88 * ===============
89 *
90 * Modified by Linda Kaufman, Bell Labs.
91 *
92 * =====================================================================
93 *
94 * .. Parameters ..
95 DOUBLE PRECISION ZERO, ONE
96 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
97 * ..
98 * .. Local Scalars ..
99 LOGICAL INITQ, UPPER, WANTQ
100 INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
101 $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
102 $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
103 DOUBLE PRECISION TEMP
104 * ..
105 * .. External Subroutines ..
106 EXTERNAL DLAR2V, DLARGV, DLARTG, DLARTV, DLASET, DROT,
107 $ XERBLA
108 * ..
109 * .. Intrinsic Functions ..
110 INTRINSIC MAX, MIN
111 * ..
112 * .. External Functions ..
113 LOGICAL LSAME
114 EXTERNAL LSAME
115 * ..
116 * .. Executable Statements ..
117 *
118 * Test the input parameters
119 *
120 INITQ = LSAME( VECT, 'V' )
121 WANTQ = INITQ .OR. LSAME( VECT, 'U' )
122 UPPER = LSAME( UPLO, 'U' )
123 KD1 = KD + 1
124 KDM1 = KD - 1
125 INCX = LDAB - 1
126 IQEND = 1
127 *
128 INFO = 0
129 IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
130 INFO = -1
131 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
132 INFO = -2
133 ELSE IF( N.LT.0 ) THEN
134 INFO = -3
135 ELSE IF( KD.LT.0 ) THEN
136 INFO = -4
137 ELSE IF( LDAB.LT.KD1 ) THEN
138 INFO = -6
139 ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
140 INFO = -10
141 END IF
142 IF( INFO.NE.0 ) THEN
143 CALL XERBLA( 'DSBTRD', -INFO )
144 RETURN
145 END IF
146 *
147 * Quick return if possible
148 *
149 IF( N.EQ.0 )
150 $ RETURN
151 *
152 * Initialize Q to the unit matrix, if needed
153 *
154 IF( INITQ )
155 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
156 *
157 * Wherever possible, plane rotations are generated and applied in
158 * vector operations of length NR over the index set J1:J2:KD1.
159 *
160 * The cosines and sines of the plane rotations are stored in the
161 * arrays D and WORK.
162 *
163 INCA = KD1*LDAB
164 KDN = MIN( N-1, KD )
165 IF( UPPER ) THEN
166 *
167 IF( KD.GT.1 ) THEN
168 *
169 * Reduce to tridiagonal form, working with upper triangle
170 *
171 NR = 0
172 J1 = KDN + 2
173 J2 = 1
174 *
175 DO 90 I = 1, N - 2
176 *
177 * Reduce i-th row of matrix to tridiagonal form
178 *
179 DO 80 K = KDN + 1, 2, -1
180 J1 = J1 + KDN
181 J2 = J2 + KDN
182 *
183 IF( NR.GT.0 ) THEN
184 *
185 * generate plane rotations to annihilate nonzero
186 * elements which have been created outside the band
187 *
188 CALL DLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
189 $ KD1, D( J1 ), KD1 )
190 *
191 * apply rotations from the right
192 *
193 *
194 * Dependent on the the number of diagonals either
195 * DLARTV or DROT is used
196 *
197 IF( NR.GE.2*KD-1 ) THEN
198 DO 10 L = 1, KD - 1
199 CALL DLARTV( NR, AB( L+1, J1-1 ), INCA,
200 $ AB( L, J1 ), INCA, D( J1 ),
201 $ WORK( J1 ), KD1 )
202 10 CONTINUE
203 *
204 ELSE
205 JEND = J1 + ( NR-1 )*KD1
206 DO 20 JINC = J1, JEND, KD1
207 CALL DROT( KDM1, AB( 2, JINC-1 ), 1,
208 $ AB( 1, JINC ), 1, D( JINC ),
209 $ WORK( JINC ) )
210 20 CONTINUE
211 END IF
212 END IF
213 *
214 *
215 IF( K.GT.2 ) THEN
216 IF( K.LE.N-I+1 ) THEN
217 *
218 * generate plane rotation to annihilate a(i,i+k-1)
219 * within the band
220 *
221 CALL DLARTG( AB( KD-K+3, I+K-2 ),
222 $ AB( KD-K+2, I+K-1 ), D( I+K-1 ),
223 $ WORK( I+K-1 ), TEMP )
224 AB( KD-K+3, I+K-2 ) = TEMP
225 *
226 * apply rotation from the right
227 *
228 CALL DROT( K-3, AB( KD-K+4, I+K-2 ), 1,
229 $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
230 $ WORK( I+K-1 ) )
231 END IF
232 NR = NR + 1
233 J1 = J1 - KDN - 1
234 END IF
235 *
236 * apply plane rotations from both sides to diagonal
237 * blocks
238 *
239 IF( NR.GT.0 )
240 $ CALL DLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
241 $ AB( KD, J1 ), INCA, D( J1 ),
242 $ WORK( J1 ), KD1 )
243 *
244 * apply plane rotations from the left
245 *
246 IF( NR.GT.0 ) THEN
247 IF( 2*KD-1.LT.NR ) THEN
248 *
249 * Dependent on the the number of diagonals either
250 * DLARTV or DROT is used
251 *
252 DO 30 L = 1, KD - 1
253 IF( J2+L.GT.N ) THEN
254 NRT = NR - 1
255 ELSE
256 NRT = NR
257 END IF
258 IF( NRT.GT.0 )
259 $ CALL DLARTV( NRT, AB( KD-L, J1+L ), INCA,
260 $ AB( KD-L+1, J1+L ), INCA,
261 $ D( J1 ), WORK( J1 ), KD1 )
262 30 CONTINUE
263 ELSE
264 J1END = J1 + KD1*( NR-2 )
265 IF( J1END.GE.J1 ) THEN
266 DO 40 JIN = J1, J1END, KD1
267 CALL DROT( KD-1, AB( KD-1, JIN+1 ), INCX,
268 $ AB( KD, JIN+1 ), INCX,
269 $ D( JIN ), WORK( JIN ) )
270 40 CONTINUE
271 END IF
272 LEND = MIN( KDM1, N-J2 )
273 LAST = J1END + KD1
274 IF( LEND.GT.0 )
275 $ CALL DROT( LEND, AB( KD-1, LAST+1 ), INCX,
276 $ AB( KD, LAST+1 ), INCX, D( LAST ),
277 $ WORK( LAST ) )
278 END IF
279 END IF
280 *
281 IF( WANTQ ) THEN
282 *
283 * accumulate product of plane rotations in Q
284 *
285 IF( INITQ ) THEN
286 *
287 * take advantage of the fact that Q was
288 * initially the Identity matrix
289 *
290 IQEND = MAX( IQEND, J2 )
291 I2 = MAX( 0, K-3 )
292 IQAEND = 1 + I*KD
293 IF( K.EQ.2 )
294 $ IQAEND = IQAEND + KD
295 IQAEND = MIN( IQAEND, IQEND )
296 DO 50 J = J1, J2, KD1
297 IBL = I - I2 / KDM1
298 I2 = I2 + 1
299 IQB = MAX( 1, J-IBL )
300 NQ = 1 + IQAEND - IQB
301 IQAEND = MIN( IQAEND+KD, IQEND )
302 CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
303 $ 1, D( J ), WORK( J ) )
304 50 CONTINUE
305 ELSE
306 *
307 DO 60 J = J1, J2, KD1
308 CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
309 $ D( J ), WORK( J ) )
310 60 CONTINUE
311 END IF
312 *
313 END IF
314 *
315 IF( J2+KDN.GT.N ) THEN
316 *
317 * adjust J2 to keep within the bounds of the matrix
318 *
319 NR = NR - 1
320 J2 = J2 - KDN - 1
321 END IF
322 *
323 DO 70 J = J1, J2, KD1
324 *
325 * create nonzero element a(j-1,j+kd) outside the band
326 * and store it in WORK
327 *
328 WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
329 AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
330 70 CONTINUE
331 80 CONTINUE
332 90 CONTINUE
333 END IF
334 *
335 IF( KD.GT.0 ) THEN
336 *
337 * copy off-diagonal elements to E
338 *
339 DO 100 I = 1, N - 1
340 E( I ) = AB( KD, I+1 )
341 100 CONTINUE
342 ELSE
343 *
344 * set E to zero if original matrix was diagonal
345 *
346 DO 110 I = 1, N - 1
347 E( I ) = ZERO
348 110 CONTINUE
349 END IF
350 *
351 * copy diagonal elements to D
352 *
353 DO 120 I = 1, N
354 D( I ) = AB( KD1, I )
355 120 CONTINUE
356 *
357 ELSE
358 *
359 IF( KD.GT.1 ) THEN
360 *
361 * Reduce to tridiagonal form, working with lower triangle
362 *
363 NR = 0
364 J1 = KDN + 2
365 J2 = 1
366 *
367 DO 210 I = 1, N - 2
368 *
369 * Reduce i-th column of matrix to tridiagonal form
370 *
371 DO 200 K = KDN + 1, 2, -1
372 J1 = J1 + KDN
373 J2 = J2 + KDN
374 *
375 IF( NR.GT.0 ) THEN
376 *
377 * generate plane rotations to annihilate nonzero
378 * elements which have been created outside the band
379 *
380 CALL DLARGV( NR, AB( KD1, J1-KD1 ), INCA,
381 $ WORK( J1 ), KD1, D( J1 ), KD1 )
382 *
383 * apply plane rotations from one side
384 *
385 *
386 * Dependent on the the number of diagonals either
387 * DLARTV or DROT is used
388 *
389 IF( NR.GT.2*KD-1 ) THEN
390 DO 130 L = 1, KD - 1
391 CALL DLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
392 $ AB( KD1-L+1, J1-KD1+L ), INCA,
393 $ D( J1 ), WORK( J1 ), KD1 )
394 130 CONTINUE
395 ELSE
396 JEND = J1 + KD1*( NR-1 )
397 DO 140 JINC = J1, JEND, KD1
398 CALL DROT( KDM1, AB( KD, JINC-KD ), INCX,
399 $ AB( KD1, JINC-KD ), INCX,
400 $ D( JINC ), WORK( JINC ) )
401 140 CONTINUE
402 END IF
403 *
404 END IF
405 *
406 IF( K.GT.2 ) THEN
407 IF( K.LE.N-I+1 ) THEN
408 *
409 * generate plane rotation to annihilate a(i+k-1,i)
410 * within the band
411 *
412 CALL DLARTG( AB( K-1, I ), AB( K, I ),
413 $ D( I+K-1 ), WORK( I+K-1 ), TEMP )
414 AB( K-1, I ) = TEMP
415 *
416 * apply rotation from the left
417 *
418 CALL DROT( K-3, AB( K-2, I+1 ), LDAB-1,
419 $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
420 $ WORK( I+K-1 ) )
421 END IF
422 NR = NR + 1
423 J1 = J1 - KDN - 1
424 END IF
425 *
426 * apply plane rotations from both sides to diagonal
427 * blocks
428 *
429 IF( NR.GT.0 )
430 $ CALL DLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
431 $ AB( 2, J1-1 ), INCA, D( J1 ),
432 $ WORK( J1 ), KD1 )
433 *
434 * apply plane rotations from the right
435 *
436 *
437 * Dependent on the the number of diagonals either
438 * DLARTV or DROT is used
439 *
440 IF( NR.GT.0 ) THEN
441 IF( NR.GT.2*KD-1 ) THEN
442 DO 150 L = 1, KD - 1
443 IF( J2+L.GT.N ) THEN
444 NRT = NR - 1
445 ELSE
446 NRT = NR
447 END IF
448 IF( NRT.GT.0 )
449 $ CALL DLARTV( NRT, AB( L+2, J1-1 ), INCA,
450 $ AB( L+1, J1 ), INCA, D( J1 ),
451 $ WORK( J1 ), KD1 )
452 150 CONTINUE
453 ELSE
454 J1END = J1 + KD1*( NR-2 )
455 IF( J1END.GE.J1 ) THEN
456 DO 160 J1INC = J1, J1END, KD1
457 CALL DROT( KDM1, AB( 3, J1INC-1 ), 1,
458 $ AB( 2, J1INC ), 1, D( J1INC ),
459 $ WORK( J1INC ) )
460 160 CONTINUE
461 END IF
462 LEND = MIN( KDM1, N-J2 )
463 LAST = J1END + KD1
464 IF( LEND.GT.0 )
465 $ CALL DROT( LEND, AB( 3, LAST-1 ), 1,
466 $ AB( 2, LAST ), 1, D( LAST ),
467 $ WORK( LAST ) )
468 END IF
469 END IF
470 *
471 *
472 *
473 IF( WANTQ ) THEN
474 *
475 * accumulate product of plane rotations in Q
476 *
477 IF( INITQ ) THEN
478 *
479 * take advantage of the fact that Q was
480 * initially the Identity matrix
481 *
482 IQEND = MAX( IQEND, J2 )
483 I2 = MAX( 0, K-3 )
484 IQAEND = 1 + I*KD
485 IF( K.EQ.2 )
486 $ IQAEND = IQAEND + KD
487 IQAEND = MIN( IQAEND, IQEND )
488 DO 170 J = J1, J2, KD1
489 IBL = I - I2 / KDM1
490 I2 = I2 + 1
491 IQB = MAX( 1, J-IBL )
492 NQ = 1 + IQAEND - IQB
493 IQAEND = MIN( IQAEND+KD, IQEND )
494 CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
495 $ 1, D( J ), WORK( J ) )
496 170 CONTINUE
497 ELSE
498 *
499 DO 180 J = J1, J2, KD1
500 CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
501 $ D( J ), WORK( J ) )
502 180 CONTINUE
503 END IF
504 END IF
505 *
506 IF( J2+KDN.GT.N ) THEN
507 *
508 * adjust J2 to keep within the bounds of the matrix
509 *
510 NR = NR - 1
511 J2 = J2 - KDN - 1
512 END IF
513 *
514 DO 190 J = J1, J2, KD1
515 *
516 * create nonzero element a(j+kd,j-1) outside the
517 * band and store it in WORK
518 *
519 WORK( J+KD ) = WORK( J )*AB( KD1, J )
520 AB( KD1, J ) = D( J )*AB( KD1, J )
521 190 CONTINUE
522 200 CONTINUE
523 210 CONTINUE
524 END IF
525 *
526 IF( KD.GT.0 ) THEN
527 *
528 * copy off-diagonal elements to E
529 *
530 DO 220 I = 1, N - 1
531 E( I ) = AB( 2, I )
532 220 CONTINUE
533 ELSE
534 *
535 * set E to zero if original matrix was diagonal
536 *
537 DO 230 I = 1, N - 1
538 E( I ) = ZERO
539 230 CONTINUE
540 END IF
541 *
542 * copy diagonal elements to D
543 *
544 DO 240 I = 1, N
545 D( I ) = AB( 1, I )
546 240 CONTINUE
547 END IF
548 *
549 RETURN
550 *
551 * End of DSBTRD
552 *
553 END
2 $ WORK, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER UPLO, VECT
11 INTEGER INFO, KD, LDAB, LDQ, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
15 $ WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DSBTRD reduces a real symmetric band matrix A to symmetric
22 * tridiagonal form T by an orthogonal similarity transformation:
23 * Q**T * A * Q = T.
24 *
25 * Arguments
26 * =========
27 *
28 * VECT (input) CHARACTER*1
29 * = 'N': do not form Q;
30 * = 'V': form Q;
31 * = 'U': update a matrix X, by forming X*Q.
32 *
33 * UPLO (input) CHARACTER*1
34 * = 'U': Upper triangle of A is stored;
35 * = 'L': Lower triangle of A is stored.
36 *
37 * N (input) INTEGER
38 * The order of the matrix A. N >= 0.
39 *
40 * KD (input) INTEGER
41 * The number of superdiagonals of the matrix A if UPLO = 'U',
42 * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
43 *
44 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
45 * On entry, the upper or lower triangle of the symmetric band
46 * matrix A, stored in the first KD+1 rows of the array. The
47 * j-th column of A is stored in the j-th column of the array AB
48 * as follows:
49 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
50 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
51 * On exit, the diagonal elements of AB are overwritten by the
52 * diagonal elements of the tridiagonal matrix T; if KD > 0, the
53 * elements on the first superdiagonal (if UPLO = 'U') or the
54 * first subdiagonal (if UPLO = 'L') are overwritten by the
55 * off-diagonal elements of T; the rest of AB is overwritten by
56 * values generated during the reduction.
57 *
58 * LDAB (input) INTEGER
59 * The leading dimension of the array AB. LDAB >= KD+1.
60 *
61 * D (output) DOUBLE PRECISION array, dimension (N)
62 * The diagonal elements of the tridiagonal matrix T.
63 *
64 * E (output) DOUBLE PRECISION array, dimension (N-1)
65 * The off-diagonal elements of the tridiagonal matrix T:
66 * E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
67 *
68 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
69 * On entry, if VECT = 'U', then Q must contain an N-by-N
70 * matrix X; if VECT = 'N' or 'V', then Q need not be set.
71 *
72 * On exit:
73 * if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
74 * if VECT = 'U', Q contains the product X*Q;
75 * if VECT = 'N', the array Q is not referenced.
76 *
77 * LDQ (input) INTEGER
78 * The leading dimension of the array Q.
79 * LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
80 *
81 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
82 *
83 * INFO (output) INTEGER
84 * = 0: successful exit
85 * < 0: if INFO = -i, the i-th argument had an illegal value
86 *
87 * Further Details
88 * ===============
89 *
90 * Modified by Linda Kaufman, Bell Labs.
91 *
92 * =====================================================================
93 *
94 * .. Parameters ..
95 DOUBLE PRECISION ZERO, ONE
96 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
97 * ..
98 * .. Local Scalars ..
99 LOGICAL INITQ, UPPER, WANTQ
100 INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
101 $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
102 $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
103 DOUBLE PRECISION TEMP
104 * ..
105 * .. External Subroutines ..
106 EXTERNAL DLAR2V, DLARGV, DLARTG, DLARTV, DLASET, DROT,
107 $ XERBLA
108 * ..
109 * .. Intrinsic Functions ..
110 INTRINSIC MAX, MIN
111 * ..
112 * .. External Functions ..
113 LOGICAL LSAME
114 EXTERNAL LSAME
115 * ..
116 * .. Executable Statements ..
117 *
118 * Test the input parameters
119 *
120 INITQ = LSAME( VECT, 'V' )
121 WANTQ = INITQ .OR. LSAME( VECT, 'U' )
122 UPPER = LSAME( UPLO, 'U' )
123 KD1 = KD + 1
124 KDM1 = KD - 1
125 INCX = LDAB - 1
126 IQEND = 1
127 *
128 INFO = 0
129 IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
130 INFO = -1
131 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
132 INFO = -2
133 ELSE IF( N.LT.0 ) THEN
134 INFO = -3
135 ELSE IF( KD.LT.0 ) THEN
136 INFO = -4
137 ELSE IF( LDAB.LT.KD1 ) THEN
138 INFO = -6
139 ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
140 INFO = -10
141 END IF
142 IF( INFO.NE.0 ) THEN
143 CALL XERBLA( 'DSBTRD', -INFO )
144 RETURN
145 END IF
146 *
147 * Quick return if possible
148 *
149 IF( N.EQ.0 )
150 $ RETURN
151 *
152 * Initialize Q to the unit matrix, if needed
153 *
154 IF( INITQ )
155 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
156 *
157 * Wherever possible, plane rotations are generated and applied in
158 * vector operations of length NR over the index set J1:J2:KD1.
159 *
160 * The cosines and sines of the plane rotations are stored in the
161 * arrays D and WORK.
162 *
163 INCA = KD1*LDAB
164 KDN = MIN( N-1, KD )
165 IF( UPPER ) THEN
166 *
167 IF( KD.GT.1 ) THEN
168 *
169 * Reduce to tridiagonal form, working with upper triangle
170 *
171 NR = 0
172 J1 = KDN + 2
173 J2 = 1
174 *
175 DO 90 I = 1, N - 2
176 *
177 * Reduce i-th row of matrix to tridiagonal form
178 *
179 DO 80 K = KDN + 1, 2, -1
180 J1 = J1 + KDN
181 J2 = J2 + KDN
182 *
183 IF( NR.GT.0 ) THEN
184 *
185 * generate plane rotations to annihilate nonzero
186 * elements which have been created outside the band
187 *
188 CALL DLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
189 $ KD1, D( J1 ), KD1 )
190 *
191 * apply rotations from the right
192 *
193 *
194 * Dependent on the the number of diagonals either
195 * DLARTV or DROT is used
196 *
197 IF( NR.GE.2*KD-1 ) THEN
198 DO 10 L = 1, KD - 1
199 CALL DLARTV( NR, AB( L+1, J1-1 ), INCA,
200 $ AB( L, J1 ), INCA, D( J1 ),
201 $ WORK( J1 ), KD1 )
202 10 CONTINUE
203 *
204 ELSE
205 JEND = J1 + ( NR-1 )*KD1
206 DO 20 JINC = J1, JEND, KD1
207 CALL DROT( KDM1, AB( 2, JINC-1 ), 1,
208 $ AB( 1, JINC ), 1, D( JINC ),
209 $ WORK( JINC ) )
210 20 CONTINUE
211 END IF
212 END IF
213 *
214 *
215 IF( K.GT.2 ) THEN
216 IF( K.LE.N-I+1 ) THEN
217 *
218 * generate plane rotation to annihilate a(i,i+k-1)
219 * within the band
220 *
221 CALL DLARTG( AB( KD-K+3, I+K-2 ),
222 $ AB( KD-K+2, I+K-1 ), D( I+K-1 ),
223 $ WORK( I+K-1 ), TEMP )
224 AB( KD-K+3, I+K-2 ) = TEMP
225 *
226 * apply rotation from the right
227 *
228 CALL DROT( K-3, AB( KD-K+4, I+K-2 ), 1,
229 $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
230 $ WORK( I+K-1 ) )
231 END IF
232 NR = NR + 1
233 J1 = J1 - KDN - 1
234 END IF
235 *
236 * apply plane rotations from both sides to diagonal
237 * blocks
238 *
239 IF( NR.GT.0 )
240 $ CALL DLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
241 $ AB( KD, J1 ), INCA, D( J1 ),
242 $ WORK( J1 ), KD1 )
243 *
244 * apply plane rotations from the left
245 *
246 IF( NR.GT.0 ) THEN
247 IF( 2*KD-1.LT.NR ) THEN
248 *
249 * Dependent on the the number of diagonals either
250 * DLARTV or DROT is used
251 *
252 DO 30 L = 1, KD - 1
253 IF( J2+L.GT.N ) THEN
254 NRT = NR - 1
255 ELSE
256 NRT = NR
257 END IF
258 IF( NRT.GT.0 )
259 $ CALL DLARTV( NRT, AB( KD-L, J1+L ), INCA,
260 $ AB( KD-L+1, J1+L ), INCA,
261 $ D( J1 ), WORK( J1 ), KD1 )
262 30 CONTINUE
263 ELSE
264 J1END = J1 + KD1*( NR-2 )
265 IF( J1END.GE.J1 ) THEN
266 DO 40 JIN = J1, J1END, KD1
267 CALL DROT( KD-1, AB( KD-1, JIN+1 ), INCX,
268 $ AB( KD, JIN+1 ), INCX,
269 $ D( JIN ), WORK( JIN ) )
270 40 CONTINUE
271 END IF
272 LEND = MIN( KDM1, N-J2 )
273 LAST = J1END + KD1
274 IF( LEND.GT.0 )
275 $ CALL DROT( LEND, AB( KD-1, LAST+1 ), INCX,
276 $ AB( KD, LAST+1 ), INCX, D( LAST ),
277 $ WORK( LAST ) )
278 END IF
279 END IF
280 *
281 IF( WANTQ ) THEN
282 *
283 * accumulate product of plane rotations in Q
284 *
285 IF( INITQ ) THEN
286 *
287 * take advantage of the fact that Q was
288 * initially the Identity matrix
289 *
290 IQEND = MAX( IQEND, J2 )
291 I2 = MAX( 0, K-3 )
292 IQAEND = 1 + I*KD
293 IF( K.EQ.2 )
294 $ IQAEND = IQAEND + KD
295 IQAEND = MIN( IQAEND, IQEND )
296 DO 50 J = J1, J2, KD1
297 IBL = I - I2 / KDM1
298 I2 = I2 + 1
299 IQB = MAX( 1, J-IBL )
300 NQ = 1 + IQAEND - IQB
301 IQAEND = MIN( IQAEND+KD, IQEND )
302 CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
303 $ 1, D( J ), WORK( J ) )
304 50 CONTINUE
305 ELSE
306 *
307 DO 60 J = J1, J2, KD1
308 CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
309 $ D( J ), WORK( J ) )
310 60 CONTINUE
311 END IF
312 *
313 END IF
314 *
315 IF( J2+KDN.GT.N ) THEN
316 *
317 * adjust J2 to keep within the bounds of the matrix
318 *
319 NR = NR - 1
320 J2 = J2 - KDN - 1
321 END IF
322 *
323 DO 70 J = J1, J2, KD1
324 *
325 * create nonzero element a(j-1,j+kd) outside the band
326 * and store it in WORK
327 *
328 WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
329 AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
330 70 CONTINUE
331 80 CONTINUE
332 90 CONTINUE
333 END IF
334 *
335 IF( KD.GT.0 ) THEN
336 *
337 * copy off-diagonal elements to E
338 *
339 DO 100 I = 1, N - 1
340 E( I ) = AB( KD, I+1 )
341 100 CONTINUE
342 ELSE
343 *
344 * set E to zero if original matrix was diagonal
345 *
346 DO 110 I = 1, N - 1
347 E( I ) = ZERO
348 110 CONTINUE
349 END IF
350 *
351 * copy diagonal elements to D
352 *
353 DO 120 I = 1, N
354 D( I ) = AB( KD1, I )
355 120 CONTINUE
356 *
357 ELSE
358 *
359 IF( KD.GT.1 ) THEN
360 *
361 * Reduce to tridiagonal form, working with lower triangle
362 *
363 NR = 0
364 J1 = KDN + 2
365 J2 = 1
366 *
367 DO 210 I = 1, N - 2
368 *
369 * Reduce i-th column of matrix to tridiagonal form
370 *
371 DO 200 K = KDN + 1, 2, -1
372 J1 = J1 + KDN
373 J2 = J2 + KDN
374 *
375 IF( NR.GT.0 ) THEN
376 *
377 * generate plane rotations to annihilate nonzero
378 * elements which have been created outside the band
379 *
380 CALL DLARGV( NR, AB( KD1, J1-KD1 ), INCA,
381 $ WORK( J1 ), KD1, D( J1 ), KD1 )
382 *
383 * apply plane rotations from one side
384 *
385 *
386 * Dependent on the the number of diagonals either
387 * DLARTV or DROT is used
388 *
389 IF( NR.GT.2*KD-1 ) THEN
390 DO 130 L = 1, KD - 1
391 CALL DLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
392 $ AB( KD1-L+1, J1-KD1+L ), INCA,
393 $ D( J1 ), WORK( J1 ), KD1 )
394 130 CONTINUE
395 ELSE
396 JEND = J1 + KD1*( NR-1 )
397 DO 140 JINC = J1, JEND, KD1
398 CALL DROT( KDM1, AB( KD, JINC-KD ), INCX,
399 $ AB( KD1, JINC-KD ), INCX,
400 $ D( JINC ), WORK( JINC ) )
401 140 CONTINUE
402 END IF
403 *
404 END IF
405 *
406 IF( K.GT.2 ) THEN
407 IF( K.LE.N-I+1 ) THEN
408 *
409 * generate plane rotation to annihilate a(i+k-1,i)
410 * within the band
411 *
412 CALL DLARTG( AB( K-1, I ), AB( K, I ),
413 $ D( I+K-1 ), WORK( I+K-1 ), TEMP )
414 AB( K-1, I ) = TEMP
415 *
416 * apply rotation from the left
417 *
418 CALL DROT( K-3, AB( K-2, I+1 ), LDAB-1,
419 $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
420 $ WORK( I+K-1 ) )
421 END IF
422 NR = NR + 1
423 J1 = J1 - KDN - 1
424 END IF
425 *
426 * apply plane rotations from both sides to diagonal
427 * blocks
428 *
429 IF( NR.GT.0 )
430 $ CALL DLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
431 $ AB( 2, J1-1 ), INCA, D( J1 ),
432 $ WORK( J1 ), KD1 )
433 *
434 * apply plane rotations from the right
435 *
436 *
437 * Dependent on the the number of diagonals either
438 * DLARTV or DROT is used
439 *
440 IF( NR.GT.0 ) THEN
441 IF( NR.GT.2*KD-1 ) THEN
442 DO 150 L = 1, KD - 1
443 IF( J2+L.GT.N ) THEN
444 NRT = NR - 1
445 ELSE
446 NRT = NR
447 END IF
448 IF( NRT.GT.0 )
449 $ CALL DLARTV( NRT, AB( L+2, J1-1 ), INCA,
450 $ AB( L+1, J1 ), INCA, D( J1 ),
451 $ WORK( J1 ), KD1 )
452 150 CONTINUE
453 ELSE
454 J1END = J1 + KD1*( NR-2 )
455 IF( J1END.GE.J1 ) THEN
456 DO 160 J1INC = J1, J1END, KD1
457 CALL DROT( KDM1, AB( 3, J1INC-1 ), 1,
458 $ AB( 2, J1INC ), 1, D( J1INC ),
459 $ WORK( J1INC ) )
460 160 CONTINUE
461 END IF
462 LEND = MIN( KDM1, N-J2 )
463 LAST = J1END + KD1
464 IF( LEND.GT.0 )
465 $ CALL DROT( LEND, AB( 3, LAST-1 ), 1,
466 $ AB( 2, LAST ), 1, D( LAST ),
467 $ WORK( LAST ) )
468 END IF
469 END IF
470 *
471 *
472 *
473 IF( WANTQ ) THEN
474 *
475 * accumulate product of plane rotations in Q
476 *
477 IF( INITQ ) THEN
478 *
479 * take advantage of the fact that Q was
480 * initially the Identity matrix
481 *
482 IQEND = MAX( IQEND, J2 )
483 I2 = MAX( 0, K-3 )
484 IQAEND = 1 + I*KD
485 IF( K.EQ.2 )
486 $ IQAEND = IQAEND + KD
487 IQAEND = MIN( IQAEND, IQEND )
488 DO 170 J = J1, J2, KD1
489 IBL = I - I2 / KDM1
490 I2 = I2 + 1
491 IQB = MAX( 1, J-IBL )
492 NQ = 1 + IQAEND - IQB
493 IQAEND = MIN( IQAEND+KD, IQEND )
494 CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
495 $ 1, D( J ), WORK( J ) )
496 170 CONTINUE
497 ELSE
498 *
499 DO 180 J = J1, J2, KD1
500 CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
501 $ D( J ), WORK( J ) )
502 180 CONTINUE
503 END IF
504 END IF
505 *
506 IF( J2+KDN.GT.N ) THEN
507 *
508 * adjust J2 to keep within the bounds of the matrix
509 *
510 NR = NR - 1
511 J2 = J2 - KDN - 1
512 END IF
513 *
514 DO 190 J = J1, J2, KD1
515 *
516 * create nonzero element a(j+kd,j-1) outside the
517 * band and store it in WORK
518 *
519 WORK( J+KD ) = WORK( J )*AB( KD1, J )
520 AB( KD1, J ) = D( J )*AB( KD1, J )
521 190 CONTINUE
522 200 CONTINUE
523 210 CONTINUE
524 END IF
525 *
526 IF( KD.GT.0 ) THEN
527 *
528 * copy off-diagonal elements to E
529 *
530 DO 220 I = 1, N - 1
531 E( I ) = AB( 2, I )
532 220 CONTINUE
533 ELSE
534 *
535 * set E to zero if original matrix was diagonal
536 *
537 DO 230 I = 1, N - 1
538 E( I ) = ZERO
539 230 CONTINUE
540 END IF
541 *
542 * copy diagonal elements to D
543 *
544 DO 240 I = 1, N
545 D( I ) = AB( 1, I )
546 240 CONTINUE
547 END IF
548 *
549 RETURN
550 *
551 * End of DSBTRD
552 *
553 END