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