![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CSYTRI2X 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CSYTRI2X + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytri2x.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytri2x.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytri2x.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, LDA, N, NB 00026 * .. 00027 * .. Array Arguments .. 00028 * INTEGER IPIV( * ) 00029 * COMPLEX A( LDA, * ), WORK( N+NB+1,* ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> CSYTRI2X computes the inverse of a real symmetric indefinite matrix 00039 *> A using the factorization A = U*D*U**T or A = L*D*L**T computed by 00040 *> CSYTRF. 00041 *> \endverbatim 00042 * 00043 * Arguments: 00044 * ========== 00045 * 00046 *> \param[in] UPLO 00047 *> \verbatim 00048 *> UPLO is CHARACTER*1 00049 *> Specifies whether the details of the factorization are stored 00050 *> as an upper or lower triangular matrix. 00051 *> = 'U': Upper triangular, form is A = U*D*U**T; 00052 *> = 'L': Lower triangular, form is A = L*D*L**T. 00053 *> \endverbatim 00054 *> 00055 *> \param[in] N 00056 *> \verbatim 00057 *> N is INTEGER 00058 *> The order of the matrix A. N >= 0. 00059 *> \endverbatim 00060 *> 00061 *> \param[in,out] A 00062 *> \verbatim 00063 *> A is COMPLEX array, dimension (LDA,N) 00064 *> On entry, the NNB diagonal matrix D and the multipliers 00065 *> used to obtain the factor U or L as computed by CSYTRF. 00066 *> 00067 *> On exit, if INFO = 0, the (symmetric) inverse of the original 00068 *> matrix. If UPLO = 'U', the upper triangular part of the 00069 *> inverse is formed and the part of A below the diagonal is not 00070 *> referenced; if UPLO = 'L' the lower triangular part of the 00071 *> inverse is formed and the part of A above the diagonal is 00072 *> not referenced. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] LDA 00076 *> \verbatim 00077 *> LDA is INTEGER 00078 *> The leading dimension of the array A. LDA >= max(1,N). 00079 *> \endverbatim 00080 *> 00081 *> \param[in] IPIV 00082 *> \verbatim 00083 *> IPIV is INTEGER array, dimension (N) 00084 *> Details of the interchanges and the NNB structure of D 00085 *> as determined by CSYTRF. 00086 *> \endverbatim 00087 *> 00088 *> \param[out] WORK 00089 *> \verbatim 00090 *> WORK is COMPLEX array, dimension (N+NNB+1,NNB+3) 00091 *> \endverbatim 00092 *> 00093 *> \param[in] NB 00094 *> \verbatim 00095 *> NB is INTEGER 00096 *> Block size 00097 *> \endverbatim 00098 *> 00099 *> \param[out] INFO 00100 *> \verbatim 00101 *> INFO is INTEGER 00102 *> = 0: successful exit 00103 *> < 0: if INFO = -i, the i-th argument had an illegal value 00104 *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its 00105 *> inverse could not be computed. 00106 *> \endverbatim 00107 * 00108 * Authors: 00109 * ======== 00110 * 00111 *> \author Univ. of Tennessee 00112 *> \author Univ. of California Berkeley 00113 *> \author Univ. of Colorado Denver 00114 *> \author NAG Ltd. 00115 * 00116 *> \date November 2011 00117 * 00118 *> \ingroup complexSYcomputational 00119 * 00120 * ===================================================================== 00121 SUBROUTINE CSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) 00122 * 00123 * -- LAPACK computational routine (version 3.4.0) -- 00124 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00125 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00126 * November 2011 00127 * 00128 * .. Scalar Arguments .. 00129 CHARACTER UPLO 00130 INTEGER INFO, LDA, N, NB 00131 * .. 00132 * .. Array Arguments .. 00133 INTEGER IPIV( * ) 00134 COMPLEX A( LDA, * ), WORK( N+NB+1,* ) 00135 * .. 00136 * 00137 * ===================================================================== 00138 * 00139 * .. Parameters .. 00140 COMPLEX ONE, ZERO 00141 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), 00142 $ ZERO = ( 0.0E+0, 0.0E+0 ) ) 00143 * .. 00144 * .. Local Scalars .. 00145 LOGICAL UPPER 00146 INTEGER I, IINFO, IP, K, CUT, NNB 00147 INTEGER COUNT 00148 INTEGER J, U11, INVD 00149 00150 COMPLEX AK, AKKP1, AKP1, D, T 00151 COMPLEX U01_I_J, U01_IP1_J 00152 COMPLEX U11_I_J, U11_IP1_J 00153 * .. 00154 * .. External Functions .. 00155 LOGICAL LSAME 00156 EXTERNAL LSAME 00157 * .. 00158 * .. External Subroutines .. 00159 EXTERNAL CSYCONV, XERBLA, CTRTRI 00160 EXTERNAL CGEMM, CTRMM, CSYSWAPR 00161 * .. 00162 * .. Intrinsic Functions .. 00163 INTRINSIC MAX 00164 * .. 00165 * .. Executable Statements .. 00166 * 00167 * Test the input parameters. 00168 * 00169 INFO = 0 00170 UPPER = LSAME( UPLO, 'U' ) 00171 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00172 INFO = -1 00173 ELSE IF( N.LT.0 ) THEN 00174 INFO = -2 00175 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00176 INFO = -4 00177 END IF 00178 * 00179 * Quick return if possible 00180 * 00181 * 00182 IF( INFO.NE.0 ) THEN 00183 CALL XERBLA( 'CSYTRI2X', -INFO ) 00184 RETURN 00185 END IF 00186 IF( N.EQ.0 ) 00187 $ RETURN 00188 * 00189 * Convert A 00190 * Workspace got Non-diag elements of D 00191 * 00192 CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 00193 * 00194 * Check that the diagonal matrix D is nonsingular. 00195 * 00196 IF( UPPER ) THEN 00197 * 00198 * Upper triangular storage: examine D from bottom to top 00199 * 00200 DO INFO = N, 1, -1 00201 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 00202 $ RETURN 00203 END DO 00204 ELSE 00205 * 00206 * Lower triangular storage: examine D from top to bottom. 00207 * 00208 DO INFO = 1, N 00209 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 00210 $ RETURN 00211 END DO 00212 END IF 00213 INFO = 0 00214 * 00215 * Splitting Workspace 00216 * U01 is a block (N,NB+1) 00217 * The first element of U01 is in WORK(1,1) 00218 * U11 is a block (NB+1,NB+1) 00219 * The first element of U11 is in WORK(N+1,1) 00220 U11 = N 00221 * INVD is a block (N,2) 00222 * The first element of INVD is in WORK(1,INVD) 00223 INVD = NB+2 00224 00225 IF( UPPER ) THEN 00226 * 00227 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T. 00228 * 00229 CALL CTRTRI( UPLO, 'U', N, A, LDA, INFO ) 00230 * 00231 * inv(D) and inv(D)*inv(U) 00232 * 00233 K=1 00234 DO WHILE ( K .LE. N ) 00235 IF( IPIV( K ).GT.0 ) THEN 00236 * 1 x 1 diagonal NNB 00237 WORK(K,INVD) = ONE / A( K, K ) 00238 WORK(K,INVD+1) = 0 00239 K=K+1 00240 ELSE 00241 * 2 x 2 diagonal NNB 00242 T = WORK(K+1,1) 00243 AK = A( K, K ) / T 00244 AKP1 = A( K+1, K+1 ) / T 00245 AKKP1 = WORK(K+1,1) / T 00246 D = T*( AK*AKP1-ONE ) 00247 WORK(K,INVD) = AKP1 / D 00248 WORK(K+1,INVD+1) = AK / D 00249 WORK(K,INVD+1) = -AKKP1 / D 00250 WORK(K+1,INVD) = -AKKP1 / D 00251 K=K+2 00252 END IF 00253 END DO 00254 * 00255 * inv(U**T) = (inv(U))**T 00256 * 00257 * inv(U**T)*inv(D)*inv(U) 00258 * 00259 CUT=N 00260 DO WHILE (CUT .GT. 0) 00261 NNB=NB 00262 IF (CUT .LE. NNB) THEN 00263 NNB=CUT 00264 ELSE 00265 COUNT = 0 00266 * count negative elements, 00267 DO I=CUT+1-NNB,CUT 00268 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 00269 END DO 00270 * need a even number for a clear cut 00271 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 00272 END IF 00273 00274 CUT=CUT-NNB 00275 * 00276 * U01 Block 00277 * 00278 DO I=1,CUT 00279 DO J=1,NNB 00280 WORK(I,J)=A(I,CUT+J) 00281 END DO 00282 END DO 00283 * 00284 * U11 Block 00285 * 00286 DO I=1,NNB 00287 WORK(U11+I,I)=ONE 00288 DO J=1,I-1 00289 WORK(U11+I,J)=ZERO 00290 END DO 00291 DO J=I+1,NNB 00292 WORK(U11+I,J)=A(CUT+I,CUT+J) 00293 END DO 00294 END DO 00295 * 00296 * invD*U01 00297 * 00298 I=1 00299 DO WHILE (I .LE. CUT) 00300 IF (IPIV(I) > 0) THEN 00301 DO J=1,NNB 00302 WORK(I,J)=WORK(I,INVD)*WORK(I,J) 00303 END DO 00304 I=I+1 00305 ELSE 00306 DO J=1,NNB 00307 U01_I_J = WORK(I,J) 00308 U01_IP1_J = WORK(I+1,J) 00309 WORK(I,J)=WORK(I,INVD)*U01_I_J+ 00310 $ WORK(I,INVD+1)*U01_IP1_J 00311 WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+ 00312 $ WORK(I+1,INVD+1)*U01_IP1_J 00313 END DO 00314 I=I+2 00315 END IF 00316 END DO 00317 * 00318 * invD1*U11 00319 * 00320 I=1 00321 DO WHILE (I .LE. NNB) 00322 IF (IPIV(CUT+I) > 0) THEN 00323 DO J=I,NNB 00324 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 00325 END DO 00326 I=I+1 00327 ELSE 00328 DO J=I,NNB 00329 U11_I_J = WORK(U11+I,J) 00330 U11_IP1_J = WORK(U11+I+1,J) 00331 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 00332 $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J) 00333 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+ 00334 $ WORK(CUT+I+1,INVD+1)*U11_IP1_J 00335 END DO 00336 I=I+2 00337 END IF 00338 END DO 00339 * 00340 * U11**T*invD1*U11->U11 00341 * 00342 CALL CTRMM('L','U','T','U',NNB, NNB, 00343 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 00344 * 00345 DO I=1,NNB 00346 DO J=I,NNB 00347 A(CUT+I,CUT+J)=WORK(U11+I,J) 00348 END DO 00349 END DO 00350 * 00351 * U01**T*invD*U01->A(CUT+I,CUT+J) 00352 * 00353 CALL CGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA, 00354 $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 00355 * 00356 * U11 = U11**T*invD1*U11 + U01**T*invD*U01 00357 * 00358 DO I=1,NNB 00359 DO J=I,NNB 00360 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 00361 END DO 00362 END DO 00363 * 00364 * U01 = U00**T*invD0*U01 00365 * 00366 CALL CTRMM('L',UPLO,'T','U',CUT, NNB, 00367 $ ONE,A,LDA,WORK,N+NB+1) 00368 00369 * 00370 * Update U01 00371 * 00372 DO I=1,CUT 00373 DO J=1,NNB 00374 A(I,CUT+J)=WORK(I,J) 00375 END DO 00376 END DO 00377 * 00378 * Next Block 00379 * 00380 END DO 00381 * 00382 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T 00383 * 00384 I=1 00385 DO WHILE ( I .LE. N ) 00386 IF( IPIV(I) .GT. 0 ) THEN 00387 IP=IPIV(I) 00388 IF (I .LT. IP) CALL CSYSWAPR( UPLO, N, A, LDA, I ,IP ) 00389 IF (I .GT. IP) CALL CSYSWAPR( UPLO, N, A, LDA, IP ,I ) 00390 ELSE 00391 IP=-IPIV(I) 00392 I=I+1 00393 IF ( (I-1) .LT. IP) 00394 $ CALL CSYSWAPR( UPLO, N, A, LDA, I-1 ,IP ) 00395 IF ( (I-1) .GT. IP) 00396 $ CALL CSYSWAPR( UPLO, N, A, LDA, IP ,I-1 ) 00397 ENDIF 00398 I=I+1 00399 END DO 00400 ELSE 00401 * 00402 * LOWER... 00403 * 00404 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T. 00405 * 00406 CALL CTRTRI( UPLO, 'U', N, A, LDA, INFO ) 00407 * 00408 * inv(D) and inv(D)*inv(U) 00409 * 00410 K=N 00411 DO WHILE ( K .GE. 1 ) 00412 IF( IPIV( K ).GT.0 ) THEN 00413 * 1 x 1 diagonal NNB 00414 WORK(K,INVD) = ONE / A( K, K ) 00415 WORK(K,INVD+1) = 0 00416 K=K-1 00417 ELSE 00418 * 2 x 2 diagonal NNB 00419 T = WORK(K-1,1) 00420 AK = A( K-1, K-1 ) / T 00421 AKP1 = A( K, K ) / T 00422 AKKP1 = WORK(K-1,1) / T 00423 D = T*( AK*AKP1-ONE ) 00424 WORK(K-1,INVD) = AKP1 / D 00425 WORK(K,INVD) = AK / D 00426 WORK(K,INVD+1) = -AKKP1 / D 00427 WORK(K-1,INVD+1) = -AKKP1 / D 00428 K=K-2 00429 END IF 00430 END DO 00431 * 00432 * inv(U**T) = (inv(U))**T 00433 * 00434 * inv(U**T)*inv(D)*inv(U) 00435 * 00436 CUT=0 00437 DO WHILE (CUT .LT. N) 00438 NNB=NB 00439 IF (CUT + NNB .GE. N) THEN 00440 NNB=N-CUT 00441 ELSE 00442 COUNT = 0 00443 * count negative elements, 00444 DO I=CUT+1,CUT+NNB 00445 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 00446 END DO 00447 * need a even number for a clear cut 00448 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 00449 END IF 00450 * L21 Block 00451 DO I=1,N-CUT-NNB 00452 DO J=1,NNB 00453 WORK(I,J)=A(CUT+NNB+I,CUT+J) 00454 END DO 00455 END DO 00456 * L11 Block 00457 DO I=1,NNB 00458 WORK(U11+I,I)=ONE 00459 DO J=I+1,NNB 00460 WORK(U11+I,J)=ZERO 00461 END DO 00462 DO J=1,I-1 00463 WORK(U11+I,J)=A(CUT+I,CUT+J) 00464 END DO 00465 END DO 00466 * 00467 * invD*L21 00468 * 00469 I=N-CUT-NNB 00470 DO WHILE (I .GE. 1) 00471 IF (IPIV(CUT+NNB+I) > 0) THEN 00472 DO J=1,NNB 00473 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J) 00474 END DO 00475 I=I-1 00476 ELSE 00477 DO J=1,NNB 00478 U01_I_J = WORK(I,J) 00479 U01_IP1_J = WORK(I-1,J) 00480 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+ 00481 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J 00482 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+ 00483 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J 00484 END DO 00485 I=I-2 00486 END IF 00487 END DO 00488 * 00489 * invD1*L11 00490 * 00491 I=NNB 00492 DO WHILE (I .GE. 1) 00493 IF (IPIV(CUT+I) > 0) THEN 00494 DO J=1,NNB 00495 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 00496 END DO 00497 I=I-1 00498 ELSE 00499 DO J=1,NNB 00500 U11_I_J = WORK(U11+I,J) 00501 U11_IP1_J = WORK(U11+I-1,J) 00502 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 00503 $ WORK(CUT+I,INVD+1)*U11_IP1_J 00504 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+ 00505 $ WORK(CUT+I-1,INVD)*U11_IP1_J 00506 END DO 00507 I=I-2 00508 END IF 00509 END DO 00510 * 00511 * L11**T*invD1*L11->L11 00512 * 00513 CALL CTRMM('L',UPLO,'T','U',NNB, NNB, 00514 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 00515 * 00516 DO I=1,NNB 00517 DO J=1,I 00518 A(CUT+I,CUT+J)=WORK(U11+I,J) 00519 END DO 00520 END DO 00521 * 00522 IF ( (CUT+NNB) .LT. N ) THEN 00523 * 00524 * L21**T*invD2*L21->A(CUT+I,CUT+J) 00525 * 00526 CALL CGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1) 00527 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 00528 00529 * 00530 * L11 = L11**T*invD1*L11 + U01**T*invD*U01 00531 * 00532 DO I=1,NNB 00533 DO J=1,I 00534 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 00535 END DO 00536 END DO 00537 * 00538 * L01 = L22**T*invD2*L21 00539 * 00540 CALL CTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB, 00541 $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1) 00542 00543 * Update L21 00544 DO I=1,N-CUT-NNB 00545 DO J=1,NNB 00546 A(CUT+NNB+I,CUT+J)=WORK(I,J) 00547 END DO 00548 END DO 00549 ELSE 00550 * 00551 * L11 = L11**T*invD1*L11 00552 * 00553 DO I=1,NNB 00554 DO J=1,I 00555 A(CUT+I,CUT+J)=WORK(U11+I,J) 00556 END DO 00557 END DO 00558 END IF 00559 * 00560 * Next Block 00561 * 00562 CUT=CUT+NNB 00563 END DO 00564 * 00565 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T 00566 * 00567 I=N 00568 DO WHILE ( I .GE. 1 ) 00569 IF( IPIV(I) .GT. 0 ) THEN 00570 IP=IPIV(I) 00571 IF (I .LT. IP) CALL CSYSWAPR( UPLO, N, A, LDA, I ,IP ) 00572 IF (I .GT. IP) CALL CSYSWAPR( UPLO, N, A, LDA, IP ,I ) 00573 ELSE 00574 IP=-IPIV(I) 00575 IF ( I .LT. IP) CALL CSYSWAPR( UPLO, N, A, LDA, I ,IP ) 00576 IF ( I .GT. IP) CALL CSYSWAPR( UPLO, N, A, LDA, IP ,I ) 00577 I=I-1 00578 ENDIF 00579 I=I-1 00580 END DO 00581 END IF 00582 * 00583 RETURN 00584 * 00585 * End of CSYTRI2X 00586 * 00587 END 00588