![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSYTRI2X 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SSYTRI2X + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytri2x.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytri2x.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytri2x.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SSYTRI2X( 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 * REAL A( LDA, * ), WORK( N+NB+1,* ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> SSYTRI2X 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 *> SSYTRF. 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 REAL 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 SSYTRF. 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 SSYTRF. 00086 *> \endverbatim 00087 *> 00088 *> \param[out] WORK 00089 *> \verbatim 00090 *> WORK is REAL 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 realSYcomputational 00119 * 00120 * ===================================================================== 00121 SUBROUTINE SSYTRI2X( 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 REAL A( LDA, * ), WORK( N+NB+1,* ) 00135 * .. 00136 * 00137 * ===================================================================== 00138 * 00139 * .. Parameters .. 00140 REAL ONE, ZERO 00141 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+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 REAL AK, AKKP1, AKP1, D, T 00150 REAL U01_I_J, U01_IP1_J 00151 REAL U11_I_J, U11_IP1_J 00152 * .. 00153 * .. External Functions .. 00154 LOGICAL LSAME 00155 EXTERNAL LSAME 00156 * .. 00157 * .. External Subroutines .. 00158 EXTERNAL SSYCONV, XERBLA, STRTRI 00159 EXTERNAL SGEMM, STRMM, SSYSWAPR 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( 'SSYTRI2X', -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 SSYCONV( 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 STRTRI( 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 STRMM('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 SGEMM('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 * U11 = U11**T*invD1*U11 + U01**T*invD*U01 00356 * 00357 DO I=1,NNB 00358 DO J=I,NNB 00359 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 00360 END DO 00361 END DO 00362 * 00363 * U01 = U00**T*invD0*U01 00364 * 00365 CALL STRMM('L',UPLO,'T','U',CUT, NNB, 00366 $ ONE,A,LDA,WORK,N+NB+1) 00367 00368 * 00369 * Update U01 00370 * 00371 DO I=1,CUT 00372 DO J=1,NNB 00373 A(I,CUT+J)=WORK(I,J) 00374 END DO 00375 END DO 00376 * 00377 * Next Block 00378 * 00379 END DO 00380 * 00381 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T 00382 * 00383 I=1 00384 DO WHILE ( I .LE. N ) 00385 IF( IPIV(I) .GT. 0 ) THEN 00386 IP=IPIV(I) 00387 IF (I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP ) 00388 IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I ) 00389 ELSE 00390 IP=-IPIV(I) 00391 I=I+1 00392 IF ( (I-1) .LT. IP) 00393 $ CALL SSYSWAPR( UPLO, N, A, LDA, I-1 ,IP ) 00394 IF ( (I-1) .GT. IP) 00395 $ CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I-1 ) 00396 ENDIF 00397 I=I+1 00398 END DO 00399 ELSE 00400 * 00401 * LOWER... 00402 * 00403 * invA = P * inv(U**T)*inv(D)*inv(U)*P**T. 00404 * 00405 CALL STRTRI( UPLO, 'U', N, A, LDA, INFO ) 00406 * 00407 * inv(D) and inv(D)*inv(U) 00408 * 00409 K=N 00410 DO WHILE ( K .GE. 1 ) 00411 IF( IPIV( K ).GT.0 ) THEN 00412 * 1 x 1 diagonal NNB 00413 WORK(K,INVD) = ONE / A( K, K ) 00414 WORK(K,INVD+1) = 0 00415 K=K-1 00416 ELSE 00417 * 2 x 2 diagonal NNB 00418 T = WORK(K-1,1) 00419 AK = A( K-1, K-1 ) / T 00420 AKP1 = A( K, K ) / T 00421 AKKP1 = WORK(K-1,1) / T 00422 D = T*( AK*AKP1-ONE ) 00423 WORK(K-1,INVD) = AKP1 / D 00424 WORK(K,INVD) = AK / D 00425 WORK(K,INVD+1) = -AKKP1 / D 00426 WORK(K-1,INVD+1) = -AKKP1 / D 00427 K=K-2 00428 END IF 00429 END DO 00430 * 00431 * inv(U**T) = (inv(U))**T 00432 * 00433 * inv(U**T)*inv(D)*inv(U) 00434 * 00435 CUT=0 00436 DO WHILE (CUT .LT. N) 00437 NNB=NB 00438 IF (CUT + NNB .GT. N) THEN 00439 NNB=N-CUT 00440 ELSE 00441 COUNT = 0 00442 * count negative elements, 00443 DO I=CUT+1,CUT+NNB 00444 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 00445 END DO 00446 * need a even number for a clear cut 00447 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 00448 END IF 00449 * L21 Block 00450 DO I=1,N-CUT-NNB 00451 DO J=1,NNB 00452 WORK(I,J)=A(CUT+NNB+I,CUT+J) 00453 END DO 00454 END DO 00455 * L11 Block 00456 DO I=1,NNB 00457 WORK(U11+I,I)=ONE 00458 DO J=I+1,NNB 00459 WORK(U11+I,J)=ZERO 00460 END DO 00461 DO J=1,I-1 00462 WORK(U11+I,J)=A(CUT+I,CUT+J) 00463 END DO 00464 END DO 00465 * 00466 * invD*L21 00467 * 00468 I=N-CUT-NNB 00469 DO WHILE (I .GE. 1) 00470 IF (IPIV(CUT+NNB+I) > 0) THEN 00471 DO J=1,NNB 00472 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J) 00473 END DO 00474 I=I-1 00475 ELSE 00476 DO J=1,NNB 00477 U01_I_J = WORK(I,J) 00478 U01_IP1_J = WORK(I-1,J) 00479 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+ 00480 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J 00481 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+ 00482 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J 00483 END DO 00484 I=I-2 00485 END IF 00486 END DO 00487 * 00488 * invD1*L11 00489 * 00490 I=NNB 00491 DO WHILE (I .GE. 1) 00492 IF (IPIV(CUT+I) > 0) THEN 00493 DO J=1,NNB 00494 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 00495 END DO 00496 I=I-1 00497 ELSE 00498 DO J=1,NNB 00499 U11_I_J = WORK(U11+I,J) 00500 U11_IP1_J = WORK(U11+I-1,J) 00501 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 00502 $ WORK(CUT+I,INVD+1)*U11_IP1_J 00503 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+ 00504 $ WORK(CUT+I-1,INVD)*U11_IP1_J 00505 END DO 00506 I=I-2 00507 END IF 00508 END DO 00509 * 00510 * L11**T*invD1*L11->L11 00511 * 00512 CALL STRMM('L',UPLO,'T','U',NNB, NNB, 00513 $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 00514 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 SGEMM('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 STRMM('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 * 00545 DO I=1,N-CUT-NNB 00546 DO J=1,NNB 00547 A(CUT+NNB+I,CUT+J)=WORK(I,J) 00548 END DO 00549 END DO 00550 00551 ELSE 00552 * 00553 * L11 = L11**T*invD1*L11 00554 * 00555 DO I=1,NNB 00556 DO J=1,I 00557 A(CUT+I,CUT+J)=WORK(U11+I,J) 00558 END DO 00559 END DO 00560 END IF 00561 * 00562 * Next Block 00563 * 00564 CUT=CUT+NNB 00565 END DO 00566 * 00567 * Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T 00568 * 00569 I=N 00570 DO WHILE ( I .GE. 1 ) 00571 IF( IPIV(I) .GT. 0 ) THEN 00572 IP=IPIV(I) 00573 IF (I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP ) 00574 IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I ) 00575 ELSE 00576 IP=-IPIV(I) 00577 IF ( I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP ) 00578 IF ( I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I ) 00579 I=I-1 00580 ENDIF 00581 I=I-1 00582 END DO 00583 END IF 00584 * 00585 RETURN 00586 * 00587 * End of SSYTRI2X 00588 * 00589 END 00590