|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CSYTF2( UPLO, N, A, LDA, IPIV, 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, LDA, N 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IPIV( * ) 00014 COMPLEX A( LDA, * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * CSYTF2 computes the factorization of a complex symmetric matrix A 00021 * using the Bunch-Kaufman diagonal pivoting method: 00022 * 00023 * A = U*D*U' or A = L*D*L' 00024 * 00025 * where U (or L) is a product of permutation and unit upper (lower) 00026 * triangular matrices, U' is the transpose of U, and D is symmetric and 00027 * block diagonal with 1-by-1 and 2-by-2 diagonal blocks. 00028 * 00029 * This is the unblocked version of the algorithm, calling Level 2 BLAS. 00030 * 00031 * Arguments 00032 * ========= 00033 * 00034 * UPLO (input) CHARACTER*1 00035 * Specifies whether the upper or lower triangular part of the 00036 * symmetric matrix A is stored: 00037 * = 'U': Upper triangular 00038 * = 'L': Lower triangular 00039 * 00040 * N (input) INTEGER 00041 * The order of the matrix A. N >= 0. 00042 * 00043 * A (input/output) COMPLEX array, dimension (LDA,N) 00044 * On entry, the symmetric matrix A. If UPLO = 'U', the leading 00045 * n-by-n upper triangular part of A contains the upper 00046 * triangular part of the matrix A, and the strictly lower 00047 * triangular part of A is not referenced. If UPLO = 'L', the 00048 * leading n-by-n lower triangular part of A contains the lower 00049 * triangular part of the matrix A, and the strictly upper 00050 * triangular part of A is not referenced. 00051 * 00052 * On exit, the block diagonal matrix D and the multipliers used 00053 * to obtain the factor U or L (see below for further details). 00054 * 00055 * LDA (input) INTEGER 00056 * The leading dimension of the array A. LDA >= max(1,N). 00057 * 00058 * IPIV (output) INTEGER array, dimension (N) 00059 * Details of the interchanges and the block structure of D. 00060 * If IPIV(k) > 0, then rows and columns k and IPIV(k) were 00061 * interchanged and D(k,k) is a 1-by-1 diagonal block. 00062 * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 00063 * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 00064 * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 00065 * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 00066 * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 00067 * 00068 * INFO (output) INTEGER 00069 * = 0: successful exit 00070 * < 0: if INFO = -k, the k-th argument had an illegal value 00071 * > 0: if INFO = k, D(k,k) is exactly zero. The factorization 00072 * has been completed, but the block diagonal matrix D is 00073 * exactly singular, and division by zero will occur if it 00074 * is used to solve a system of equations. 00075 * 00076 * Further Details 00077 * =============== 00078 * 00079 * 09-29-06 - patch from 00080 * Bobby Cheng, MathWorks 00081 * 00082 * Replace l.209 and l.377 00083 * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN 00084 * by 00085 * IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN 00086 * 00087 * 1-96 - Based on modifications by J. Lewis, Boeing Computer Services 00088 * Company 00089 * 00090 * If UPLO = 'U', then A = U*D*U', where 00091 * U = P(n)*U(n)* ... *P(k)U(k)* ..., 00092 * i.e., U is a product of terms P(k)*U(k), where k decreases from n to 00093 * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 00094 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 00095 * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such 00096 * that if the diagonal block D(k) is of order s (s = 1 or 2), then 00097 * 00098 * ( I v 0 ) k-s 00099 * U(k) = ( 0 I 0 ) s 00100 * ( 0 0 I ) n-k 00101 * k-s s n-k 00102 * 00103 * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). 00104 * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), 00105 * and A(k,k), and v overwrites A(1:k-2,k-1:k). 00106 * 00107 * If UPLO = 'L', then A = L*D*L', where 00108 * L = P(1)*L(1)* ... *P(k)*L(k)* ..., 00109 * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to 00110 * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 00111 * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as 00112 * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such 00113 * that if the diagonal block D(k) is of order s (s = 1 or 2), then 00114 * 00115 * ( I 0 0 ) k-1 00116 * L(k) = ( 0 I 0 ) s 00117 * ( 0 v I ) n-k-s+1 00118 * k-1 s n-k-s+1 00119 * 00120 * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). 00121 * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), 00122 * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). 00123 * 00124 * ===================================================================== 00125 * 00126 * .. Parameters .. 00127 REAL ZERO, ONE 00128 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00129 REAL EIGHT, SEVTEN 00130 PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 ) 00131 COMPLEX CONE 00132 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00133 * .. 00134 * .. Local Scalars .. 00135 LOGICAL UPPER 00136 INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP 00137 REAL ABSAKK, ALPHA, COLMAX, ROWMAX 00138 COMPLEX D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, Z 00139 * .. 00140 * .. External Functions .. 00141 LOGICAL LSAME, SISNAN 00142 INTEGER ICAMAX 00143 EXTERNAL LSAME, ICAMAX, SISNAN 00144 * .. 00145 * .. External Subroutines .. 00146 EXTERNAL CSCAL, CSWAP, CSYR, XERBLA 00147 * .. 00148 * .. Intrinsic Functions .. 00149 INTRINSIC ABS, AIMAG, MAX, REAL, SQRT 00150 * .. 00151 * .. Statement Functions .. 00152 REAL CABS1 00153 * .. 00154 * .. Statement Function definitions .. 00155 CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) ) 00156 * .. 00157 * .. Executable Statements .. 00158 * 00159 * Test the input parameters. 00160 * 00161 INFO = 0 00162 UPPER = LSAME( UPLO, 'U' ) 00163 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00164 INFO = -1 00165 ELSE IF( N.LT.0 ) THEN 00166 INFO = -2 00167 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00168 INFO = -4 00169 END IF 00170 IF( INFO.NE.0 ) THEN 00171 CALL XERBLA( 'CSYTF2', -INFO ) 00172 RETURN 00173 END IF 00174 * 00175 * Initialize ALPHA for use in choosing pivot block size. 00176 * 00177 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT 00178 * 00179 IF( UPPER ) THEN 00180 * 00181 * Factorize A as U*D*U' using the upper triangle of A 00182 * 00183 * K is the main loop index, decreasing from N to 1 in steps of 00184 * 1 or 2 00185 * 00186 K = N 00187 10 CONTINUE 00188 * 00189 * If K < 1, exit from loop 00190 * 00191 IF( K.LT.1 ) 00192 $ GO TO 70 00193 KSTEP = 1 00194 * 00195 * Determine rows and columns to be interchanged and whether 00196 * a 1-by-1 or 2-by-2 pivot block will be used 00197 * 00198 ABSAKK = CABS1( A( K, K ) ) 00199 * 00200 * IMAX is the row-index of the largest off-diagonal element in 00201 * column K, and COLMAX is its absolute value 00202 * 00203 IF( K.GT.1 ) THEN 00204 IMAX = ICAMAX( K-1, A( 1, K ), 1 ) 00205 COLMAX = CABS1( A( IMAX, K ) ) 00206 ELSE 00207 COLMAX = ZERO 00208 END IF 00209 * 00210 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. SISNAN(ABSAKK) ) THEN 00211 * 00212 * Column K is zero or NaN: set INFO and continue 00213 * 00214 IF( INFO.EQ.0 ) 00215 $ INFO = K 00216 KP = K 00217 ELSE 00218 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 00219 * 00220 * no interchange, use 1-by-1 pivot block 00221 * 00222 KP = K 00223 ELSE 00224 * 00225 * JMAX is the column-index of the largest off-diagonal 00226 * element in row IMAX, and ROWMAX is its absolute value 00227 * 00228 JMAX = IMAX + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA ) 00229 ROWMAX = CABS1( A( IMAX, JMAX ) ) 00230 IF( IMAX.GT.1 ) THEN 00231 JMAX = ICAMAX( IMAX-1, A( 1, IMAX ), 1 ) 00232 ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) ) 00233 END IF 00234 * 00235 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 00236 * 00237 * no interchange, use 1-by-1 pivot block 00238 * 00239 KP = K 00240 ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN 00241 * 00242 * interchange rows and columns K and IMAX, use 1-by-1 00243 * pivot block 00244 * 00245 KP = IMAX 00246 ELSE 00247 * 00248 * interchange rows and columns K-1 and IMAX, use 2-by-2 00249 * pivot block 00250 * 00251 KP = IMAX 00252 KSTEP = 2 00253 END IF 00254 END IF 00255 * 00256 KK = K - KSTEP + 1 00257 IF( KP.NE.KK ) THEN 00258 * 00259 * Interchange rows and columns KK and KP in the leading 00260 * submatrix A(1:k,1:k) 00261 * 00262 CALL CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 ) 00263 CALL CSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ), 00264 $ LDA ) 00265 T = A( KK, KK ) 00266 A( KK, KK ) = A( KP, KP ) 00267 A( KP, KP ) = T 00268 IF( KSTEP.EQ.2 ) THEN 00269 T = A( K-1, K ) 00270 A( K-1, K ) = A( KP, K ) 00271 A( KP, K ) = T 00272 END IF 00273 END IF 00274 * 00275 * Update the leading submatrix 00276 * 00277 IF( KSTEP.EQ.1 ) THEN 00278 * 00279 * 1-by-1 pivot block D(k): column k now holds 00280 * 00281 * W(k) = U(k)*D(k) 00282 * 00283 * where U(k) is the k-th column of U 00284 * 00285 * Perform a rank-1 update of A(1:k-1,1:k-1) as 00286 * 00287 * A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)' 00288 * 00289 R1 = CONE / A( K, K ) 00290 CALL CSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA ) 00291 * 00292 * Store U(k) in column k 00293 * 00294 CALL CSCAL( K-1, R1, A( 1, K ), 1 ) 00295 ELSE 00296 * 00297 * 2-by-2 pivot block D(k): columns k and k-1 now hold 00298 * 00299 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) 00300 * 00301 * where U(k) and U(k-1) are the k-th and (k-1)-th columns 00302 * of U 00303 * 00304 * Perform a rank-2 update of A(1:k-2,1:k-2) as 00305 * 00306 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )' 00307 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )' 00308 * 00309 IF( K.GT.2 ) THEN 00310 * 00311 D12 = A( K-1, K ) 00312 D22 = A( K-1, K-1 ) / D12 00313 D11 = A( K, K ) / D12 00314 T = CONE / ( D11*D22-CONE ) 00315 D12 = T / D12 00316 * 00317 DO 30 J = K - 2, 1, -1 00318 WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) ) 00319 WK = D12*( D22*A( J, K )-A( J, K-1 ) ) 00320 DO 20 I = J, 1, -1 00321 A( I, J ) = A( I, J ) - A( I, K )*WK - 00322 $ A( I, K-1 )*WKM1 00323 20 CONTINUE 00324 A( J, K ) = WK 00325 A( J, K-1 ) = WKM1 00326 30 CONTINUE 00327 * 00328 END IF 00329 * 00330 END IF 00331 END IF 00332 * 00333 * Store details of the interchanges in IPIV 00334 * 00335 IF( KSTEP.EQ.1 ) THEN 00336 IPIV( K ) = KP 00337 ELSE 00338 IPIV( K ) = -KP 00339 IPIV( K-1 ) = -KP 00340 END IF 00341 * 00342 * Decrease K and return to the start of the main loop 00343 * 00344 K = K - KSTEP 00345 GO TO 10 00346 * 00347 ELSE 00348 * 00349 * Factorize A as L*D*L' using the lower triangle of A 00350 * 00351 * K is the main loop index, increasing from 1 to N in steps of 00352 * 1 or 2 00353 * 00354 K = 1 00355 40 CONTINUE 00356 * 00357 * If K > N, exit from loop 00358 * 00359 IF( K.GT.N ) 00360 $ GO TO 70 00361 KSTEP = 1 00362 * 00363 * Determine rows and columns to be interchanged and whether 00364 * a 1-by-1 or 2-by-2 pivot block will be used 00365 * 00366 ABSAKK = CABS1( A( K, K ) ) 00367 * 00368 * IMAX is the row-index of the largest off-diagonal element in 00369 * column K, and COLMAX is its absolute value 00370 * 00371 IF( K.LT.N ) THEN 00372 IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 ) 00373 COLMAX = CABS1( A( IMAX, K ) ) 00374 ELSE 00375 COLMAX = ZERO 00376 END IF 00377 * 00378 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO .OR. SISNAN(ABSAKK) ) THEN 00379 * 00380 * Column K is zero or NaN: set INFO and continue 00381 * 00382 IF( INFO.EQ.0 ) 00383 $ INFO = K 00384 KP = K 00385 ELSE 00386 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN 00387 * 00388 * no interchange, use 1-by-1 pivot block 00389 * 00390 KP = K 00391 ELSE 00392 * 00393 * JMAX is the column-index of the largest off-diagonal 00394 * element in row IMAX, and ROWMAX is its absolute value 00395 * 00396 JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA ) 00397 ROWMAX = CABS1( A( IMAX, JMAX ) ) 00398 IF( IMAX.LT.N ) THEN 00399 JMAX = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 ) 00400 ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) ) 00401 END IF 00402 * 00403 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN 00404 * 00405 * no interchange, use 1-by-1 pivot block 00406 * 00407 KP = K 00408 ELSE IF( CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN 00409 * 00410 * interchange rows and columns K and IMAX, use 1-by-1 00411 * pivot block 00412 * 00413 KP = IMAX 00414 ELSE 00415 * 00416 * interchange rows and columns K+1 and IMAX, use 2-by-2 00417 * pivot block 00418 * 00419 KP = IMAX 00420 KSTEP = 2 00421 END IF 00422 END IF 00423 * 00424 KK = K + KSTEP - 1 00425 IF( KP.NE.KK ) THEN 00426 * 00427 * Interchange rows and columns KK and KP in the trailing 00428 * submatrix A(k:n,k:n) 00429 * 00430 IF( KP.LT.N ) 00431 $ CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 ) 00432 CALL CSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ), 00433 $ LDA ) 00434 T = A( KK, KK ) 00435 A( KK, KK ) = A( KP, KP ) 00436 A( KP, KP ) = T 00437 IF( KSTEP.EQ.2 ) THEN 00438 T = A( K+1, K ) 00439 A( K+1, K ) = A( KP, K ) 00440 A( KP, K ) = T 00441 END IF 00442 END IF 00443 * 00444 * Update the trailing submatrix 00445 * 00446 IF( KSTEP.EQ.1 ) THEN 00447 * 00448 * 1-by-1 pivot block D(k): column k now holds 00449 * 00450 * W(k) = L(k)*D(k) 00451 * 00452 * where L(k) is the k-th column of L 00453 * 00454 IF( K.LT.N ) THEN 00455 * 00456 * Perform a rank-1 update of A(k+1:n,k+1:n) as 00457 * 00458 * A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)' 00459 * 00460 R1 = CONE / A( K, K ) 00461 CALL CSYR( UPLO, N-K, -R1, A( K+1, K ), 1, 00462 $ A( K+1, K+1 ), LDA ) 00463 * 00464 * Store L(k) in column K 00465 * 00466 CALL CSCAL( N-K, R1, A( K+1, K ), 1 ) 00467 END IF 00468 ELSE 00469 * 00470 * 2-by-2 pivot block D(k) 00471 * 00472 IF( K.LT.N-1 ) THEN 00473 * 00474 * Perform a rank-2 update of A(k+2:n,k+2:n) as 00475 * 00476 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )' 00477 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )' 00478 * 00479 * where L(k) and L(k+1) are the k-th and (k+1)-th 00480 * columns of L 00481 * 00482 D21 = A( K+1, K ) 00483 D11 = A( K+1, K+1 ) / D21 00484 D22 = A( K, K ) / D21 00485 T = CONE / ( D11*D22-CONE ) 00486 D21 = T / D21 00487 * 00488 DO 60 J = K + 2, N 00489 WK = D21*( D11*A( J, K )-A( J, K+1 ) ) 00490 WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) ) 00491 DO 50 I = J, N 00492 A( I, J ) = A( I, J ) - A( I, K )*WK - 00493 $ A( I, K+1 )*WKP1 00494 50 CONTINUE 00495 A( J, K ) = WK 00496 A( J, K+1 ) = WKP1 00497 60 CONTINUE 00498 END IF 00499 END IF 00500 END IF 00501 * 00502 * Store details of the interchanges in IPIV 00503 * 00504 IF( KSTEP.EQ.1 ) THEN 00505 IPIV( K ) = KP 00506 ELSE 00507 IPIV( K ) = -KP 00508 IPIV( K+1 ) = -KP 00509 END IF 00510 * 00511 * Increase K and return to the start of the main loop 00512 * 00513 K = K + KSTEP 00514 GO TO 40 00515 * 00516 END IF 00517 * 00518 70 CONTINUE 00519 RETURN 00520 * 00521 * End of CSYTF2 00522 * 00523 END
1.7.3 
