![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLASYF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZLASYF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlasyf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlasyf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlasyf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZLASYF( 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 *> ZLASYF computes a partial factorization of a complex symmetric matrix 00039 *> A using the Bunch-Kaufman diagonal pivoting method. The partial 00040 *> factorization has the form: 00041 *> 00042 *> A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or: 00043 *> ( 0 U22 ) ( 0 D ) ( U12**T U22**T ) 00044 *> 00045 *> A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) 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**T denotes the transpose of U. 00051 *> 00052 *> ZLASYF is an auxiliary routine called by ZSYTRF. 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 *> symmetric 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 symmetric 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 complex16SYcomputational 00156 * 00157 * ===================================================================== 00158 SUBROUTINE ZLASYF( 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 DOUBLE PRECISION EIGHT, SEVTEN 00180 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 ) 00181 COMPLEX*16 CONE 00182 PARAMETER ( CONE = ( 1.0D+0, 0.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, ROWMAX 00188 COMPLEX*16 D11, D21, D22, R1, T, Z 00189 * .. 00190 * .. External Functions .. 00191 LOGICAL LSAME 00192 INTEGER IZAMAX 00193 EXTERNAL LSAME, IZAMAX 00194 * .. 00195 * .. External Subroutines .. 00196 EXTERNAL ZCOPY, ZGEMM, ZGEMV, ZSCAL, ZSWAP 00197 * .. 00198 * .. Intrinsic Functions .. 00199 INTRINSIC ABS, DBLE, 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 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, A( 1, K ), 1, W( 1, KW ), 1 ) 00237 IF( K.LT.N ) 00238 $ CALL ZGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA, 00239 $ W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 ) 00240 * 00241 KSTEP = 1 00242 * 00243 * Determine rows and columns to be interchanged and whether 00244 * a 1-by-1 or 2-by-2 pivot block will be used 00245 * 00246 ABSAKK = CABS1( W( K, KW ) ) 00247 * 00248 * IMAX is the row-index of the largest off-diagonal element in 00249 * column K, and COLMAX is its absolute value 00250 * 00251 IF( K.GT.1 ) THEN 00252 IMAX = IZAMAX( K-1, W( 1, KW ), 1 ) 00253 COLMAX = CABS1( W( IMAX, KW ) ) 00254 ELSE 00255 COLMAX = ZERO 00256 END IF 00257 * 00258 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 00259 * 00260 * Column K is zero: set INFO and continue 00261 * 00262 IF( INFO.EQ.0 ) 00263 $ INFO = K 00264 KP = K 00265 ELSE 00266 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 00267 * 00268 * no interchange, use 1-by-1 pivot block 00269 * 00270 KP = K 00271 ELSE 00272 * 00273 * Copy column IMAX to column KW-1 of W and update it 00274 * 00275 CALL ZCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 ) 00276 CALL ZCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA, 00277 $ W( IMAX+1, KW-1 ), 1 ) 00278 IF( K.LT.N ) 00279 $ CALL ZGEMV( 'No transpose', K, N-K, -CONE, 00280 $ A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW, 00281 $ CONE, W( 1, KW-1 ), 1 ) 00282 * 00283 * JMAX is the column-index of the largest off-diagonal 00284 * element in row IMAX, and ROWMAX is its absolute value 00285 * 00286 JMAX = IMAX + IZAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 ) 00287 ROWMAX = CABS1( W( JMAX, KW-1 ) ) 00288 IF( IMAX.GT.1 ) THEN 00289 JMAX = IZAMAX( IMAX-1, W( 1, KW-1 ), 1 ) 00290 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) ) 00291 END IF 00292 * 00293 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 00294 * 00295 * no interchange, use 1-by-1 pivot block 00296 * 00297 KP = K 00298 ELSE IF( CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN 00299 * 00300 * interchange rows and columns K and IMAX, use 1-by-1 00301 * pivot block 00302 * 00303 KP = IMAX 00304 * 00305 * copy column KW-1 of W to column KW 00306 * 00307 CALL ZCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 ) 00308 ELSE 00309 * 00310 * interchange rows and columns K-1 and IMAX, use 2-by-2 00311 * pivot block 00312 * 00313 KP = IMAX 00314 KSTEP = 2 00315 END IF 00316 END IF 00317 * 00318 KK = K - KSTEP + 1 00319 KKW = NB + KK - N 00320 * 00321 * Updated column KP is already stored in column KKW of W 00322 * 00323 IF( KP.NE.KK ) THEN 00324 * 00325 * Copy non-updated column KK to column KP 00326 * 00327 A( KP, K ) = A( KK, K ) 00328 CALL ZCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ), 00329 $ LDA ) 00330 CALL ZCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 ) 00331 * 00332 * Interchange rows KK and KP in last KK columns of A and W 00333 * 00334 CALL ZSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA ) 00335 CALL ZSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ), 00336 $ LDW ) 00337 END IF 00338 * 00339 IF( KSTEP.EQ.1 ) THEN 00340 * 00341 * 1-by-1 pivot block D(k): column KW of W now holds 00342 * 00343 * W(k) = U(k)*D(k) 00344 * 00345 * where U(k) is the k-th column of U 00346 * 00347 * Store U(k) in column k of A 00348 * 00349 CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 ) 00350 R1 = CONE / A( K, K ) 00351 CALL ZSCAL( K-1, R1, A( 1, K ), 1 ) 00352 ELSE 00353 * 00354 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now 00355 * hold 00356 * 00357 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 00358 * 00359 * where U(k) and U(k-1) are the k-th and (k-1)-th columns 00360 * of U 00361 * 00362 IF( K.GT.2 ) THEN 00363 * 00364 * Store U(k) and U(k-1) in columns k and k-1 of A 00365 * 00366 D21 = W( K-1, KW ) 00367 D11 = W( K, KW ) / D21 00368 D22 = W( K-1, KW-1 ) / D21 00369 T = CONE / ( D11*D22-CONE ) 00370 D21 = T / D21 00371 DO 20 J = 1, K - 2 00372 A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) ) 00373 A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) ) 00374 20 CONTINUE 00375 END IF 00376 * 00377 * Copy D(k) to A 00378 * 00379 A( K-1, K-1 ) = W( K-1, KW-1 ) 00380 A( K-1, K ) = W( K-1, KW ) 00381 A( K, K ) = W( K, KW ) 00382 END IF 00383 END IF 00384 * 00385 * Store details of the interchanges in IPIV 00386 * 00387 IF( KSTEP.EQ.1 ) THEN 00388 IPIV( K ) = KP 00389 ELSE 00390 IPIV( K ) = -KP 00391 IPIV( K-1 ) = -KP 00392 END IF 00393 * 00394 * Decrease K and return to the start of the main loop 00395 * 00396 K = K - KSTEP 00397 GO TO 10 00398 * 00399 30 CONTINUE 00400 * 00401 * Update the upper triangle of A11 (= A(1:k,1:k)) as 00402 * 00403 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T 00404 * 00405 * computing blocks of NB columns at a time 00406 * 00407 DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB 00408 JB = MIN( NB, K-J+1 ) 00409 * 00410 * Update the upper triangle of the diagonal block 00411 * 00412 DO 40 JJ = J, J + JB - 1 00413 CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE, 00414 $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE, 00415 $ A( J, JJ ), 1 ) 00416 40 CONTINUE 00417 * 00418 * Update the rectangular superdiagonal block 00419 * 00420 CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K, 00421 $ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, 00422 $ CONE, A( 1, J ), LDA ) 00423 50 CONTINUE 00424 * 00425 * Put U12 in standard form by partially undoing the interchanges 00426 * in columns k+1:n 00427 * 00428 J = K + 1 00429 60 CONTINUE 00430 JJ = J 00431 JP = IPIV( J ) 00432 IF( JP.LT.0 ) THEN 00433 JP = -JP 00434 J = J + 1 00435 END IF 00436 J = J + 1 00437 IF( JP.NE.JJ .AND. J.LE.N ) 00438 $ CALL ZSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA ) 00439 IF( J.LE.N ) 00440 $ GO TO 60 00441 * 00442 * Set KB to the number of columns factorized 00443 * 00444 KB = N - K 00445 * 00446 ELSE 00447 * 00448 * Factorize the leading columns of A using the lower triangle 00449 * of A and working forwards, and compute the matrix W = L21*D 00450 * for use in updating A22 00451 * 00452 * K is the main loop index, increasing from 1 in steps of 1 or 2 00453 * 00454 K = 1 00455 70 CONTINUE 00456 * 00457 * Exit from loop 00458 * 00459 IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N ) 00460 $ GO TO 90 00461 * 00462 * Copy column K of A to column K of W and update it 00463 * 00464 CALL ZCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 ) 00465 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA, 00466 $ W( K, 1 ), LDW, CONE, W( K, K ), 1 ) 00467 * 00468 KSTEP = 1 00469 * 00470 * Determine rows and columns to be interchanged and whether 00471 * a 1-by-1 or 2-by-2 pivot block will be used 00472 * 00473 ABSAKK = CABS1( W( K, K ) ) 00474 * 00475 * IMAX is the row-index of the largest off-diagonal element in 00476 * column K, and COLMAX is its absolute value 00477 * 00478 IF( K.LT.N ) THEN 00479 IMAX = K + IZAMAX( N-K, W( K+1, K ), 1 ) 00480 COLMAX = CABS1( W( IMAX, K ) ) 00481 ELSE 00482 COLMAX = ZERO 00483 END IF 00484 * 00485 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 00486 * 00487 * Column K is zero: set INFO and continue 00488 * 00489 IF( INFO.EQ.0 ) 00490 $ INFO = K 00491 KP = K 00492 ELSE 00493 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 00494 * 00495 * no interchange, use 1-by-1 pivot block 00496 * 00497 KP = K 00498 ELSE 00499 * 00500 * Copy column IMAX to column K+1 of W and update it 00501 * 00502 CALL ZCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 ) 00503 CALL ZCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ), 00504 $ 1 ) 00505 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), 00506 $ LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ), 00507 $ 1 ) 00508 * 00509 * JMAX is the column-index of the largest off-diagonal 00510 * element in row IMAX, and ROWMAX is its absolute value 00511 * 00512 JMAX = K - 1 + IZAMAX( IMAX-K, W( K, K+1 ), 1 ) 00513 ROWMAX = CABS1( W( JMAX, K+1 ) ) 00514 IF( IMAX.LT.N ) THEN 00515 JMAX = IMAX + IZAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 ) 00516 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) ) 00517 END IF 00518 * 00519 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 00520 * 00521 * no interchange, use 1-by-1 pivot block 00522 * 00523 KP = K 00524 ELSE IF( CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN 00525 * 00526 * interchange rows and columns K and IMAX, use 1-by-1 00527 * pivot block 00528 * 00529 KP = IMAX 00530 * 00531 * copy column K+1 of W to column K 00532 * 00533 CALL ZCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 ) 00534 ELSE 00535 * 00536 * interchange rows and columns K+1 and IMAX, use 2-by-2 00537 * pivot block 00538 * 00539 KP = IMAX 00540 KSTEP = 2 00541 END IF 00542 END IF 00543 * 00544 KK = K + KSTEP - 1 00545 * 00546 * Updated column KP is already stored in column KK of W 00547 * 00548 IF( KP.NE.KK ) THEN 00549 * 00550 * Copy non-updated column KK to column KP 00551 * 00552 A( KP, K ) = A( KK, K ) 00553 CALL ZCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA ) 00554 CALL ZCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 ) 00555 * 00556 * Interchange rows KK and KP in first KK columns of A and W 00557 * 00558 CALL ZSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA ) 00559 CALL ZSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW ) 00560 END IF 00561 * 00562 IF( KSTEP.EQ.1 ) THEN 00563 * 00564 * 1-by-1 pivot block D(k): column k of W now holds 00565 * 00566 * W(k) = L(k)*D(k) 00567 * 00568 * where L(k) is the k-th column of L 00569 * 00570 * Store L(k) in column k of A 00571 * 00572 CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 ) 00573 IF( K.LT.N ) THEN 00574 R1 = CONE / A( K, K ) 00575 CALL ZSCAL( N-K, R1, A( K+1, K ), 1 ) 00576 END IF 00577 ELSE 00578 * 00579 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold 00580 * 00581 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) 00582 * 00583 * where L(k) and L(k+1) are the k-th and (k+1)-th columns 00584 * of L 00585 * 00586 IF( K.LT.N-1 ) THEN 00587 * 00588 * Store L(k) and L(k+1) in columns k and k+1 of A 00589 * 00590 D21 = W( K+1, K ) 00591 D11 = W( K+1, K+1 ) / D21 00592 D22 = W( K, K ) / D21 00593 T = CONE / ( D11*D22-CONE ) 00594 D21 = T / D21 00595 DO 80 J = K + 2, N 00596 A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) ) 00597 A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) ) 00598 80 CONTINUE 00599 END IF 00600 * 00601 * Copy D(k) to A 00602 * 00603 A( K, K ) = W( K, K ) 00604 A( K+1, K ) = W( K+1, K ) 00605 A( K+1, K+1 ) = W( K+1, K+1 ) 00606 END IF 00607 END IF 00608 * 00609 * Store details of the interchanges in IPIV 00610 * 00611 IF( KSTEP.EQ.1 ) THEN 00612 IPIV( K ) = KP 00613 ELSE 00614 IPIV( K ) = -KP 00615 IPIV( K+1 ) = -KP 00616 END IF 00617 * 00618 * Increase K and return to the start of the main loop 00619 * 00620 K = K + KSTEP 00621 GO TO 70 00622 * 00623 90 CONTINUE 00624 * 00625 * Update the lower triangle of A22 (= A(k:n,k:n)) as 00626 * 00627 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T 00628 * 00629 * computing blocks of NB columns at a time 00630 * 00631 DO 110 J = K, N, NB 00632 JB = MIN( NB, N-J+1 ) 00633 * 00634 * Update the lower triangle of the diagonal block 00635 * 00636 DO 100 JJ = J, J + JB - 1 00637 CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE, 00638 $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE, 00639 $ A( JJ, JJ ), 1 ) 00640 100 CONTINUE 00641 * 00642 * Update the rectangular subdiagonal block 00643 * 00644 IF( J+JB.LE.N ) 00645 $ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB, 00646 $ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ), 00647 $ LDW, CONE, A( J+JB, J ), LDA ) 00648 110 CONTINUE 00649 * 00650 * Put L21 in standard form by partially undoing the interchanges 00651 * in columns 1:k-1 00652 * 00653 J = K - 1 00654 120 CONTINUE 00655 JJ = J 00656 JP = IPIV( J ) 00657 IF( JP.LT.0 ) THEN 00658 JP = -JP 00659 J = J - 1 00660 END IF 00661 J = J - 1 00662 IF( JP.NE.JJ .AND. J.GE.1 ) 00663 $ CALL ZSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA ) 00664 IF( J.GE.1 ) 00665 $ GO TO 120 00666 * 00667 * Set KB to the number of columns factorized 00668 * 00669 KB = K - 1 00670 * 00671 END IF 00672 RETURN 00673 * 00674 * End of ZLASYF 00675 * 00676 END