![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHBTRD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHBTRD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbtrd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbtrd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbtrd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, 00022 * WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER UPLO, VECT 00026 * INTEGER INFO, KD, LDAB, LDQ, N 00027 * .. 00028 * .. Array Arguments .. 00029 * DOUBLE PRECISION D( * ), E( * ) 00030 * COMPLEX*16 AB( LDAB, * ), Q( LDQ, * ), WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> ZHBTRD reduces a complex Hermitian band matrix A to real symmetric 00040 *> tridiagonal form T by a unitary similarity transformation: 00041 *> Q**H * A * Q = T. 00042 *> \endverbatim 00043 * 00044 * Arguments: 00045 * ========== 00046 * 00047 *> \param[in] VECT 00048 *> \verbatim 00049 *> VECT is CHARACTER*1 00050 *> = 'N': do not form Q; 00051 *> = 'V': form Q; 00052 *> = 'U': update a matrix X, by forming X*Q. 00053 *> \endverbatim 00054 *> 00055 *> \param[in] UPLO 00056 *> \verbatim 00057 *> UPLO is CHARACTER*1 00058 *> = 'U': Upper triangle of A is stored; 00059 *> = 'L': Lower triangle of A is stored. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] N 00063 *> \verbatim 00064 *> N is INTEGER 00065 *> The order of the matrix A. N >= 0. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] KD 00069 *> \verbatim 00070 *> KD is INTEGER 00071 *> The number of superdiagonals of the matrix A if UPLO = 'U', 00072 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0. 00073 *> \endverbatim 00074 *> 00075 *> \param[in,out] AB 00076 *> \verbatim 00077 *> AB is COMPLEX*16 array, dimension (LDAB,N) 00078 *> On entry, the upper or lower triangle of the Hermitian band 00079 *> matrix A, stored in the first KD+1 rows of the array. The 00080 *> j-th column of A is stored in the j-th column of the array AB 00081 *> as follows: 00082 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 00083 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 00084 *> On exit, the diagonal elements of AB are overwritten by the 00085 *> diagonal elements of the tridiagonal matrix T; if KD > 0, the 00086 *> elements on the first superdiagonal (if UPLO = 'U') or the 00087 *> first subdiagonal (if UPLO = 'L') are overwritten by the 00088 *> off-diagonal elements of T; the rest of AB is overwritten by 00089 *> values generated during the reduction. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] LDAB 00093 *> \verbatim 00094 *> LDAB is INTEGER 00095 *> The leading dimension of the array AB. LDAB >= KD+1. 00096 *> \endverbatim 00097 *> 00098 *> \param[out] D 00099 *> \verbatim 00100 *> D is DOUBLE PRECISION array, dimension (N) 00101 *> The diagonal elements of the tridiagonal matrix T. 00102 *> \endverbatim 00103 *> 00104 *> \param[out] E 00105 *> \verbatim 00106 *> E is DOUBLE PRECISION array, dimension (N-1) 00107 *> The off-diagonal elements of the tridiagonal matrix T: 00108 *> E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'. 00109 *> \endverbatim 00110 *> 00111 *> \param[in,out] Q 00112 *> \verbatim 00113 *> Q is COMPLEX*16 array, dimension (LDQ,N) 00114 *> On entry, if VECT = 'U', then Q must contain an N-by-N 00115 *> matrix X; if VECT = 'N' or 'V', then Q need not be set. 00116 *> 00117 *> On exit: 00118 *> if VECT = 'V', Q contains the N-by-N unitary matrix Q; 00119 *> if VECT = 'U', Q contains the product X*Q; 00120 *> if VECT = 'N', the array Q is not referenced. 00121 *> \endverbatim 00122 *> 00123 *> \param[in] LDQ 00124 *> \verbatim 00125 *> LDQ is INTEGER 00126 *> The leading dimension of the array Q. 00127 *> LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'. 00128 *> \endverbatim 00129 *> 00130 *> \param[out] WORK 00131 *> \verbatim 00132 *> WORK is COMPLEX*16 array, dimension (N) 00133 *> \endverbatim 00134 *> 00135 *> \param[out] INFO 00136 *> \verbatim 00137 *> INFO is INTEGER 00138 *> = 0: successful exit 00139 *> < 0: if INFO = -i, the i-th argument had an illegal value 00140 *> \endverbatim 00141 * 00142 * Authors: 00143 * ======== 00144 * 00145 *> \author Univ. of Tennessee 00146 *> \author Univ. of California Berkeley 00147 *> \author Univ. of Colorado Denver 00148 *> \author NAG Ltd. 00149 * 00150 *> \date November 2011 00151 * 00152 *> \ingroup complex16OTHERcomputational 00153 * 00154 *> \par Further Details: 00155 * ===================== 00156 *> 00157 *> \verbatim 00158 *> 00159 *> Modified by Linda Kaufman, Bell Labs. 00160 *> \endverbatim 00161 *> 00162 * ===================================================================== 00163 SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, 00164 $ WORK, INFO ) 00165 * 00166 * -- LAPACK computational routine (version 3.4.0) -- 00167 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00168 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00169 * November 2011 00170 * 00171 * .. Scalar Arguments .. 00172 CHARACTER UPLO, VECT 00173 INTEGER INFO, KD, LDAB, LDQ, N 00174 * .. 00175 * .. Array Arguments .. 00176 DOUBLE PRECISION D( * ), E( * ) 00177 COMPLEX*16 AB( LDAB, * ), Q( LDQ, * ), WORK( * ) 00178 * .. 00179 * 00180 * ===================================================================== 00181 * 00182 * .. Parameters .. 00183 DOUBLE PRECISION ZERO 00184 PARAMETER ( ZERO = 0.0D+0 ) 00185 COMPLEX*16 CZERO, CONE 00186 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00187 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00188 * .. 00189 * .. Local Scalars .. 00190 LOGICAL INITQ, UPPER, WANTQ 00191 INTEGER I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J, 00192 $ J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1, 00193 $ KDM1, KDN, L, LAST, LEND, NQ, NR, NRT 00194 DOUBLE PRECISION ABST 00195 COMPLEX*16 T, TEMP 00196 * .. 00197 * .. External Subroutines .. 00198 EXTERNAL XERBLA, ZLACGV, ZLAR2V, ZLARGV, ZLARTG, ZLARTV, 00199 $ ZLASET, ZROT, ZSCAL 00200 * .. 00201 * .. Intrinsic Functions .. 00202 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN 00203 * .. 00204 * .. External Functions .. 00205 LOGICAL LSAME 00206 EXTERNAL LSAME 00207 * .. 00208 * .. Executable Statements .. 00209 * 00210 * Test the input parameters 00211 * 00212 INITQ = LSAME( VECT, 'V' ) 00213 WANTQ = INITQ .OR. LSAME( VECT, 'U' ) 00214 UPPER = LSAME( UPLO, 'U' ) 00215 KD1 = KD + 1 00216 KDM1 = KD - 1 00217 INCX = LDAB - 1 00218 IQEND = 1 00219 * 00220 INFO = 0 00221 IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN 00222 INFO = -1 00223 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00224 INFO = -2 00225 ELSE IF( N.LT.0 ) THEN 00226 INFO = -3 00227 ELSE IF( KD.LT.0 ) THEN 00228 INFO = -4 00229 ELSE IF( LDAB.LT.KD1 ) THEN 00230 INFO = -6 00231 ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN 00232 INFO = -10 00233 END IF 00234 IF( INFO.NE.0 ) THEN 00235 CALL XERBLA( 'ZHBTRD', -INFO ) 00236 RETURN 00237 END IF 00238 * 00239 * Quick return if possible 00240 * 00241 IF( N.EQ.0 ) 00242 $ RETURN 00243 * 00244 * Initialize Q to the unit matrix, if needed 00245 * 00246 IF( INITQ ) 00247 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 00248 * 00249 * Wherever possible, plane rotations are generated and applied in 00250 * vector operations of length NR over the index set J1:J2:KD1. 00251 * 00252 * The real cosines and complex sines of the plane rotations are 00253 * stored in the arrays D and WORK. 00254 * 00255 INCA = KD1*LDAB 00256 KDN = MIN( N-1, KD ) 00257 IF( UPPER ) THEN 00258 * 00259 IF( KD.GT.1 ) THEN 00260 * 00261 * Reduce to complex Hermitian tridiagonal form, working with 00262 * the upper triangle 00263 * 00264 NR = 0 00265 J1 = KDN + 2 00266 J2 = 1 00267 * 00268 AB( KD1, 1 ) = DBLE( AB( KD1, 1 ) ) 00269 DO 90 I = 1, N - 2 00270 * 00271 * Reduce i-th row of matrix to tridiagonal form 00272 * 00273 DO 80 K = KDN + 1, 2, -1 00274 J1 = J1 + KDN 00275 J2 = J2 + KDN 00276 * 00277 IF( NR.GT.0 ) THEN 00278 * 00279 * generate plane rotations to annihilate nonzero 00280 * elements which have been created outside the band 00281 * 00282 CALL ZLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ), 00283 $ KD1, D( J1 ), KD1 ) 00284 * 00285 * apply rotations from the right 00286 * 00287 * 00288 * Dependent on the the number of diagonals either 00289 * ZLARTV or ZROT is used 00290 * 00291 IF( NR.GE.2*KD-1 ) THEN 00292 DO 10 L = 1, KD - 1 00293 CALL ZLARTV( NR, AB( L+1, J1-1 ), INCA, 00294 $ AB( L, J1 ), INCA, D( J1 ), 00295 $ WORK( J1 ), KD1 ) 00296 10 CONTINUE 00297 * 00298 ELSE 00299 JEND = J1 + ( NR-1 )*KD1 00300 DO 20 JINC = J1, JEND, KD1 00301 CALL ZROT( KDM1, AB( 2, JINC-1 ), 1, 00302 $ AB( 1, JINC ), 1, D( JINC ), 00303 $ WORK( JINC ) ) 00304 20 CONTINUE 00305 END IF 00306 END IF 00307 * 00308 * 00309 IF( K.GT.2 ) THEN 00310 IF( K.LE.N-I+1 ) THEN 00311 * 00312 * generate plane rotation to annihilate a(i,i+k-1) 00313 * within the band 00314 * 00315 CALL ZLARTG( AB( KD-K+3, I+K-2 ), 00316 $ AB( KD-K+2, I+K-1 ), D( I+K-1 ), 00317 $ WORK( I+K-1 ), TEMP ) 00318 AB( KD-K+3, I+K-2 ) = TEMP 00319 * 00320 * apply rotation from the right 00321 * 00322 CALL ZROT( K-3, AB( KD-K+4, I+K-2 ), 1, 00323 $ AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ), 00324 $ WORK( I+K-1 ) ) 00325 END IF 00326 NR = NR + 1 00327 J1 = J1 - KDN - 1 00328 END IF 00329 * 00330 * apply plane rotations from both sides to diagonal 00331 * blocks 00332 * 00333 IF( NR.GT.0 ) 00334 $ CALL ZLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ), 00335 $ AB( KD, J1 ), INCA, D( J1 ), 00336 $ WORK( J1 ), KD1 ) 00337 * 00338 * apply plane rotations from the left 00339 * 00340 IF( NR.GT.0 ) THEN 00341 CALL ZLACGV( NR, WORK( J1 ), KD1 ) 00342 IF( 2*KD-1.LT.NR ) THEN 00343 * 00344 * Dependent on the the number of diagonals either 00345 * ZLARTV or ZROT is used 00346 * 00347 DO 30 L = 1, KD - 1 00348 IF( J2+L.GT.N ) THEN 00349 NRT = NR - 1 00350 ELSE 00351 NRT = NR 00352 END IF 00353 IF( NRT.GT.0 ) 00354 $ CALL ZLARTV( NRT, AB( KD-L, J1+L ), INCA, 00355 $ AB( KD-L+1, J1+L ), INCA, 00356 $ D( J1 ), WORK( J1 ), KD1 ) 00357 30 CONTINUE 00358 ELSE 00359 J1END = J1 + KD1*( NR-2 ) 00360 IF( J1END.GE.J1 ) THEN 00361 DO 40 JIN = J1, J1END, KD1 00362 CALL ZROT( KD-1, AB( KD-1, JIN+1 ), INCX, 00363 $ AB( KD, JIN+1 ), INCX, 00364 $ D( JIN ), WORK( JIN ) ) 00365 40 CONTINUE 00366 END IF 00367 LEND = MIN( KDM1, N-J2 ) 00368 LAST = J1END + KD1 00369 IF( LEND.GT.0 ) 00370 $ CALL ZROT( LEND, AB( KD-1, LAST+1 ), INCX, 00371 $ AB( KD, LAST+1 ), INCX, D( LAST ), 00372 $ WORK( LAST ) ) 00373 END IF 00374 END IF 00375 * 00376 IF( WANTQ ) THEN 00377 * 00378 * accumulate product of plane rotations in Q 00379 * 00380 IF( INITQ ) THEN 00381 * 00382 * take advantage of the fact that Q was 00383 * initially the Identity matrix 00384 * 00385 IQEND = MAX( IQEND, J2 ) 00386 I2 = MAX( 0, K-3 ) 00387 IQAEND = 1 + I*KD 00388 IF( K.EQ.2 ) 00389 $ IQAEND = IQAEND + KD 00390 IQAEND = MIN( IQAEND, IQEND ) 00391 DO 50 J = J1, J2, KD1 00392 IBL = I - I2 / KDM1 00393 I2 = I2 + 1 00394 IQB = MAX( 1, J-IBL ) 00395 NQ = 1 + IQAEND - IQB 00396 IQAEND = MIN( IQAEND+KD, IQEND ) 00397 CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ), 00398 $ 1, D( J ), DCONJG( WORK( J ) ) ) 00399 50 CONTINUE 00400 ELSE 00401 * 00402 DO 60 J = J1, J2, KD1 00403 CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1, 00404 $ D( J ), DCONJG( WORK( J ) ) ) 00405 60 CONTINUE 00406 END IF 00407 * 00408 END IF 00409 * 00410 IF( J2+KDN.GT.N ) THEN 00411 * 00412 * adjust J2 to keep within the bounds of the matrix 00413 * 00414 NR = NR - 1 00415 J2 = J2 - KDN - 1 00416 END IF 00417 * 00418 DO 70 J = J1, J2, KD1 00419 * 00420 * create nonzero element a(j-1,j+kd) outside the band 00421 * and store it in WORK 00422 * 00423 WORK( J+KD ) = WORK( J )*AB( 1, J+KD ) 00424 AB( 1, J+KD ) = D( J )*AB( 1, J+KD ) 00425 70 CONTINUE 00426 80 CONTINUE 00427 90 CONTINUE 00428 END IF 00429 * 00430 IF( KD.GT.0 ) THEN 00431 * 00432 * make off-diagonal elements real and copy them to E 00433 * 00434 DO 100 I = 1, N - 1 00435 T = AB( KD, I+1 ) 00436 ABST = ABS( T ) 00437 AB( KD, I+1 ) = ABST 00438 E( I ) = ABST 00439 IF( ABST.NE.ZERO ) THEN 00440 T = T / ABST 00441 ELSE 00442 T = CONE 00443 END IF 00444 IF( I.LT.N-1 ) 00445 $ AB( KD, I+2 ) = AB( KD, I+2 )*T 00446 IF( WANTQ ) THEN 00447 CALL ZSCAL( N, DCONJG( T ), Q( 1, I+1 ), 1 ) 00448 END IF 00449 100 CONTINUE 00450 ELSE 00451 * 00452 * set E to zero if original matrix was diagonal 00453 * 00454 DO 110 I = 1, N - 1 00455 E( I ) = ZERO 00456 110 CONTINUE 00457 END IF 00458 * 00459 * copy diagonal elements to D 00460 * 00461 DO 120 I = 1, N 00462 D( I ) = AB( KD1, I ) 00463 120 CONTINUE 00464 * 00465 ELSE 00466 * 00467 IF( KD.GT.1 ) THEN 00468 * 00469 * Reduce to complex Hermitian tridiagonal form, working with 00470 * the lower triangle 00471 * 00472 NR = 0 00473 J1 = KDN + 2 00474 J2 = 1 00475 * 00476 AB( 1, 1 ) = DBLE( AB( 1, 1 ) ) 00477 DO 210 I = 1, N - 2 00478 * 00479 * Reduce i-th column of matrix to tridiagonal form 00480 * 00481 DO 200 K = KDN + 1, 2, -1 00482 J1 = J1 + KDN 00483 J2 = J2 + KDN 00484 * 00485 IF( NR.GT.0 ) THEN 00486 * 00487 * generate plane rotations to annihilate nonzero 00488 * elements which have been created outside the band 00489 * 00490 CALL ZLARGV( NR, AB( KD1, J1-KD1 ), INCA, 00491 $ WORK( J1 ), KD1, D( J1 ), KD1 ) 00492 * 00493 * apply plane rotations from one side 00494 * 00495 * 00496 * Dependent on the the number of diagonals either 00497 * ZLARTV or ZROT is used 00498 * 00499 IF( NR.GT.2*KD-1 ) THEN 00500 DO 130 L = 1, KD - 1 00501 CALL ZLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA, 00502 $ AB( KD1-L+1, J1-KD1+L ), INCA, 00503 $ D( J1 ), WORK( J1 ), KD1 ) 00504 130 CONTINUE 00505 ELSE 00506 JEND = J1 + KD1*( NR-1 ) 00507 DO 140 JINC = J1, JEND, KD1 00508 CALL ZROT( KDM1, AB( KD, JINC-KD ), INCX, 00509 $ AB( KD1, JINC-KD ), INCX, 00510 $ D( JINC ), WORK( JINC ) ) 00511 140 CONTINUE 00512 END IF 00513 * 00514 END IF 00515 * 00516 IF( K.GT.2 ) THEN 00517 IF( K.LE.N-I+1 ) THEN 00518 * 00519 * generate plane rotation to annihilate a(i+k-1,i) 00520 * within the band 00521 * 00522 CALL ZLARTG( AB( K-1, I ), AB( K, I ), 00523 $ D( I+K-1 ), WORK( I+K-1 ), TEMP ) 00524 AB( K-1, I ) = TEMP 00525 * 00526 * apply rotation from the left 00527 * 00528 CALL ZROT( K-3, AB( K-2, I+1 ), LDAB-1, 00529 $ AB( K-1, I+1 ), LDAB-1, D( I+K-1 ), 00530 $ WORK( I+K-1 ) ) 00531 END IF 00532 NR = NR + 1 00533 J1 = J1 - KDN - 1 00534 END IF 00535 * 00536 * apply plane rotations from both sides to diagonal 00537 * blocks 00538 * 00539 IF( NR.GT.0 ) 00540 $ CALL ZLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ), 00541 $ AB( 2, J1-1 ), INCA, D( J1 ), 00542 $ WORK( J1 ), KD1 ) 00543 * 00544 * apply plane rotations from the right 00545 * 00546 * 00547 * Dependent on the the number of diagonals either 00548 * ZLARTV or ZROT is used 00549 * 00550 IF( NR.GT.0 ) THEN 00551 CALL ZLACGV( NR, WORK( J1 ), KD1 ) 00552 IF( NR.GT.2*KD-1 ) THEN 00553 DO 150 L = 1, KD - 1 00554 IF( J2+L.GT.N ) THEN 00555 NRT = NR - 1 00556 ELSE 00557 NRT = NR 00558 END IF 00559 IF( NRT.GT.0 ) 00560 $ CALL ZLARTV( NRT, AB( L+2, J1-1 ), INCA, 00561 $ AB( L+1, J1 ), INCA, D( J1 ), 00562 $ WORK( J1 ), KD1 ) 00563 150 CONTINUE 00564 ELSE 00565 J1END = J1 + KD1*( NR-2 ) 00566 IF( J1END.GE.J1 ) THEN 00567 DO 160 J1INC = J1, J1END, KD1 00568 CALL ZROT( KDM1, AB( 3, J1INC-1 ), 1, 00569 $ AB( 2, J1INC ), 1, D( J1INC ), 00570 $ WORK( J1INC ) ) 00571 160 CONTINUE 00572 END IF 00573 LEND = MIN( KDM1, N-J2 ) 00574 LAST = J1END + KD1 00575 IF( LEND.GT.0 ) 00576 $ CALL ZROT( LEND, AB( 3, LAST-1 ), 1, 00577 $ AB( 2, LAST ), 1, D( LAST ), 00578 $ WORK( LAST ) ) 00579 END IF 00580 END IF 00581 * 00582 * 00583 * 00584 IF( WANTQ ) THEN 00585 * 00586 * accumulate product of plane rotations in Q 00587 * 00588 IF( INITQ ) THEN 00589 * 00590 * take advantage of the fact that Q was 00591 * initially the Identity matrix 00592 * 00593 IQEND = MAX( IQEND, J2 ) 00594 I2 = MAX( 0, K-3 ) 00595 IQAEND = 1 + I*KD 00596 IF( K.EQ.2 ) 00597 $ IQAEND = IQAEND + KD 00598 IQAEND = MIN( IQAEND, IQEND ) 00599 DO 170 J = J1, J2, KD1 00600 IBL = I - I2 / KDM1 00601 I2 = I2 + 1 00602 IQB = MAX( 1, J-IBL ) 00603 NQ = 1 + IQAEND - IQB 00604 IQAEND = MIN( IQAEND+KD, IQEND ) 00605 CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ), 00606 $ 1, D( J ), WORK( J ) ) 00607 170 CONTINUE 00608 ELSE 00609 * 00610 DO 180 J = J1, J2, KD1 00611 CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1, 00612 $ D( J ), WORK( J ) ) 00613 180 CONTINUE 00614 END IF 00615 END IF 00616 * 00617 IF( J2+KDN.GT.N ) THEN 00618 * 00619 * adjust J2 to keep within the bounds of the matrix 00620 * 00621 NR = NR - 1 00622 J2 = J2 - KDN - 1 00623 END IF 00624 * 00625 DO 190 J = J1, J2, KD1 00626 * 00627 * create nonzero element a(j+kd,j-1) outside the 00628 * band and store it in WORK 00629 * 00630 WORK( J+KD ) = WORK( J )*AB( KD1, J ) 00631 AB( KD1, J ) = D( J )*AB( KD1, J ) 00632 190 CONTINUE 00633 200 CONTINUE 00634 210 CONTINUE 00635 END IF 00636 * 00637 IF( KD.GT.0 ) THEN 00638 * 00639 * make off-diagonal elements real and copy them to E 00640 * 00641 DO 220 I = 1, N - 1 00642 T = AB( 2, I ) 00643 ABST = ABS( T ) 00644 AB( 2, I ) = ABST 00645 E( I ) = ABST 00646 IF( ABST.NE.ZERO ) THEN 00647 T = T / ABST 00648 ELSE 00649 T = CONE 00650 END IF 00651 IF( I.LT.N-1 ) 00652 $ AB( 2, I+1 ) = AB( 2, I+1 )*T 00653 IF( WANTQ ) THEN 00654 CALL ZSCAL( N, T, Q( 1, I+1 ), 1 ) 00655 END IF 00656 220 CONTINUE 00657 ELSE 00658 * 00659 * set E to zero if original matrix was diagonal 00660 * 00661 DO 230 I = 1, N - 1 00662 E( I ) = ZERO 00663 230 CONTINUE 00664 END IF 00665 * 00666 * copy diagonal elements to D 00667 * 00668 DO 240 I = 1, N 00669 D( I ) = AB( 1, I ) 00670 240 CONTINUE 00671 END IF 00672 * 00673 RETURN 00674 * 00675 * End of ZHBTRD 00676 * 00677 END