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