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