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