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