![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSYTRI2X 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DSYTRI2X + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytri2x.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytri2x.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytri2x.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DSYTRI2X( 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 * DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> DSYTRI2X 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 *> DSYTRF. 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 DOUBLE PRECISION 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 DSYTRF. 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 DSYTRF. 00086 *> \endverbatim 00087 *> 00088 *> \param[out] WORK 00089 *> \verbatim 00090 *> WORK is DOUBLE PRECISION 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 doubleSYcomputational 00119 * 00120 * ===================================================================== 00121 SUBROUTINE DSYTRI2X( 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 DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* ) 00135 * .. 00136 * 00137 * ===================================================================== 00138 * 00139 * .. Parameters .. 00140 DOUBLE PRECISION ONE, ZERO 00141 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00142 * .. 00143 * .. Local Scalars .. 00144 LOGICAL UPPER 00145 INTEGER I, IINFO, IP, K, CUT, NNB 00146 INTEGER COUNT 00147 INTEGER J, U11, INVD 00148 00149 DOUBLE PRECISION AK, AKKP1, AKP1, D, T 00150 DOUBLE PRECISION U01_I_J, U01_IP1_J 00151 DOUBLE PRECISION U11_I_J, U11_IP1_J 00152 * .. 00153 * .. External Functions .. 00154 LOGICAL LSAME 00155 EXTERNAL LSAME 00156 * .. 00157 * .. External Subroutines .. 00158 EXTERNAL DSYCONV, XERBLA, DTRTRI 00159 EXTERNAL DGEMM, DTRMM, DSYSWAPR 00160 * .. 00161 * .. Intrinsic Functions .. 00162 INTRINSIC MAX 00163 * .. 00164 * .. Executable Statements .. 00165 * 00166 * Test the input parameters. 00167 * 00168 INFO = 0 00169 UPPER = LSAME( UPLO, 'U' ) 00170 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00171 INFO = -1 00172 ELSE IF( N.LT.0 ) THEN 00173 INFO = -2 00174 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00175 INFO = -4 00176 END IF 00177 * 00178 * Quick return if possible 00179 * 00180 * 00181 IF( INFO.NE.0 ) THEN 00182 CALL XERBLA( 'DSYTRI2X', -INFO ) 00183 RETURN 00184 END IF 00185 IF( N.EQ.0 ) 00186 $ RETURN 00187 * 00188 * Convert A 00189 * Workspace got Non-diag elements of D 00190 * 00191 CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 00192 * 00193 * Check that the diagonal matrix D is nonsingular. 00194 * 00195 IF( UPPER ) THEN 00196 * 00197 * Upper triangular storage: examine D from bottom to top 00198 * 00199 DO INFO = N, 1, -1 00200 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 00201 $ RETURN 00202 END DO 00203 ELSE 00204 * 00205 * Lower triangular storage: examine D from top to bottom. 00206 * 00207 DO INFO = 1, N 00208 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 00209 $ RETURN 00210 END DO 00211 END IF 00212 INFO = 0 00213 * 00214 * Splitting Workspace 00215 * U01 is a block (N,NB+1) 00216 * The first element of U01 is in WORK(1,1) 00217 * U11 is a block (NB+1,NB+1) 00218 * The first element of U11 is in WORK(N+1,1) 00219 U11 = N 00220 * INVD is a block (N,2) 00221 * The first element of INVD is in WORK(1,INVD) 00222 INVD = NB+2 00223 00224 IF( UPPER ) THEN 00225 * 00226 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T. 00227 * 00228 CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO ) 00229 * 00230 * inv(D) and inv(D)*inv(U) 00231 * 00232 K=1 00233 DO WHILE ( K .LE. N ) 00234 IF( IPIV( K ).GT.0 ) THEN 00235 * 1 x 1 diagonal NNB 00236 WORK(K,INVD) = ONE / A( K, K ) 00237 WORK(K,INVD+1) = 0 00238 K=K+1 00239 ELSE 00240 * 2 x 2 diagonal NNB 00241 T = WORK(K+1,1) 00242 AK = A( K, K ) / T 00243 AKP1 = A( K+1, K+1 ) / T 00244 AKKP1 = WORK(K+1,1) / T 00245 D = T*( AK*AKP1-ONE ) 00246 WORK(K,INVD) = AKP1 / D 00247 WORK(K+1,INVD+1) = AK / D 00248 WORK(K,INVD+1) = -AKKP1 / D 00249 WORK(K+1,INVD) = -AKKP1 / D 00250 K=K+2 00251 END IF 00252 END DO 00253 * 00254 * inv(U**T) = (inv(U))**T 00255 * 00256 * inv(U**T)*inv(D)*inv(U) 00257 * 00258 CUT=N 00259 DO WHILE (CUT .GT. 0) 00260 NNB=NB 00261 IF (CUT .LE. NNB) THEN 00262 NNB=CUT 00263 ELSE 00264 COUNT = 0 00265 * count negative elements, 00266 DO I=CUT+1-NNB,CUT 00267 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 00268 END DO 00269 * need a even number for a clear cut 00270 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 00271 END IF 00272 00273 CUT=CUT-NNB 00274 * 00275 * U01 Block 00276 * 00277 DO I=1,CUT 00278 DO J=1,NNB 00279 WORK(I,J)=A(I,CUT+J) 00280 END DO 00281 END DO 00282 * 00283 * U11 Block 00284 * 00285 DO I=1,NNB 00286 WORK(U11+I,I)=ONE 00287 DO J=1,I-1 00288 WORK(U11+I,J)=ZERO 00289 END DO 00290 DO J=I+1,NNB 00291 WORK(U11+I,J)=A(CUT+I,CUT+J) 00292 END DO 00293 END DO 00294 * 00295 * invD*U01 00296 * 00297 I=1 00298 DO WHILE (I .LE. CUT) 00299 IF (IPIV(I) > 0) THEN 00300 DO J=1,NNB 00301 WORK(I,J)=WORK(I,INVD)*WORK(I,J) 00302 END DO 00303 I=I+1 00304 ELSE 00305 DO J=1,NNB 00306 U01_I_J = WORK(I,J) 00307 U01_IP1_J = WORK(I+1,J) 00308 WORK(I,J)=WORK(I,INVD)*U01_I_J+ 00309 $ WORK(I,INVD+1)*U01_IP1_J 00310 WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+ 00311 $ WORK(I+1,INVD+1)*U01_IP1_J 00312 END DO 00313 I=I+2 00314 END IF 00315 END DO 00316 * 00317 * invD1*U11 00318 * 00319 I=1 00320 DO WHILE (I .LE. NNB) 00321 IF (IPIV(CUT+I) > 0) THEN 00322 DO J=I,NNB 00323 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 00324 END DO 00325 I=I+1 00326 ELSE 00327 DO J=I,NNB 00328 U11_I_J = WORK(U11+I,J) 00329 U11_IP1_J = WORK(U11+I+1,J) 00330 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 00331 $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J) 00332 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+ 00333 $ WORK(CUT+I+1,INVD+1)*U11_IP1_J 00334 END DO 00335 I=I+2 00336 END IF 00337 END DO 00338 * 00339 * U11**T*invD1*U11->U11 00340 * 00341 CALL DTRMM('L','U','T','U',NNB, NNB, 00342 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 00343 * 00344 DO I=1,NNB 00345 DO J=I,NNB 00346 A(CUT+I,CUT+J)=WORK(U11+I,J) 00347 END DO 00348 END DO 00349 * 00350 * U01**T*invD*U01->A(CUT+I,CUT+J) 00351 * 00352 CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA, 00353 $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 00354 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 DTRMM('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 DSYSWAPR( UPLO, N, A, LDA, I ,IP ) 00389 IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I ) 00390 ELSE 00391 IP=-IPIV(I) 00392 I=I+1 00393 IF ( (I-1) .LT. IP) 00394 $ CALL DSYSWAPR( UPLO, N, A, LDA, I-1 ,IP ) 00395 IF ( (I-1) .GT. IP) 00396 $ CALL DSYSWAPR( 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 DTRTRI( 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 .GT. 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 DTRMM('L',UPLO,'T','U',NNB, NNB, 00514 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 00515 00516 * 00517 DO I=1,NNB 00518 DO J=1,I 00519 A(CUT+I,CUT+J)=WORK(U11+I,J) 00520 END DO 00521 END DO 00522 * 00523 IF ( (CUT+NNB) .LT. N ) THEN 00524 * 00525 * L21**T*invD2*L21->A(CUT+I,CUT+J) 00526 * 00527 CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1) 00528 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 00529 00530 * 00531 * L11 = L11**T*invD1*L11 + U01**T*invD*U01 00532 * 00533 DO I=1,NNB 00534 DO J=1,I 00535 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 00536 END DO 00537 END DO 00538 * 00539 * L01 = L22**T*invD2*L21 00540 * 00541 CALL DTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB, 00542 $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1) 00543 * 00544 * Update L21 00545 * 00546 DO I=1,N-CUT-NNB 00547 DO J=1,NNB 00548 A(CUT+NNB+I,CUT+J)=WORK(I,J) 00549 END DO 00550 END DO 00551 00552 ELSE 00553 * 00554 * L11 = L11**T*invD1*L11 00555 * 00556 DO I=1,NNB 00557 DO J=1,I 00558 A(CUT+I,CUT+J)=WORK(U11+I,J) 00559 END DO 00560 END DO 00561 END IF 00562 * 00563 * Next Block 00564 * 00565 CUT=CUT+NNB 00566 END DO 00567 * 00568 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T 00569 * 00570 I=N 00571 DO WHILE ( I .GE. 1 ) 00572 IF( IPIV(I) .GT. 0 ) THEN 00573 IP=IPIV(I) 00574 IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP ) 00575 IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I ) 00576 ELSE 00577 IP=-IPIV(I) 00578 IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP ) 00579 IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP, I ) 00580 I=I-1 00581 ENDIF 00582 I=I-1 00583 END DO 00584 END IF 00585 * 00586 RETURN 00587 * 00588 * End of DSYTRI2X 00589 * 00590 END 00591