|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CHBTRD( 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 REAL D( * ), E( * ) 00015 COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CHBTRD reduces a complex Hermitian band matrix A to real symmetric 00022 * tridiagonal form T by a unitary similarity transformation: 00023 * Q**H * 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) COMPLEX array, dimension (LDAB,N) 00045 * On entry, the upper or lower triangle of the Hermitian 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) REAL array, dimension (N) 00062 * The diagonal elements of the tridiagonal matrix T. 00063 * 00064 * E (output) REAL 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) COMPLEX 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 unitary 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) COMPLEX 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 REAL ZERO 00096 PARAMETER ( ZERO = 0.0E+0 ) 00097 COMPLEX CZERO, CONE 00098 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00099 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00100 * .. 00101 * .. Local Scalars .. 00102 LOGICAL INITQ, UPPER, WANTQ 00103 INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J, 00104 $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1, 00105 $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT 00106 REAL ABST 00107 COMPLEX T, TEMP 00108 * .. 00109 * .. External Subroutines .. 00110 EXTERNAL CLACGV, CLAR2V, CLARGV, CLARTG, CLARTV, CLASET, 00111 $ CROT, CSCAL, XERBLA 00112 * .. 00113 * .. Intrinsic Functions .. 00114 INTRINSIC ABS, CONJG, MAX, MIN, REAL 00115 * .. 00116 * .. External Functions .. 00117 LOGICAL LSAME 00118 EXTERNAL LSAME 00119 * .. 00120 * .. Executable Statements .. 00121 * 00122 * Test the input parameters 00123 * 00124 INITQ = LSAME( VECT, 'V' ) 00125 WANTQ = INITQ .OR. LSAME( VECT, 'U' ) 00126 UPPER = LSAME( UPLO, 'U' ) 00127 KD1 = KD + 1 00128 KDM1 = KD - 1 00129 INCX = LDAB - 1 00130 IQEND = 1 00131 * 00132 INFO = 0 00133 IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN 00134 INFO = -1 00135 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00136 INFO = -2 00137 ELSE IF( N.LT.0 ) THEN 00138 INFO = -3 00139 ELSE IF( KD.LT.0 ) THEN 00140 INFO = -4 00141 ELSE IF( LDAB.LT.KD1 ) THEN 00142 INFO = -6 00143 ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN 00144 INFO = -10 00145 END IF 00146 IF( INFO.NE.0 ) THEN 00147 CALL XERBLA( 'CHBTRD', -INFO ) 00148 RETURN 00149 END IF 00150 * 00151 * Quick return if possible 00152 * 00153 IF( N.EQ.0 ) 00154 $ RETURN 00155 * 00156 * Initialize Q to the unit matrix, if needed 00157 * 00158 IF( INITQ ) 00159 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 00160 * 00161 * Wherever possible, plane rotations are generated and applied in 00162 * vector operations of length NR over the index set J1:J2:KD1. 00163 * 00164 * The real cosines and complex sines of the plane rotations are 00165 * stored in the arrays D and WORK. 00166 * 00167 INCA = KD1*LDAB 00168 KDN = MIN( N-1, KD ) 00169 IF( UPPER ) THEN 00170 * 00171 IF( KD.GT.1 ) THEN 00172 * 00173 * Reduce to complex Hermitian tridiagonal form, working with 00174 * the upper triangle 00175 * 00176 NR = 0 00177 J1 = KDN + 2 00178 J2 = 1 00179 * 00180 AB( KD1, 1 ) = REAL( AB( KD1, 1 ) ) 00181 DO 90 I = 1, N - 2 00182 * 00183 * Reduce i-th row of matrix to tridiagonal form 00184 * 00185 DO 80 K = KDN + 1, 2, -1 00186 J1 = J1 + KDN 00187 J2 = J2 + KDN 00188 * 00189 IF( NR.GT.0 ) THEN 00190 * 00191 * generate plane rotations to annihilate nonzero 00192 * elements which have been created outside the band 00193 * 00194 CALL CLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ), 00195 $ KD1, D( J1 ), KD1 ) 00196 * 00197 * apply rotations from the right 00198 * 00199 * 00200 * Dependent on the the number of diagonals either 00201 * CLARTV or CROT is used 00202 * 00203 IF( NR.GE.2*KD-1 ) THEN 00204 DO 10 L = 1, KD - 1 00205 CALL CLARTV( NR, AB( L+1, J1-1 ), INCA, 00206 $ AB( L, J1 ), INCA, D( J1 ), 00207 $ WORK( J1 ), KD1 ) 00208 10 CONTINUE 00209 * 00210 ELSE 00211 JEND = J1 + ( NR-1 )*KD1 00212 DO 20 JINC = J1, JEND, KD1 00213 CALL CROT( KDM1, AB( 2, JINC-1 ), 1, 00214 $ AB( 1, JINC ), 1, D( JINC ), 00215 $ WORK( JINC ) ) 00216 20 CONTINUE 00217 END IF 00218 END IF 00219 * 00220 * 00221 IF( K.GT.2 ) THEN 00222 IF( K.LE.N-I+1 ) THEN 00223 * 00224 * generate plane rotation to annihilate a(i,i+k-1) 00225 * within the band 00226 * 00227 CALL CLARTG( AB( KD-K+3, I+K-2 ), 00228 $ AB( KD-K+2, I+K-1 ), D( I+K-1 ), 00229 $ WORK( I+K-1 ), TEMP ) 00230 AB( KD-K+3, I+K-2 ) = TEMP 00231 * 00232 * apply rotation from the right 00233 * 00234 CALL CROT( K-3, AB( KD-K+4, I+K-2 ), 1, 00235 $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ), 00236 $ WORK( I+K-1 ) ) 00237 END IF 00238 NR = NR + 1 00239 J1 = J1 - KDN - 1 00240 END IF 00241 * 00242 * apply plane rotations from both sides to diagonal 00243 * blocks 00244 * 00245 IF( NR.GT.0 ) 00246 $ CALL CLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ), 00247 $ AB( KD, J1 ), INCA, D( J1 ), 00248 $ WORK( J1 ), KD1 ) 00249 * 00250 * apply plane rotations from the left 00251 * 00252 IF( NR.GT.0 ) THEN 00253 CALL CLACGV( NR, WORK( J1 ), KD1 ) 00254 IF( 2*KD-1.LT.NR ) THEN 00255 * 00256 * Dependent on the the number of diagonals either 00257 * CLARTV or CROT is used 00258 * 00259 DO 30 L = 1, KD - 1 00260 IF( J2+L.GT.N ) THEN 00261 NRT = NR - 1 00262 ELSE 00263 NRT = NR 00264 END IF 00265 IF( NRT.GT.0 ) 00266 $ CALL CLARTV( NRT, AB( KD-L, J1+L ), INCA, 00267 $ AB( KD-L+1, J1+L ), INCA, 00268 $ D( J1 ), WORK( J1 ), KD1 ) 00269 30 CONTINUE 00270 ELSE 00271 J1END = J1 + KD1*( NR-2 ) 00272 IF( J1END.GE.J1 ) THEN 00273 DO 40 JIN = J1, J1END, KD1 00274 CALL CROT( KD-1, AB( KD-1, JIN+1 ), INCX, 00275 $ AB( KD, JIN+1 ), INCX, 00276 $ D( JIN ), WORK( JIN ) ) 00277 40 CONTINUE 00278 END IF 00279 LEND = MIN( KDM1, N-J2 ) 00280 LAST = J1END + KD1 00281 IF( LEND.GT.0 ) 00282 $ CALL CROT( LEND, AB( KD-1, LAST+1 ), INCX, 00283 $ AB( KD, LAST+1 ), INCX, D( LAST ), 00284 $ WORK( LAST ) ) 00285 END IF 00286 END IF 00287 * 00288 IF( WANTQ ) THEN 00289 * 00290 * accumulate product of plane rotations in Q 00291 * 00292 IF( INITQ ) THEN 00293 * 00294 * take advantage of the fact that Q was 00295 * initially the Identity matrix 00296 * 00297 IQEND = MAX( IQEND, J2 ) 00298 I2 = MAX( 0, K-3 ) 00299 IQAEND = 1 + I*KD 00300 IF( K.EQ.2 ) 00301 $ IQAEND = IQAEND + KD 00302 IQAEND = MIN( IQAEND, IQEND ) 00303 DO 50 J = J1, J2, KD1 00304 IBL = I - I2 / KDM1 00305 I2 = I2 + 1 00306 IQB = MAX( 1, J-IBL ) 00307 NQ = 1 + IQAEND - IQB 00308 IQAEND = MIN( IQAEND+KD, IQEND ) 00309 CALL CROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ), 00310 $ 1, D( J ), CONJG( WORK( J ) ) ) 00311 50 CONTINUE 00312 ELSE 00313 * 00314 DO 60 J = J1, J2, KD1 00315 CALL CROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1, 00316 $ D( J ), CONJG( WORK( J ) ) ) 00317 60 CONTINUE 00318 END IF 00319 * 00320 END IF 00321 * 00322 IF( J2+KDN.GT.N ) THEN 00323 * 00324 * adjust J2 to keep within the bounds of the matrix 00325 * 00326 NR = NR - 1 00327 J2 = J2 - KDN - 1 00328 END IF 00329 * 00330 DO 70 J = J1, J2, KD1 00331 * 00332 * create nonzero element a(j-1,j+kd) outside the band 00333 * and store it in WORK 00334 * 00335 WORK( J+KD ) = WORK( J )*AB( 1, J+KD ) 00336 AB( 1, J+KD ) = D( J )*AB( 1, J+KD ) 00337 70 CONTINUE 00338 80 CONTINUE 00339 90 CONTINUE 00340 END IF 00341 * 00342 IF( KD.GT.0 ) THEN 00343 * 00344 * make off-diagonal elements real and copy them to E 00345 * 00346 DO 100 I = 1, N - 1 00347 T = AB( KD, I+1 ) 00348 ABST = ABS( T ) 00349 AB( KD, I+1 ) = ABST 00350 E( I ) = ABST 00351 IF( ABST.NE.ZERO ) THEN 00352 T = T / ABST 00353 ELSE 00354 T = CONE 00355 END IF 00356 IF( I.LT.N-1 ) 00357 $ AB( KD, I+2 ) = AB( KD, I+2 )*T 00358 IF( WANTQ ) THEN 00359 CALL CSCAL( N, CONJG( T ), Q( 1, I+1 ), 1 ) 00360 END IF 00361 100 CONTINUE 00362 ELSE 00363 * 00364 * set E to zero if original matrix was diagonal 00365 * 00366 DO 110 I = 1, N - 1 00367 E( I ) = ZERO 00368 110 CONTINUE 00369 END IF 00370 * 00371 * copy diagonal elements to D 00372 * 00373 DO 120 I = 1, N 00374 D( I ) = AB( KD1, I ) 00375 120 CONTINUE 00376 * 00377 ELSE 00378 * 00379 IF( KD.GT.1 ) THEN 00380 * 00381 * Reduce to complex Hermitian tridiagonal form, working with 00382 * the lower triangle 00383 * 00384 NR = 0 00385 J1 = KDN + 2 00386 J2 = 1 00387 * 00388 AB( 1, 1 ) = REAL( AB( 1, 1 ) ) 00389 DO 210 I = 1, N - 2 00390 * 00391 * Reduce i-th column of matrix to tridiagonal form 00392 * 00393 DO 200 K = KDN + 1, 2, -1 00394 J1 = J1 + KDN 00395 J2 = J2 + KDN 00396 * 00397 IF( NR.GT.0 ) THEN 00398 * 00399 * generate plane rotations to annihilate nonzero 00400 * elements which have been created outside the band 00401 * 00402 CALL CLARGV( NR, AB( KD1, J1-KD1 ), INCA, 00403 $ WORK( J1 ), KD1, D( J1 ), KD1 ) 00404 * 00405 * apply plane rotations from one side 00406 * 00407 * 00408 * Dependent on the the number of diagonals either 00409 * CLARTV or CROT is used 00410 * 00411 IF( NR.GT.2*KD-1 ) THEN 00412 DO 130 L = 1, KD - 1 00413 CALL CLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA, 00414 $ AB( KD1-L+1, J1-KD1+L ), INCA, 00415 $ D( J1 ), WORK( J1 ), KD1 ) 00416 130 CONTINUE 00417 ELSE 00418 JEND = J1 + KD1*( NR-1 ) 00419 DO 140 JINC = J1, JEND, KD1 00420 CALL CROT( KDM1, AB( KD, JINC-KD ), INCX, 00421 $ AB( KD1, JINC-KD ), INCX, 00422 $ D( JINC ), WORK( JINC ) ) 00423 140 CONTINUE 00424 END IF 00425 * 00426 END IF 00427 * 00428 IF( K.GT.2 ) THEN 00429 IF( K.LE.N-I+1 ) THEN 00430 * 00431 * generate plane rotation to annihilate a(i+k-1,i) 00432 * within the band 00433 * 00434 CALL CLARTG( AB( K-1, I ), AB( K, I ), 00435 $ D( I+K-1 ), WORK( I+K-1 ), TEMP ) 00436 AB( K-1, I ) = TEMP 00437 * 00438 * apply rotation from the left 00439 * 00440 CALL CROT( K-3, AB( K-2, I+1 ), LDAB-1, 00441 $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ), 00442 $ WORK( I+K-1 ) ) 00443 END IF 00444 NR = NR + 1 00445 J1 = J1 - KDN - 1 00446 END IF 00447 * 00448 * apply plane rotations from both sides to diagonal 00449 * blocks 00450 * 00451 IF( NR.GT.0 ) 00452 $ CALL CLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ), 00453 $ AB( 2, J1-1 ), INCA, D( J1 ), 00454 $ WORK( J1 ), KD1 ) 00455 * 00456 * apply plane rotations from the right 00457 * 00458 * 00459 * Dependent on the the number of diagonals either 00460 * CLARTV or CROT is used 00461 * 00462 IF( NR.GT.0 ) THEN 00463 CALL CLACGV( NR, WORK( J1 ), KD1 ) 00464 IF( NR.GT.2*KD-1 ) THEN 00465 DO 150 L = 1, KD - 1 00466 IF( J2+L.GT.N ) THEN 00467 NRT = NR - 1 00468 ELSE 00469 NRT = NR 00470 END IF 00471 IF( NRT.GT.0 ) 00472 $ CALL CLARTV( NRT, AB( L+2, J1-1 ), INCA, 00473 $ AB( L+1, J1 ), INCA, D( J1 ), 00474 $ WORK( J1 ), KD1 ) 00475 150 CONTINUE 00476 ELSE 00477 J1END = J1 + KD1*( NR-2 ) 00478 IF( J1END.GE.J1 ) THEN 00479 DO 160 J1INC = J1, J1END, KD1 00480 CALL CROT( KDM1, AB( 3, J1INC-1 ), 1, 00481 $ AB( 2, J1INC ), 1, D( J1INC ), 00482 $ WORK( J1INC ) ) 00483 160 CONTINUE 00484 END IF 00485 LEND = MIN( KDM1, N-J2 ) 00486 LAST = J1END + KD1 00487 IF( LEND.GT.0 ) 00488 $ CALL CROT( LEND, AB( 3, LAST-1 ), 1, 00489 $ AB( 2, LAST ), 1, D( LAST ), 00490 $ WORK( LAST ) ) 00491 END IF 00492 END IF 00493 * 00494 * 00495 * 00496 IF( WANTQ ) THEN 00497 * 00498 * accumulate product of plane rotations in Q 00499 * 00500 IF( INITQ ) THEN 00501 * 00502 * take advantage of the fact that Q was 00503 * initially the Identity matrix 00504 * 00505 IQEND = MAX( IQEND, J2 ) 00506 I2 = MAX( 0, K-3 ) 00507 IQAEND = 1 + I*KD 00508 IF( K.EQ.2 ) 00509 $ IQAEND = IQAEND + KD 00510 IQAEND = MIN( IQAEND, IQEND ) 00511 DO 170 J = J1, J2, KD1 00512 IBL = I - I2 / KDM1 00513 I2 = I2 + 1 00514 IQB = MAX( 1, J-IBL ) 00515 NQ = 1 + IQAEND - IQB 00516 IQAEND = MIN( IQAEND+KD, IQEND ) 00517 CALL CROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ), 00518 $ 1, D( J ), WORK( J ) ) 00519 170 CONTINUE 00520 ELSE 00521 * 00522 DO 180 J = J1, J2, KD1 00523 CALL CROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1, 00524 $ D( J ), WORK( J ) ) 00525 180 CONTINUE 00526 END IF 00527 END IF 00528 * 00529 IF( J2+KDN.GT.N ) THEN 00530 * 00531 * adjust J2 to keep within the bounds of the matrix 00532 * 00533 NR = NR - 1 00534 J2 = J2 - KDN - 1 00535 END IF 00536 * 00537 DO 190 J = J1, J2, KD1 00538 * 00539 * create nonzero element a(j+kd,j-1) outside the 00540 * band and store it in WORK 00541 * 00542 WORK( J+KD ) = WORK( J )*AB( KD1, J ) 00543 AB( KD1, J ) = D( J )*AB( KD1, J ) 00544 190 CONTINUE 00545 200 CONTINUE 00546 210 CONTINUE 00547 END IF 00548 * 00549 IF( KD.GT.0 ) THEN 00550 * 00551 * make off-diagonal elements real and copy them to E 00552 * 00553 DO 220 I = 1, N - 1 00554 T = AB( 2, I ) 00555 ABST = ABS( T ) 00556 AB( 2, I ) = ABST 00557 E( I ) = ABST 00558 IF( ABST.NE.ZERO ) THEN 00559 T = T / ABST 00560 ELSE 00561 T = CONE 00562 END IF 00563 IF( I.LT.N-1 ) 00564 $ AB( 2, I+1 ) = AB( 2, I+1 )*T 00565 IF( WANTQ ) THEN 00566 CALL CSCAL( N, T, Q( 1, I+1 ), 1 ) 00567 END IF 00568 220 CONTINUE 00569 ELSE 00570 * 00571 * set E to zero if original matrix was diagonal 00572 * 00573 DO 230 I = 1, N - 1 00574 E( I ) = ZERO 00575 230 CONTINUE 00576 END IF 00577 * 00578 * copy diagonal elements to D 00579 * 00580 DO 240 I = 1, N 00581 D( I ) = AB( 1, I ) 00582 240 CONTINUE 00583 END IF 00584 * 00585 RETURN 00586 * 00587 * End of CHBTRD 00588 * 00589 END
1.7.3 
