![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHETRI2X 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHETRI2X + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetri2x.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetri2x.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetri2x.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHETRI2X( 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*16 A( LDA, * ), WORK( N+NB+1,* ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> ZHETRI2X computes the inverse of a COMPLEX*16 Hermitian indefinite matrix 00039 *> A using the factorization A = U*D*U**H or A = L*D*L**H computed by 00040 *> ZHETRF. 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**H; 00052 *> = 'L': Lower triangular, form is A = L*D*L**H. 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*16 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 ZHETRF. 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 ZHETRF. 00086 *> \endverbatim 00087 *> 00088 *> \param[out] WORK 00089 *> \verbatim 00090 *> WORK is COMPLEX*16 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 complex16HEcomputational 00119 * 00120 * ===================================================================== 00121 SUBROUTINE ZHETRI2X( 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*16 A( LDA, * ), WORK( N+NB+1,* ) 00135 * .. 00136 * 00137 * ===================================================================== 00138 * 00139 * .. Parameters .. 00140 REAL ONE 00141 COMPLEX*16 CONE, ZERO 00142 PARAMETER ( ONE = 1.0D+0, 00143 $ CONE = ( 1.0D+0, 0.0D+0 ), 00144 $ ZERO = ( 0.0D+0, 0.0D+0 ) ) 00145 * .. 00146 * .. Local Scalars .. 00147 LOGICAL UPPER 00148 INTEGER I, IINFO, IP, K, CUT, NNB 00149 INTEGER COUNT 00150 INTEGER J, U11, INVD 00151 00152 COMPLEX*16 AK, AKKP1, AKP1, D, T 00153 COMPLEX*16 U01_I_J, U01_IP1_J 00154 COMPLEX*16 U11_I_J, U11_IP1_J 00155 * .. 00156 * .. External Functions .. 00157 LOGICAL LSAME 00158 EXTERNAL LSAME 00159 * .. 00160 * .. External Subroutines .. 00161 EXTERNAL ZSYCONV, XERBLA, ZTRTRI 00162 EXTERNAL ZGEMM, ZTRMM, ZHESWAPR 00163 * .. 00164 * .. Intrinsic Functions .. 00165 INTRINSIC MAX 00166 * .. 00167 * .. Executable Statements .. 00168 * 00169 * Test the input parameters. 00170 * 00171 INFO = 0 00172 UPPER = LSAME( UPLO, 'U' ) 00173 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00174 INFO = -1 00175 ELSE IF( N.LT.0 ) THEN 00176 INFO = -2 00177 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00178 INFO = -4 00179 END IF 00180 * 00181 * Quick return if possible 00182 * 00183 * 00184 IF( INFO.NE.0 ) THEN 00185 CALL XERBLA( 'ZHETRI2X', -INFO ) 00186 RETURN 00187 END IF 00188 IF( N.EQ.0 ) 00189 $ RETURN 00190 * 00191 * Convert A 00192 * Workspace got Non-diag elements of D 00193 * 00194 CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) 00195 * 00196 * Check that the diagonal matrix D is nonsingular. 00197 * 00198 IF( UPPER ) THEN 00199 * 00200 * Upper triangular storage: examine D from bottom to top 00201 * 00202 DO INFO = N, 1, -1 00203 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 00204 $ RETURN 00205 END DO 00206 ELSE 00207 * 00208 * Lower triangular storage: examine D from top to bottom. 00209 * 00210 DO INFO = 1, N 00211 IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) 00212 $ RETURN 00213 END DO 00214 END IF 00215 INFO = 0 00216 * 00217 * Splitting Workspace 00218 * U01 is a block (N,NB+1) 00219 * The first element of U01 is in WORK(1,1) 00220 * U11 is a block (NB+1,NB+1) 00221 * The first element of U11 is in WORK(N+1,1) 00222 U11 = N 00223 * INVD is a block (N,2) 00224 * The first element of INVD is in WORK(1,INVD) 00225 INVD = NB+2 00226 00227 IF( UPPER ) THEN 00228 * 00229 * invA = P * inv(U**H)*inv(D)*inv(U)*P**H. 00230 * 00231 CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) 00232 * 00233 * inv(D) and inv(D)*inv(U) 00234 * 00235 K=1 00236 DO WHILE ( K .LE. N ) 00237 IF( IPIV( K ).GT.0 ) THEN 00238 * 1 x 1 diagonal NNB 00239 WORK(K,INVD) = ONE / REAL ( A( K, K ) ) 00240 WORK(K,INVD+1) = 0 00241 K=K+1 00242 ELSE 00243 * 2 x 2 diagonal NNB 00244 T = ABS ( WORK(K+1,1) ) 00245 AK = REAL ( A( K, K ) ) / T 00246 AKP1 = REAL ( A( K+1, K+1 ) ) / T 00247 AKKP1 = WORK(K+1,1) / T 00248 D = T*( AK*AKP1-ONE ) 00249 WORK(K,INVD) = AKP1 / D 00250 WORK(K+1,INVD+1) = AK / D 00251 WORK(K,INVD+1) = -AKKP1 / D 00252 WORK(K+1,INVD) = DCONJG (WORK(K,INVD+1) ) 00253 K=K+2 00254 END IF 00255 END DO 00256 * 00257 * inv(U**H) = (inv(U))**H 00258 * 00259 * inv(U**H)*inv(D)*inv(U) 00260 * 00261 CUT=N 00262 DO WHILE (CUT .GT. 0) 00263 NNB=NB 00264 IF (CUT .LE. NNB) THEN 00265 NNB=CUT 00266 ELSE 00267 COUNT = 0 00268 * count negative elements, 00269 DO I=CUT+1-NNB,CUT 00270 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 00271 END DO 00272 * need a even number for a clear cut 00273 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 00274 END IF 00275 00276 CUT=CUT-NNB 00277 * 00278 * U01 Block 00279 * 00280 DO I=1,CUT 00281 DO J=1,NNB 00282 WORK(I,J)=A(I,CUT+J) 00283 END DO 00284 END DO 00285 * 00286 * U11 Block 00287 * 00288 DO I=1,NNB 00289 WORK(U11+I,I)=CONE 00290 DO J=1,I-1 00291 WORK(U11+I,J)=ZERO 00292 END DO 00293 DO J=I+1,NNB 00294 WORK(U11+I,J)=A(CUT+I,CUT+J) 00295 END DO 00296 END DO 00297 * 00298 * invD*U01 00299 * 00300 I=1 00301 DO WHILE (I .LE. CUT) 00302 IF (IPIV(I) > 0) THEN 00303 DO J=1,NNB 00304 WORK(I,J)=WORK(I,INVD)*WORK(I,J) 00305 END DO 00306 I=I+1 00307 ELSE 00308 DO J=1,NNB 00309 U01_I_J = WORK(I,J) 00310 U01_IP1_J = WORK(I+1,J) 00311 WORK(I,J)=WORK(I,INVD)*U01_I_J+ 00312 $ WORK(I,INVD+1)*U01_IP1_J 00313 WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+ 00314 $ WORK(I+1,INVD+1)*U01_IP1_J 00315 END DO 00316 I=I+2 00317 END IF 00318 END DO 00319 * 00320 * invD1*U11 00321 * 00322 I=1 00323 DO WHILE (I .LE. NNB) 00324 IF (IPIV(CUT+I) > 0) THEN 00325 DO J=I,NNB 00326 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 00327 END DO 00328 I=I+1 00329 ELSE 00330 DO J=I,NNB 00331 U11_I_J = WORK(U11+I,J) 00332 U11_IP1_J = WORK(U11+I+1,J) 00333 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 00334 $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J) 00335 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+ 00336 $ WORK(CUT+I+1,INVD+1)*U11_IP1_J 00337 END DO 00338 I=I+2 00339 END IF 00340 END DO 00341 * 00342 * U11**H*invD1*U11->U11 00343 * 00344 CALL ZTRMM('L','U','C','U',NNB, NNB, 00345 $ CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 00346 * 00347 DO I=1,NNB 00348 DO J=I,NNB 00349 A(CUT+I,CUT+J)=WORK(U11+I,J) 00350 END DO 00351 END DO 00352 * 00353 * U01**H*invD*U01->A(CUT+I,CUT+J) 00354 * 00355 CALL ZGEMM('C','N',NNB,NNB,CUT,CONE,A(1,CUT+1),LDA, 00356 $ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 00357 * 00358 * U11 = U11**H*invD1*U11 + U01**H*invD*U01 00359 * 00360 DO I=1,NNB 00361 DO J=I,NNB 00362 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 00363 END DO 00364 END DO 00365 * 00366 * U01 = U00**H*invD0*U01 00367 * 00368 CALL ZTRMM('L',UPLO,'C','U',CUT, NNB, 00369 $ CONE,A,LDA,WORK,N+NB+1) 00370 00371 * 00372 * Update U01 00373 * 00374 DO I=1,CUT 00375 DO J=1,NNB 00376 A(I,CUT+J)=WORK(I,J) 00377 END DO 00378 END DO 00379 * 00380 * Next Block 00381 * 00382 END DO 00383 * 00384 * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H 00385 * 00386 I=1 00387 DO WHILE ( I .LE. N ) 00388 IF( IPIV(I) .GT. 0 ) THEN 00389 IP=IPIV(I) 00390 IF (I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP ) 00391 IF (I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I ) 00392 ELSE 00393 IP=-IPIV(I) 00394 I=I+1 00395 IF ( (I-1) .LT. IP) 00396 $ CALL ZHESWAPR( UPLO, N, A, LDA, I-1 ,IP ) 00397 IF ( (I-1) .GT. IP) 00398 $ CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I-1 ) 00399 ENDIF 00400 I=I+1 00401 END DO 00402 ELSE 00403 * 00404 * LOWER... 00405 * 00406 * invA = P * inv(U**H)*inv(D)*inv(U)*P**H. 00407 * 00408 CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) 00409 * 00410 * inv(D) and inv(D)*inv(U) 00411 * 00412 K=N 00413 DO WHILE ( K .GE. 1 ) 00414 IF( IPIV( K ).GT.0 ) THEN 00415 * 1 x 1 diagonal NNB 00416 WORK(K,INVD) = ONE / REAL ( A( K, K ) ) 00417 WORK(K,INVD+1) = 0 00418 K=K-1 00419 ELSE 00420 * 2 x 2 diagonal NNB 00421 T = ABS ( WORK(K-1,1) ) 00422 AK = REAL ( A( K-1, K-1 ) ) / T 00423 AKP1 = REAL ( A( K, K ) ) / T 00424 AKKP1 = WORK(K-1,1) / T 00425 D = T*( AK*AKP1-ONE ) 00426 WORK(K-1,INVD) = AKP1 / D 00427 WORK(K,INVD) = AK / D 00428 WORK(K,INVD+1) = -AKKP1 / D 00429 WORK(K-1,INVD+1) = DCONJG (WORK(K,INVD+1) ) 00430 K=K-2 00431 END IF 00432 END DO 00433 * 00434 * inv(U**H) = (inv(U))**H 00435 * 00436 * inv(U**H)*inv(D)*inv(U) 00437 * 00438 CUT=0 00439 DO WHILE (CUT .LT. N) 00440 NNB=NB 00441 IF (CUT + NNB .GE. N) THEN 00442 NNB=N-CUT 00443 ELSE 00444 COUNT = 0 00445 * count negative elements, 00446 DO I=CUT+1,CUT+NNB 00447 IF (IPIV(I) .LT. 0) COUNT=COUNT+1 00448 END DO 00449 * need a even number for a clear cut 00450 IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 00451 END IF 00452 * L21 Block 00453 DO I=1,N-CUT-NNB 00454 DO J=1,NNB 00455 WORK(I,J)=A(CUT+NNB+I,CUT+J) 00456 END DO 00457 END DO 00458 * L11 Block 00459 DO I=1,NNB 00460 WORK(U11+I,I)=CONE 00461 DO J=I+1,NNB 00462 WORK(U11+I,J)=ZERO 00463 END DO 00464 DO J=1,I-1 00465 WORK(U11+I,J)=A(CUT+I,CUT+J) 00466 END DO 00467 END DO 00468 * 00469 * invD*L21 00470 * 00471 I=N-CUT-NNB 00472 DO WHILE (I .GE. 1) 00473 IF (IPIV(CUT+NNB+I) > 0) THEN 00474 DO J=1,NNB 00475 WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J) 00476 END DO 00477 I=I-1 00478 ELSE 00479 DO J=1,NNB 00480 U01_I_J = WORK(I,J) 00481 U01_IP1_J = WORK(I-1,J) 00482 WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+ 00483 $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J 00484 WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+ 00485 $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J 00486 END DO 00487 I=I-2 00488 END IF 00489 END DO 00490 * 00491 * invD1*L11 00492 * 00493 I=NNB 00494 DO WHILE (I .GE. 1) 00495 IF (IPIV(CUT+I) > 0) THEN 00496 DO J=1,NNB 00497 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) 00498 END DO 00499 I=I-1 00500 ELSE 00501 DO J=1,NNB 00502 U11_I_J = WORK(U11+I,J) 00503 U11_IP1_J = WORK(U11+I-1,J) 00504 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + 00505 $ WORK(CUT+I,INVD+1)*U11_IP1_J 00506 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+ 00507 $ WORK(CUT+I-1,INVD)*U11_IP1_J 00508 END DO 00509 I=I-2 00510 END IF 00511 END DO 00512 * 00513 * L11**H*invD1*L11->L11 00514 * 00515 CALL ZTRMM('L',UPLO,'C','U',NNB, NNB, 00516 $ CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) 00517 * 00518 DO I=1,NNB 00519 DO J=1,I 00520 A(CUT+I,CUT+J)=WORK(U11+I,J) 00521 END DO 00522 END DO 00523 * 00524 IF ( (CUT+NNB) .LT. N ) THEN 00525 * 00526 * L21**H*invD2*L21->A(CUT+I,CUT+J) 00527 * 00528 CALL ZGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1) 00529 $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) 00530 00531 * 00532 * L11 = L11**H*invD1*L11 + U01**H*invD*U01 00533 * 00534 DO I=1,NNB 00535 DO J=1,I 00536 A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) 00537 END DO 00538 END DO 00539 * 00540 * L01 = L22**H*invD2*L21 00541 * 00542 CALL ZTRMM('L',UPLO,'C','U', N-NNB-CUT, NNB, 00543 $ CONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1) 00544 00545 * Update L21 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 ELSE 00552 * 00553 * L11 = L11**H*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**H: P * inv(U**H)*inv(D)*inv(U) *P**H 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 ZHESWAPR( UPLO, N, A, LDA, I ,IP ) 00574 IF (I .GT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I ) 00575 ELSE 00576 IP=-IPIV(I) 00577 IF ( I .LT. IP) CALL ZHESWAPR( UPLO, N, A, LDA, I ,IP ) 00578 IF ( I .GT. IP) CALL ZHESWAPR( 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 ZHETRI2X 00588 * 00589 END 00590