|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER INFO, KB, LDA, LDW, N, NB 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IPIV( * ) 00014 COMPLEX*16 A( LDA, * ), W( LDW, * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * ZLAHEF computes a partial factorization of a complex Hermitian 00021 * matrix A using the Bunch-Kaufman diagonal pivoting method. The 00022 * partial factorization has the form: 00023 * 00024 * A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or: 00025 * ( 0 U22 ) ( 0 D ) ( U12' U22' ) 00026 * 00027 * A = ( L11 0 ) ( D 0 ) ( L11' L21' ) if UPLO = 'L' 00028 * ( L21 I ) ( 0 A22 ) ( 0 I ) 00029 * 00030 * where the order of D is at most NB. The actual order is returned in 00031 * the argument KB, and is either NB or NB-1, or N if N <= NB. 00032 * Note that U' denotes the conjugate transpose of U. 00033 * 00034 * ZLAHEF is an auxiliary routine called by ZHETRF. It uses blocked code 00035 * (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or 00036 * A22 (if UPLO = 'L'). 00037 * 00038 * Arguments 00039 * ========= 00040 * 00041 * UPLO (input) CHARACTER*1 00042 * Specifies whether the upper or lower triangular part of the 00043 * Hermitian matrix A is stored: 00044 * = 'U': Upper triangular 00045 * = 'L': Lower triangular 00046 * 00047 * N (input) INTEGER 00048 * The order of the matrix A. N >= 0. 00049 * 00050 * NB (input) INTEGER 00051 * The maximum number of columns of the matrix A that should be 00052 * factored. NB should be at least 2 to allow for 2-by-2 pivot 00053 * blocks. 00054 * 00055 * KB (output) INTEGER 00056 * The number of columns of A that were actually factored. 00057 * KB is either NB-1 or NB, or N if N <= NB. 00058 * 00059 * A (input/output) COMPLEX*16 array, dimension (LDA,N) 00060 * On entry, the Hermitian matrix A. If UPLO = 'U', the leading 00061 * n-by-n upper triangular part of A contains the upper 00062 * triangular part of the matrix A, and the strictly lower 00063 * triangular part of A is not referenced. If UPLO = 'L', the 00064 * leading n-by-n lower triangular part of A contains the lower 00065 * triangular part of the matrix A, and the strictly upper 00066 * triangular part of A is not referenced. 00067 * On exit, A contains details of the partial factorization. 00068 * 00069 * LDA (input) INTEGER 00070 * The leading dimension of the array A. LDA >= max(1,N). 00071 * 00072 * IPIV (output) INTEGER array, dimension (N) 00073 * Details of the interchanges and the block structure of D. 00074 * If UPLO = 'U', only the last KB elements of IPIV are set; 00075 * if UPLO = 'L', only the first KB elements are set. 00076 * 00077 * If IPIV(k) > 0, then rows and columns k and IPIV(k) were 00078 * interchanged and D(k,k) is a 1-by-1 diagonal block. 00079 * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 00080 * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 00081 * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 00082 * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 00083 * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 00084 * 00085 * W (workspace) COMPLEX*16 array, dimension (LDW,NB) 00086 * 00087 * LDW (input) INTEGER 00088 * The leading dimension of the array W. LDW >= max(1,N). 00089 * 00090 * INFO (output) INTEGER 00091 * = 0: successful exit 00092 * > 0: if INFO = k, D(k,k) is exactly zero. The factorization 00093 * has been completed, but the block diagonal matrix D is 00094 * exactly singular. 00095 * 00096 * ===================================================================== 00097 * 00098 * .. Parameters .. 00099 DOUBLE PRECISION ZERO, ONE 00100 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00101 COMPLEX*16 CONE 00102 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 00103 DOUBLE PRECISION EIGHT, SEVTEN 00104 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 ) 00105 * .. 00106 * .. Local Scalars .. 00107 INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP, 00108 $ KSTEP, KW 00109 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, R1, ROWMAX, T 00110 COMPLEX*16 D11, D21, D22, Z 00111 * .. 00112 * .. External Functions .. 00113 LOGICAL LSAME 00114 INTEGER IZAMAX 00115 EXTERNAL LSAME, IZAMAX 00116 * .. 00117 * .. External Subroutines .. 00118 EXTERNAL ZCOPY, ZDSCAL, ZGEMM, ZGEMV, ZLACGV, ZSWAP 00119 * .. 00120 * .. Intrinsic Functions .. 00121 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT 00122 * .. 00123 * .. Statement Functions .. 00124 DOUBLE PRECISION CABS1 00125 * .. 00126 * .. Statement Function definitions .. 00127 CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) ) 00128 * .. 00129 * .. Executable Statements .. 00130 * 00131 INFO = 0 00132 * 00133 * Initialize ALPHA for use in choosing pivot block size. 00134 * 00135 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 00136 * 00137 IF( LSAME( UPLO, 'U' ) ) THEN 00138 * 00139 * Factorize the trailing columns of A using the upper triangle 00140 * of A and working backwards, and compute the matrix W = U12*D 00141 * for use in updating A11 (note that conjg(W) is actually stored) 00142 * 00143 * K is the main loop index, decreasing from N in steps of 1 or 2 00144 * 00145 * KW is the column of W which corresponds to column K of A 00146 * 00147 K = N 00148 10 CONTINUE 00149 KW = NB + K - N 00150 * 00151 * Exit from loop 00152 * 00153 IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 ) 00154 $ GO TO 30 00155 * 00156 * Copy column K of A to column KW of W and update it 00157 * 00158 CALL ZCOPY( K-1, A( 1, K ), 1, W( 1, KW ), 1 ) 00159 W( K, KW ) = DBLE( A( K, K ) ) 00160 IF( K.LT.N ) THEN 00161 CALL ZGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA, 00162 $ W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 ) 00163 W( K, KW ) = DBLE( W( K, KW ) ) 00164 END IF 00165 * 00166 KSTEP = 1 00167 * 00168 * Determine rows and columns to be interchanged and whether 00169 * a 1-by-1 or 2-by-2 pivot block will be used 00170 * 00171 ABSAKK = ABS( DBLE( W( K, KW ) ) ) 00172 * 00173 * IMAX is the row-index of the largest off-diagonal element in 00174 * column K, and COLMAX is its absolute value 00175 * 00176 IF( K.GT.1 ) THEN 00177 IMAX = IZAMAX( K-1, W( 1, KW ), 1 ) 00178 COLMAX = CABS1( W( IMAX, KW ) ) 00179 ELSE 00180 COLMAX = ZERO 00181 END IF 00182 * 00183 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 00184 * 00185 * Column K is zero: set INFO and continue 00186 * 00187 IF( INFO.EQ.0 ) 00188 $ INFO = K 00189 KP = K 00190 A( K, K ) = DBLE( A( K, K ) ) 00191 ELSE 00192 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 00193 * 00194 * no interchange, use 1-by-1 pivot block 00195 * 00196 KP = K 00197 ELSE 00198 * 00199 * Copy column IMAX to column KW-1 of W and update it 00200 * 00201 CALL ZCOPY( IMAX-1, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 ) 00202 W( IMAX, KW-1 ) = DBLE( A( IMAX, IMAX ) ) 00203 CALL ZCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA, 00204 $ W( IMAX+1, KW-1 ), 1 ) 00205 CALL ZLACGV( K-IMAX, W( IMAX+1, KW-1 ), 1 ) 00206 IF( K.LT.N ) THEN 00207 CALL ZGEMV( 'No transpose', K, N-K, -CONE, 00208 $ A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW, 00209 $ CONE, W( 1, KW-1 ), 1 ) 00210 W( IMAX, KW-1 ) = DBLE( W( IMAX, KW-1 ) ) 00211 END IF 00212 * 00213 * JMAX is the column-index of the largest off-diagonal 00214 * element in row IMAX, and ROWMAX is its absolute value 00215 * 00216 JMAX = IMAX + IZAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 ) 00217 ROWMAX = CABS1( W( JMAX, KW-1 ) ) 00218 IF( IMAX.GT.1 ) THEN 00219 JMAX = IZAMAX( IMAX-1, W( 1, KW-1 ), 1 ) 00220 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) ) 00221 END IF 00222 * 00223 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 00224 * 00225 * no interchange, use 1-by-1 pivot block 00226 * 00227 KP = K 00228 ELSE IF( ABS( DBLE( W( IMAX, KW-1 ) ) ).GE.ALPHA*ROWMAX ) 00229 $ THEN 00230 * 00231 * interchange rows and columns K and IMAX, use 1-by-1 00232 * pivot block 00233 * 00234 KP = IMAX 00235 * 00236 * copy column KW-1 of W to column KW 00237 * 00238 CALL ZCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 ) 00239 ELSE 00240 * 00241 * interchange rows and columns K-1 and IMAX, use 2-by-2 00242 * pivot block 00243 * 00244 KP = IMAX 00245 KSTEP = 2 00246 END IF 00247 END IF 00248 * 00249 KK = K - KSTEP + 1 00250 KKW = NB + KK - N 00251 * 00252 * Updated column KP is already stored in column KKW of W 00253 * 00254 IF( KP.NE.KK ) THEN 00255 * 00256 * Copy non-updated column KK to column KP 00257 * 00258 A( KP, KP ) = DBLE( A( KK, KK ) ) 00259 CALL ZCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ), 00260 $ LDA ) 00261 CALL ZLACGV( KK-1-KP, A( KP, KP+1 ), LDA ) 00262 CALL ZCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 ) 00263 * 00264 * Interchange rows KK and KP in last KK columns of A and W 00265 * 00266 IF( KK.LT.N ) 00267 $ CALL ZSWAP( N-KK, A( KK, KK+1 ), LDA, A( KP, KK+1 ), 00268 $ LDA ) 00269 CALL ZSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ), 00270 $ LDW ) 00271 END IF 00272 * 00273 IF( KSTEP.EQ.1 ) THEN 00274 * 00275 * 1-by-1 pivot block D(k): column KW of W now holds 00276 * 00277 * W(k) = U(k)*D(k) 00278 * 00279 * where U(k) is the k-th column of U 00280 * 00281 * Store U(k) in column k of A 00282 * 00283 CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 ) 00284 R1 = ONE / DBLE( A( K, K ) ) 00285 CALL ZDSCAL( K-1, R1, A( 1, K ), 1 ) 00286 * 00287 * Conjugate W(k) 00288 * 00289 CALL ZLACGV( K-1, W( 1, KW ), 1 ) 00290 ELSE 00291 * 00292 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now 00293 * hold 00294 * 00295 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 00296 * 00297 * where U(k) and U(k-1) are the k-th and (k-1)-th columns 00298 * of U 00299 * 00300 IF( K.GT.2 ) THEN 00301 * 00302 * Store U(k) and U(k-1) in columns k and k-1 of A 00303 * 00304 D21 = W( K-1, KW ) 00305 D11 = W( K, KW ) / DCONJG( D21 ) 00306 D22 = W( K-1, KW-1 ) / D21 00307 T = ONE / ( DBLE( D11*D22 )-ONE ) 00308 D21 = T / D21 00309 DO 20 J = 1, K - 2 00310 A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) ) 00311 A( J, K ) = DCONJG( D21 )* 00312 $ ( D22*W( J, KW )-W( J, KW-1 ) ) 00313 20 CONTINUE 00314 END IF 00315 * 00316 * Copy D(k) to A 00317 * 00318 A( K-1, K-1 ) = W( K-1, KW-1 ) 00319 A( K-1, K ) = W( K-1, KW ) 00320 A( K, K ) = W( K, KW ) 00321 * 00322 * Conjugate W(k) and W(k-1) 00323 * 00324 CALL ZLACGV( K-1, W( 1, KW ), 1 ) 00325 CALL ZLACGV( K-2, W( 1, KW-1 ), 1 ) 00326 END IF 00327 END IF 00328 * 00329 * Store details of the interchanges in IPIV 00330 * 00331 IF( KSTEP.EQ.1 ) THEN 00332 IPIV( K ) = KP 00333 ELSE 00334 IPIV( K ) = -KP 00335 IPIV( K-1 ) = -KP 00336 END IF 00337 * 00338 * Decrease K and return to the start of the main loop 00339 * 00340 K = K - KSTEP 00341 GO TO 10 00342 * 00343 30 CONTINUE 00344 * 00345 * Update the upper triangle of A11 (= A(1:k,1:k)) as 00346 * 00347 * A11 := A11 - U12*D*U12' = A11 - U12*W' 00348 * 00349 * computing blocks of NB columns at a time (note that conjg(W) is 00350 * actually stored) 00351 * 00352 DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB 00353 JB = MIN( NB, K-J+1 ) 00354 * 00355 * Update the upper triangle of the diagonal block 00356 * 00357 DO 40 JJ = J, J + JB - 1 00358 A( JJ, JJ ) = DBLE( A( JJ, JJ ) ) 00359 CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE, 00360 $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE, 00361 $ A( J, JJ ), 1 ) 00362 A( JJ, JJ ) = DBLE( A( JJ, JJ ) ) 00363 40 CONTINUE 00364 * 00365 * Update the rectangular superdiagonal block 00366 * 00367 CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K, 00368 $ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, 00369 $ CONE, A( 1, J ), LDA ) 00370 50 CONTINUE 00371 * 00372 * Put U12 in standard form by partially undoing the interchanges 00373 * in columns k+1:n 00374 * 00375 J = K + 1 00376 60 CONTINUE 00377 JJ = J 00378 JP = IPIV( J ) 00379 IF( JP.LT.0 ) THEN 00380 JP = -JP 00381 J = J + 1 00382 END IF 00383 J = J + 1 00384 IF( JP.NE.JJ .AND. J.LE.N ) 00385 $ CALL ZSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA ) 00386 IF( J.LE.N ) 00387 $ GO TO 60 00388 * 00389 * Set KB to the number of columns factorized 00390 * 00391 KB = N - K 00392 * 00393 ELSE 00394 * 00395 * Factorize the leading columns of A using the lower triangle 00396 * of A and working forwards, and compute the matrix W = L21*D 00397 * for use in updating A22 (note that conjg(W) is actually stored) 00398 * 00399 * K is the main loop index, increasing from 1 in steps of 1 or 2 00400 * 00401 K = 1 00402 70 CONTINUE 00403 * 00404 * Exit from loop 00405 * 00406 IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N ) 00407 $ GO TO 90 00408 * 00409 * Copy column K of A to column K of W and update it 00410 * 00411 W( K, K ) = DBLE( A( K, K ) ) 00412 IF( K.LT.N ) 00413 $ CALL ZCOPY( N-K, A( K+1, K ), 1, W( K+1, K ), 1 ) 00414 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA, 00415 $ W( K, 1 ), LDW, CONE, W( K, K ), 1 ) 00416 W( K, K ) = DBLE( W( K, K ) ) 00417 * 00418 KSTEP = 1 00419 * 00420 * Determine rows and columns to be interchanged and whether 00421 * a 1-by-1 or 2-by-2 pivot block will be used 00422 * 00423 ABSAKK = ABS( DBLE( W( K, K ) ) ) 00424 * 00425 * IMAX is the row-index of the largest off-diagonal element in 00426 * column K, and COLMAX is its absolute value 00427 * 00428 IF( K.LT.N ) THEN 00429 IMAX = K + IZAMAX( N-K, W( K+1, K ), 1 ) 00430 COLMAX = CABS1( W( IMAX, K ) ) 00431 ELSE 00432 COLMAX = ZERO 00433 END IF 00434 * 00435 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 00436 * 00437 * Column K is zero: set INFO and continue 00438 * 00439 IF( INFO.EQ.0 ) 00440 $ INFO = K 00441 KP = K 00442 A( K, K ) = DBLE( A( K, K ) ) 00443 ELSE 00444 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 00445 * 00446 * no interchange, use 1-by-1 pivot block 00447 * 00448 KP = K 00449 ELSE 00450 * 00451 * Copy column IMAX to column K+1 of W and update it 00452 * 00453 CALL ZCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 ) 00454 CALL ZLACGV( IMAX-K, W( K, K+1 ), 1 ) 00455 W( IMAX, K+1 ) = DBLE( A( IMAX, IMAX ) ) 00456 IF( IMAX.LT.N ) 00457 $ CALL ZCOPY( N-IMAX, A( IMAX+1, IMAX ), 1, 00458 $ W( IMAX+1, K+1 ), 1 ) 00459 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), 00460 $ LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ), 00461 $ 1 ) 00462 W( IMAX, K+1 ) = DBLE( W( IMAX, K+1 ) ) 00463 * 00464 * JMAX is the column-index of the largest off-diagonal 00465 * element in row IMAX, and ROWMAX is its absolute value 00466 * 00467 JMAX = K - 1 + IZAMAX( IMAX-K, W( K, K+1 ), 1 ) 00468 ROWMAX = CABS1( W( JMAX, K+1 ) ) 00469 IF( IMAX.LT.N ) THEN 00470 JMAX = IMAX + IZAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 ) 00471 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) ) 00472 END IF 00473 * 00474 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 00475 * 00476 * no interchange, use 1-by-1 pivot block 00477 * 00478 KP = K 00479 ELSE IF( ABS( DBLE( W( IMAX, K+1 ) ) ).GE.ALPHA*ROWMAX ) 00480 $ THEN 00481 * 00482 * interchange rows and columns K and IMAX, use 1-by-1 00483 * pivot block 00484 * 00485 KP = IMAX 00486 * 00487 * copy column K+1 of W to column K 00488 * 00489 CALL ZCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 ) 00490 ELSE 00491 * 00492 * interchange rows and columns K+1 and IMAX, use 2-by-2 00493 * pivot block 00494 * 00495 KP = IMAX 00496 KSTEP = 2 00497 END IF 00498 END IF 00499 * 00500 KK = K + KSTEP - 1 00501 * 00502 * Updated column KP is already stored in column KK of W 00503 * 00504 IF( KP.NE.KK ) THEN 00505 * 00506 * Copy non-updated column KK to column KP 00507 * 00508 A( KP, KP ) = DBLE( A( KK, KK ) ) 00509 CALL ZCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ), 00510 $ LDA ) 00511 CALL ZLACGV( KP-KK-1, A( KP, KK+1 ), LDA ) 00512 IF( KP.LT.N ) 00513 $ CALL ZCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 ) 00514 * 00515 * Interchange rows KK and KP in first KK columns of A and W 00516 * 00517 CALL ZSWAP( KK-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA ) 00518 CALL ZSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW ) 00519 END IF 00520 * 00521 IF( KSTEP.EQ.1 ) THEN 00522 * 00523 * 1-by-1 pivot block D(k): column k of W now holds 00524 * 00525 * W(k) = L(k)*D(k) 00526 * 00527 * where L(k) is the k-th column of L 00528 * 00529 * Store L(k) in column k of A 00530 * 00531 CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 ) 00532 IF( K.LT.N ) THEN 00533 R1 = ONE / DBLE( A( K, K ) ) 00534 CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 ) 00535 * 00536 * Conjugate W(k) 00537 * 00538 CALL ZLACGV( N-K, W( K+1, K ), 1 ) 00539 END IF 00540 ELSE 00541 * 00542 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold 00543 * 00544 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) 00545 * 00546 * where L(k) and L(k+1) are the k-th and (k+1)-th columns 00547 * of L 00548 * 00549 IF( K.LT.N-1 ) THEN 00550 * 00551 * Store L(k) and L(k+1) in columns k and k+1 of A 00552 * 00553 D21 = W( K+1, K ) 00554 D11 = W( K+1, K+1 ) / D21 00555 D22 = W( K, K ) / DCONJG( D21 ) 00556 T = ONE / ( DBLE( D11*D22 )-ONE ) 00557 D21 = T / D21 00558 DO 80 J = K + 2, N 00559 A( J, K ) = DCONJG( D21 )* 00560 $ ( D11*W( J, K )-W( J, K+1 ) ) 00561 A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) ) 00562 80 CONTINUE 00563 END IF 00564 * 00565 * Copy D(k) to A 00566 * 00567 A( K, K ) = W( K, K ) 00568 A( K+1, K ) = W( K+1, K ) 00569 A( K+1, K+1 ) = W( K+1, K+1 ) 00570 * 00571 * Conjugate W(k) and W(k+1) 00572 * 00573 CALL ZLACGV( N-K, W( K+1, K ), 1 ) 00574 CALL ZLACGV( N-K-1, W( K+2, K+1 ), 1 ) 00575 END IF 00576 END IF 00577 * 00578 * Store details of the interchanges in IPIV 00579 * 00580 IF( KSTEP.EQ.1 ) THEN 00581 IPIV( K ) = KP 00582 ELSE 00583 IPIV( K ) = -KP 00584 IPIV( K+1 ) = -KP 00585 END IF 00586 * 00587 * Increase K and return to the start of the main loop 00588 * 00589 K = K + KSTEP 00590 GO TO 70 00591 * 00592 90 CONTINUE 00593 * 00594 * Update the lower triangle of A22 (= A(k:n,k:n)) as 00595 * 00596 * A22 := A22 - L21*D*L21' = A22 - L21*W' 00597 * 00598 * computing blocks of NB columns at a time (note that conjg(W) is 00599 * actually stored) 00600 * 00601 DO 110 J = K, N, NB 00602 JB = MIN( NB, N-J+1 ) 00603 * 00604 * Update the lower triangle of the diagonal block 00605 * 00606 DO 100 JJ = J, J + JB - 1 00607 A( JJ, JJ ) = DBLE( A( JJ, JJ ) ) 00608 CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE, 00609 $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE, 00610 $ A( JJ, JJ ), 1 ) 00611 A( JJ, JJ ) = DBLE( A( JJ, JJ ) ) 00612 100 CONTINUE 00613 * 00614 * Update the rectangular subdiagonal block 00615 * 00616 IF( J+JB.LE.N ) 00617 $ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB, 00618 $ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ), 00619 $ LDW, CONE, A( J+JB, J ), LDA ) 00620 110 CONTINUE 00621 * 00622 * Put L21 in standard form by partially undoing the interchanges 00623 * in columns 1:k-1 00624 * 00625 J = K - 1 00626 120 CONTINUE 00627 JJ = J 00628 JP = IPIV( J ) 00629 IF( JP.LT.0 ) THEN 00630 JP = -JP 00631 J = J - 1 00632 END IF 00633 J = J - 1 00634 IF( JP.NE.JJ .AND. J.GE.1 ) 00635 $ CALL ZSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA ) 00636 IF( J.GE.1 ) 00637 $ GO TO 120 00638 * 00639 * Set KB to the number of columns factorized 00640 * 00641 KB = K - 1 00642 * 00643 END IF 00644 RETURN 00645 * 00646 * End of ZLAHEF 00647 * 00648 END
1.7.3 
