|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, 00002 $ WORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER UPLO, VECT 00011 INTEGER INFO, KD, LDAB, LDQ, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ), 00015 $ WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DSBTRD reduces a real symmetric band matrix A to symmetric 00022 * tridiagonal form T by an orthogonal similarity transformation: 00023 * Q**T * A * Q = T. 00024 * 00025 * Arguments 00026 * ========= 00027 * 00028 * VECT (input) CHARACTER*1 00029 * = 'N': do not form Q; 00030 * = 'V': form Q; 00031 * = 'U': update a matrix X, by forming X*Q. 00032 * 00033 * UPLO (input) CHARACTER*1 00034 * = 'U': Upper triangle of A is stored; 00035 * = 'L': Lower triangle of A is stored. 00036 * 00037 * N (input) INTEGER 00038 * The order of the matrix A. N >= 0. 00039 * 00040 * KD (input) INTEGER 00041 * The number of superdiagonals of the matrix A if UPLO = 'U', 00042 * or the number of subdiagonals if UPLO = 'L'. KD >= 0. 00043 * 00044 * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) 00045 * On entry, the upper or lower triangle of the symmetric band 00046 * matrix A, stored in the first KD+1 rows of the array. The 00047 * j-th column of A is stored in the j-th column of the array AB 00048 * as follows: 00049 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 00050 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 00051 * On exit, the diagonal elements of AB are overwritten by the 00052 * diagonal elements of the tridiagonal matrix T; if KD > 0, the 00053 * elements on the first superdiagonal (if UPLO = 'U') or the 00054 * first subdiagonal (if UPLO = 'L') are overwritten by the 00055 * off-diagonal elements of T; the rest of AB is overwritten by 00056 * values generated during the reduction. 00057 * 00058 * LDAB (input) INTEGER 00059 * The leading dimension of the array AB. LDAB >= KD+1. 00060 * 00061 * D (output) DOUBLE PRECISION array, dimension (N) 00062 * The diagonal elements of the tridiagonal matrix T. 00063 * 00064 * E (output) DOUBLE PRECISION array, dimension (N-1) 00065 * The off-diagonal elements of the tridiagonal matrix T: 00066 * E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'. 00067 * 00068 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) 00069 * On entry, if VECT = 'U', then Q must contain an N-by-N 00070 * matrix X; if VECT = 'N' or 'V', then Q need not be set. 00071 * 00072 * On exit: 00073 * if VECT = 'V', Q contains the N-by-N orthogonal matrix Q; 00074 * if VECT = 'U', Q contains the product X*Q; 00075 * if VECT = 'N', the array Q is not referenced. 00076 * 00077 * LDQ (input) INTEGER 00078 * The leading dimension of the array Q. 00079 * LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'. 00080 * 00081 * WORK (workspace) DOUBLE PRECISION array, dimension (N) 00082 * 00083 * INFO (output) INTEGER 00084 * = 0: successful exit 00085 * < 0: if INFO = -i, the i-th argument had an illegal value 00086 * 00087 * Further Details 00088 * =============== 00089 * 00090 * Modified by Linda Kaufman, Bell Labs. 00091 * 00092 * ===================================================================== 00093 * 00094 * .. Parameters .. 00095 DOUBLE PRECISION ZERO, ONE 00096 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00097 * .. 00098 * .. Local Scalars .. 00099 LOGICAL INITQ, UPPER, WANTQ 00100 INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J, 00101 $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1, 00102 $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT 00103 DOUBLE PRECISION TEMP 00104 * .. 00105 * .. External Subroutines .. 00106 EXTERNAL DLAR2V, DLARGV, DLARTG, DLARTV, DLASET, DROT, 00107 $ XERBLA 00108 * .. 00109 * .. Intrinsic Functions .. 00110 INTRINSIC MAX, MIN 00111 * .. 00112 * .. External Functions .. 00113 LOGICAL LSAME 00114 EXTERNAL LSAME 00115 * .. 00116 * .. Executable Statements .. 00117 * 00118 * Test the input parameters 00119 * 00120 INITQ = LSAME( VECT, 'V' ) 00121 WANTQ = INITQ .OR. LSAME( VECT, 'U' ) 00122 UPPER = LSAME( UPLO, 'U' ) 00123 KD1 = KD + 1 00124 KDM1 = KD - 1 00125 INCX = LDAB - 1 00126 IQEND = 1 00127 * 00128 INFO = 0 00129 IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN 00130 INFO = -1 00131 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00132 INFO = -2 00133 ELSE IF( N.LT.0 ) THEN 00134 INFO = -3 00135 ELSE IF( KD.LT.0 ) THEN 00136 INFO = -4 00137 ELSE IF( LDAB.LT.KD1 ) THEN 00138 INFO = -6 00139 ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN 00140 INFO = -10 00141 END IF 00142 IF( INFO.NE.0 ) THEN 00143 CALL XERBLA( 'DSBTRD', -INFO ) 00144 RETURN 00145 END IF 00146 * 00147 * Quick return if possible 00148 * 00149 IF( N.EQ.0 ) 00150 $ RETURN 00151 * 00152 * Initialize Q to the unit matrix, if needed 00153 * 00154 IF( INITQ ) 00155 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) 00156 * 00157 * Wherever possible, plane rotations are generated and applied in 00158 * vector operations of length NR over the index set J1:J2:KD1. 00159 * 00160 * The cosines and sines of the plane rotations are stored in the 00161 * arrays D and WORK. 00162 * 00163 INCA = KD1*LDAB 00164 KDN = MIN( N-1, KD ) 00165 IF( UPPER ) THEN 00166 * 00167 IF( KD.GT.1 ) THEN 00168 * 00169 * Reduce to tridiagonal form, working with upper triangle 00170 * 00171 NR = 0 00172 J1 = KDN + 2 00173 J2 = 1 00174 * 00175 DO 90 I = 1, N - 2 00176 * 00177 * Reduce i-th row of matrix to tridiagonal form 00178 * 00179 DO 80 K = KDN + 1, 2, -1 00180 J1 = J1 + KDN 00181 J2 = J2 + KDN 00182 * 00183 IF( NR.GT.0 ) THEN 00184 * 00185 * generate plane rotations to annihilate nonzero 00186 * elements which have been created outside the band 00187 * 00188 CALL DLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ), 00189 $ KD1, D( J1 ), KD1 ) 00190 * 00191 * apply rotations from the right 00192 * 00193 * 00194 * Dependent on the the number of diagonals either 00195 * DLARTV or DROT is used 00196 * 00197 IF( NR.GE.2*KD-1 ) THEN 00198 DO 10 L = 1, KD - 1 00199 CALL DLARTV( NR, AB( L+1, J1-1 ), INCA, 00200 $ AB( L, J1 ), INCA, D( J1 ), 00201 $ WORK( J1 ), KD1 ) 00202 10 CONTINUE 00203 * 00204 ELSE 00205 JEND = J1 + ( NR-1 )*KD1 00206 DO 20 JINC = J1, JEND, KD1 00207 CALL DROT( KDM1, AB( 2, JINC-1 ), 1, 00208 $ AB( 1, JINC ), 1, D( JINC ), 00209 $ WORK( JINC ) ) 00210 20 CONTINUE 00211 END IF 00212 END IF 00213 * 00214 * 00215 IF( K.GT.2 ) THEN 00216 IF( K.LE.N-I+1 ) THEN 00217 * 00218 * generate plane rotation to annihilate a(i,i+k-1) 00219 * within the band 00220 * 00221 CALL DLARTG( AB( KD-K+3, I+K-2 ), 00222 $ AB( KD-K+2, I+K-1 ), D( I+K-1 ), 00223 $ WORK( I+K-1 ), TEMP ) 00224 AB( KD-K+3, I+K-2 ) = TEMP 00225 * 00226 * apply rotation from the right 00227 * 00228 CALL DROT( K-3, AB( KD-K+4, I+K-2 ), 1, 00229 $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ), 00230 $ WORK( I+K-1 ) ) 00231 END IF 00232 NR = NR + 1 00233 J1 = J1 - KDN - 1 00234 END IF 00235 * 00236 * apply plane rotations from both sides to diagonal 00237 * blocks 00238 * 00239 IF( NR.GT.0 ) 00240 $ CALL DLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ), 00241 $ AB( KD, J1 ), INCA, D( J1 ), 00242 $ WORK( J1 ), KD1 ) 00243 * 00244 * apply plane rotations from the left 00245 * 00246 IF( NR.GT.0 ) THEN 00247 IF( 2*KD-1.LT.NR ) THEN 00248 * 00249 * Dependent on the the number of diagonals either 00250 * DLARTV or DROT is used 00251 * 00252 DO 30 L = 1, KD - 1 00253 IF( J2+L.GT.N ) THEN 00254 NRT = NR - 1 00255 ELSE 00256 NRT = NR 00257 END IF 00258 IF( NRT.GT.0 ) 00259 $ CALL DLARTV( NRT, AB( KD-L, J1+L ), INCA, 00260 $ AB( KD-L+1, J1+L ), INCA, 00261 $ D( J1 ), WORK( J1 ), KD1 ) 00262 30 CONTINUE 00263 ELSE 00264 J1END = J1 + KD1*( NR-2 ) 00265 IF( J1END.GE.J1 ) THEN 00266 DO 40 JIN = J1, J1END, KD1 00267 CALL DROT( KD-1, AB( KD-1, JIN+1 ), INCX, 00268 $ AB( KD, JIN+1 ), INCX, 00269 $ D( JIN ), WORK( JIN ) ) 00270 40 CONTINUE 00271 END IF 00272 LEND = MIN( KDM1, N-J2 ) 00273 LAST = J1END + KD1 00274 IF( LEND.GT.0 ) 00275 $ CALL DROT( LEND, AB( KD-1, LAST+1 ), INCX, 00276 $ AB( KD, LAST+1 ), INCX, D( LAST ), 00277 $ WORK( LAST ) ) 00278 END IF 00279 END IF 00280 * 00281 IF( WANTQ ) THEN 00282 * 00283 * accumulate product of plane rotations in Q 00284 * 00285 IF( INITQ ) THEN 00286 * 00287 * take advantage of the fact that Q was 00288 * initially the Identity matrix 00289 * 00290 IQEND = MAX( IQEND, J2 ) 00291 I2 = MAX( 0, K-3 ) 00292 IQAEND = 1 + I*KD 00293 IF( K.EQ.2 ) 00294 $ IQAEND = IQAEND + KD 00295 IQAEND = MIN( IQAEND, IQEND ) 00296 DO 50 J = J1, J2, KD1 00297 IBL = I - I2 / KDM1 00298 I2 = I2 + 1 00299 IQB = MAX( 1, J-IBL ) 00300 NQ = 1 + IQAEND - IQB 00301 IQAEND = MIN( IQAEND+KD, IQEND ) 00302 CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ), 00303 $ 1, D( J ), WORK( J ) ) 00304 50 CONTINUE 00305 ELSE 00306 * 00307 DO 60 J = J1, J2, KD1 00308 CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1, 00309 $ D( J ), WORK( J ) ) 00310 60 CONTINUE 00311 END IF 00312 * 00313 END IF 00314 * 00315 IF( J2+KDN.GT.N ) THEN 00316 * 00317 * adjust J2 to keep within the bounds of the matrix 00318 * 00319 NR = NR - 1 00320 J2 = J2 - KDN - 1 00321 END IF 00322 * 00323 DO 70 J = J1, J2, KD1 00324 * 00325 * create nonzero element a(j-1,j+kd) outside the band 00326 * and store it in WORK 00327 * 00328 WORK( J+KD ) = WORK( J )*AB( 1, J+KD ) 00329 AB( 1, J+KD ) = D( J )*AB( 1, J+KD ) 00330 70 CONTINUE 00331 80 CONTINUE 00332 90 CONTINUE 00333 END IF 00334 * 00335 IF( KD.GT.0 ) THEN 00336 * 00337 * copy off-diagonal elements to E 00338 * 00339 DO 100 I = 1, N - 1 00340 E( I ) = AB( KD, I+1 ) 00341 100 CONTINUE 00342 ELSE 00343 * 00344 * set E to zero if original matrix was diagonal 00345 * 00346 DO 110 I = 1, N - 1 00347 E( I ) = ZERO 00348 110 CONTINUE 00349 END IF 00350 * 00351 * copy diagonal elements to D 00352 * 00353 DO 120 I = 1, N 00354 D( I ) = AB( KD1, I ) 00355 120 CONTINUE 00356 * 00357 ELSE 00358 * 00359 IF( KD.GT.1 ) THEN 00360 * 00361 * Reduce to tridiagonal form, working with lower triangle 00362 * 00363 NR = 0 00364 J1 = KDN + 2 00365 J2 = 1 00366 * 00367 DO 210 I = 1, N - 2 00368 * 00369 * Reduce i-th column of matrix to tridiagonal form 00370 * 00371 DO 200 K = KDN + 1, 2, -1 00372 J1 = J1 + KDN 00373 J2 = J2 + KDN 00374 * 00375 IF( NR.GT.0 ) THEN 00376 * 00377 * generate plane rotations to annihilate nonzero 00378 * elements which have been created outside the band 00379 * 00380 CALL DLARGV( NR, AB( KD1, J1-KD1 ), INCA, 00381 $ WORK( J1 ), KD1, D( J1 ), KD1 ) 00382 * 00383 * apply plane rotations from one side 00384 * 00385 * 00386 * Dependent on the the number of diagonals either 00387 * DLARTV or DROT is used 00388 * 00389 IF( NR.GT.2*KD-1 ) THEN 00390 DO 130 L = 1, KD - 1 00391 CALL DLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA, 00392 $ AB( KD1-L+1, J1-KD1+L ), INCA, 00393 $ D( J1 ), WORK( J1 ), KD1 ) 00394 130 CONTINUE 00395 ELSE 00396 JEND = J1 + KD1*( NR-1 ) 00397 DO 140 JINC = J1, JEND, KD1 00398 CALL DROT( KDM1, AB( KD, JINC-KD ), INCX, 00399 $ AB( KD1, JINC-KD ), INCX, 00400 $ D( JINC ), WORK( JINC ) ) 00401 140 CONTINUE 00402 END IF 00403 * 00404 END IF 00405 * 00406 IF( K.GT.2 ) THEN 00407 IF( K.LE.N-I+1 ) THEN 00408 * 00409 * generate plane rotation to annihilate a(i+k-1,i) 00410 * within the band 00411 * 00412 CALL DLARTG( AB( K-1, I ), AB( K, I ), 00413 $ D( I+K-1 ), WORK( I+K-1 ), TEMP ) 00414 AB( K-1, I ) = TEMP 00415 * 00416 * apply rotation from the left 00417 * 00418 CALL DROT( K-3, AB( K-2, I+1 ), LDAB-1, 00419 $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ), 00420 $ WORK( I+K-1 ) ) 00421 END IF 00422 NR = NR + 1 00423 J1 = J1 - KDN - 1 00424 END IF 00425 * 00426 * apply plane rotations from both sides to diagonal 00427 * blocks 00428 * 00429 IF( NR.GT.0 ) 00430 $ CALL DLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ), 00431 $ AB( 2, J1-1 ), INCA, D( J1 ), 00432 $ WORK( J1 ), KD1 ) 00433 * 00434 * apply plane rotations from the right 00435 * 00436 * 00437 * Dependent on the the number of diagonals either 00438 * DLARTV or DROT is used 00439 * 00440 IF( NR.GT.0 ) THEN 00441 IF( NR.GT.2*KD-1 ) THEN 00442 DO 150 L = 1, KD - 1 00443 IF( J2+L.GT.N ) THEN 00444 NRT = NR - 1 00445 ELSE 00446 NRT = NR 00447 END IF 00448 IF( NRT.GT.0 ) 00449 $ CALL DLARTV( NRT, AB( L+2, J1-1 ), INCA, 00450 $ AB( L+1, J1 ), INCA, D( J1 ), 00451 $ WORK( J1 ), KD1 ) 00452 150 CONTINUE 00453 ELSE 00454 J1END = J1 + KD1*( NR-2 ) 00455 IF( J1END.GE.J1 ) THEN 00456 DO 160 J1INC = J1, J1END, KD1 00457 CALL DROT( KDM1, AB( 3, J1INC-1 ), 1, 00458 $ AB( 2, J1INC ), 1, D( J1INC ), 00459 $ WORK( J1INC ) ) 00460 160 CONTINUE 00461 END IF 00462 LEND = MIN( KDM1, N-J2 ) 00463 LAST = J1END + KD1 00464 IF( LEND.GT.0 ) 00465 $ CALL DROT( LEND, AB( 3, LAST-1 ), 1, 00466 $ AB( 2, LAST ), 1, D( LAST ), 00467 $ WORK( LAST ) ) 00468 END IF 00469 END IF 00470 * 00471 * 00472 * 00473 IF( WANTQ ) THEN 00474 * 00475 * accumulate product of plane rotations in Q 00476 * 00477 IF( INITQ ) THEN 00478 * 00479 * take advantage of the fact that Q was 00480 * initially the Identity matrix 00481 * 00482 IQEND = MAX( IQEND, J2 ) 00483 I2 = MAX( 0, K-3 ) 00484 IQAEND = 1 + I*KD 00485 IF( K.EQ.2 ) 00486 $ IQAEND = IQAEND + KD 00487 IQAEND = MIN( IQAEND, IQEND ) 00488 DO 170 J = J1, J2, KD1 00489 IBL = I - I2 / KDM1 00490 I2 = I2 + 1 00491 IQB = MAX( 1, J-IBL ) 00492 NQ = 1 + IQAEND - IQB 00493 IQAEND = MIN( IQAEND+KD, IQEND ) 00494 CALL DROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ), 00495 $ 1, D( J ), WORK( J ) ) 00496 170 CONTINUE 00497 ELSE 00498 * 00499 DO 180 J = J1, J2, KD1 00500 CALL DROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1, 00501 $ D( J ), WORK( J ) ) 00502 180 CONTINUE 00503 END IF 00504 END IF 00505 * 00506 IF( J2+KDN.GT.N ) THEN 00507 * 00508 * adjust J2 to keep within the bounds of the matrix 00509 * 00510 NR = NR - 1 00511 J2 = J2 - KDN - 1 00512 END IF 00513 * 00514 DO 190 J = J1, J2, KD1 00515 * 00516 * create nonzero element a(j+kd,j-1) outside the 00517 * band and store it in WORK 00518 * 00519 WORK( J+KD ) = WORK( J )*AB( KD1, J ) 00520 AB( KD1, J ) = D( J )*AB( KD1, J ) 00521 190 CONTINUE 00522 200 CONTINUE 00523 210 CONTINUE 00524 END IF 00525 * 00526 IF( KD.GT.0 ) THEN 00527 * 00528 * copy off-diagonal elements to E 00529 * 00530 DO 220 I = 1, N - 1 00531 E( I ) = AB( 2, I ) 00532 220 CONTINUE 00533 ELSE 00534 * 00535 * set E to zero if original matrix was diagonal 00536 * 00537 DO 230 I = 1, N - 1 00538 E( I ) = ZERO 00539 230 CONTINUE 00540 END IF 00541 * 00542 * copy diagonal elements to D 00543 * 00544 DO 240 I = 1, N 00545 D( I ) = AB( 1, I ) 00546 240 CONTINUE 00547 END IF 00548 * 00549 RETURN 00550 * 00551 * End of DSBTRD 00552 * 00553 END
1.7.3 
