|
version 1.1, 2011/07/24 10:31:39
|
version 1.13, 2023/08/07 08:39:25
|
|
Line 1
|
Line 1
|
| |
*> \brief \b ZHETRI2X |
| |
* |
| |
* =========== DOCUMENTATION =========== |
| |
* |
| |
* Online html documentation available at |
| |
* http://www.netlib.org/lapack/explore-html/ |
| |
* |
| |
*> \htmlonly |
| |
*> Download ZHETRI2X + dependencies |
| |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetri2x.f"> |
| |
*> [TGZ]</a> |
| |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetri2x.f"> |
| |
*> [ZIP]</a> |
| |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetri2x.f"> |
| |
*> [TXT]</a> |
| |
*> \endhtmlonly |
| |
* |
| |
* Definition: |
| |
* =========== |
| |
* |
| |
* SUBROUTINE ZHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) |
| |
* |
| |
* .. Scalar Arguments .. |
| |
* CHARACTER UPLO |
| |
* INTEGER INFO, LDA, N, NB |
| |
* .. |
| |
* .. Array Arguments .. |
| |
* INTEGER IPIV( * ) |
| |
* COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* ) |
| |
* .. |
| |
* |
| |
* |
| |
*> \par Purpose: |
| |
* ============= |
| |
*> |
| |
*> \verbatim |
| |
*> |
| |
*> ZHETRI2X computes the inverse of a COMPLEX*16 Hermitian indefinite matrix |
| |
*> A using the factorization A = U*D*U**H or A = L*D*L**H computed by |
| |
*> ZHETRF. |
| |
*> \endverbatim |
| |
* |
| |
* Arguments: |
| |
* ========== |
| |
* |
| |
*> \param[in] UPLO |
| |
*> \verbatim |
| |
*> UPLO is CHARACTER*1 |
| |
*> Specifies whether the details of the factorization are stored |
| |
*> as an upper or lower triangular matrix. |
| |
*> = 'U': Upper triangular, form is A = U*D*U**H; |
| |
*> = 'L': Lower triangular, form is A = L*D*L**H. |
| |
*> \endverbatim |
| |
*> |
| |
*> \param[in] N |
| |
*> \verbatim |
| |
*> N is INTEGER |
| |
*> The order of the matrix A. N >= 0. |
| |
*> \endverbatim |
| |
*> |
| |
*> \param[in,out] A |
| |
*> \verbatim |
| |
*> A is COMPLEX*16 array, dimension (LDA,N) |
| |
*> On entry, the NNB diagonal matrix D and the multipliers |
| |
*> used to obtain the factor U or L as computed by ZHETRF. |
| |
*> |
| |
*> On exit, if INFO = 0, the (symmetric) inverse of the original |
| |
*> matrix. If UPLO = 'U', the upper triangular part of the |
| |
*> inverse is formed and the part of A below the diagonal is not |
| |
*> referenced; if UPLO = 'L' the lower triangular part of the |
| |
*> inverse is formed and the part of A above the diagonal is |
| |
*> not referenced. |
| |
*> \endverbatim |
| |
*> |
| |
*> \param[in] LDA |
| |
*> \verbatim |
| |
*> LDA is INTEGER |
| |
*> The leading dimension of the array A. LDA >= max(1,N). |
| |
*> \endverbatim |
| |
*> |
| |
*> \param[in] IPIV |
| |
*> \verbatim |
| |
*> IPIV is INTEGER array, dimension (N) |
| |
*> Details of the interchanges and the NNB structure of D |
| |
*> as determined by ZHETRF. |
| |
*> \endverbatim |
| |
*> |
| |
*> \param[out] WORK |
| |
*> \verbatim |
| |
*> WORK is COMPLEX*16 array, dimension (N+NB+1,NB+3) |
| |
*> \endverbatim |
| |
*> |
| |
*> \param[in] NB |
| |
*> \verbatim |
| |
*> NB is INTEGER |
| |
*> Block size |
| |
*> \endverbatim |
| |
*> |
| |
*> \param[out] INFO |
| |
*> \verbatim |
| |
*> INFO is INTEGER |
| |
*> = 0: successful exit |
| |
*> < 0: if INFO = -i, the i-th argument had an illegal value |
| |
*> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its |
| |
*> inverse could not be computed. |
| |
*> \endverbatim |
| |
* |
| |
* Authors: |
| |
* ======== |
| |
* |
| |
*> \author Univ. of Tennessee |
| |
*> \author Univ. of California Berkeley |
| |
*> \author Univ. of Colorado Denver |
| |
*> \author NAG Ltd. |
| |
* |
| |
*> \ingroup complex16HEcomputational |
| |
* |
| |
* ===================================================================== |
| SUBROUTINE ZHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) |
SUBROUTINE ZHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) |
| * |
* |
| * -- LAPACK routine (version 3.3.1) -- |
* -- LAPACK computational routine -- |
| * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
| * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
| * -- April 2011 -- |
|
| * |
|
| * -- Written by Julie Langou of the Univ. of TN -- |
|
| * |
* |
| * .. Scalar Arguments .. |
* .. Scalar Arguments .. |
| CHARACTER UPLO |
CHARACTER UPLO |
|
Line 16
|
Line 131
|
| COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* ) |
COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* ) |
| * .. |
* .. |
| * |
* |
| * Purpose |
|
| * ======= |
|
| * |
|
| * ZHETRI2X computes the inverse of a COMPLEX*16 Hermitian indefinite matrix |
|
| * A using the factorization A = U*D*U**H or A = L*D*L**H computed by |
|
| * ZHETRF. |
|
| * |
|
| * Arguments |
|
| * ========= |
|
| * |
|
| * UPLO (input) CHARACTER*1 |
|
| * Specifies whether the details of the factorization are stored |
|
| * as an upper or lower triangular matrix. |
|
| * = 'U': Upper triangular, form is A = U*D*U**H; |
|
| * = 'L': Lower triangular, form is A = L*D*L**H. |
|
| * |
|
| * N (input) INTEGER |
|
| * The order of the matrix A. N >= 0. |
|
| * |
|
| * A (input/output) COMPLEX*16 array, dimension (LDA,N) |
|
| * On entry, the NNB diagonal matrix D and the multipliers |
|
| * used to obtain the factor U or L as computed by ZHETRF. |
|
| * |
|
| * On exit, if INFO = 0, the (symmetric) inverse of the original |
|
| * matrix. If UPLO = 'U', the upper triangular part of the |
|
| * inverse is formed and the part of A below the diagonal is not |
|
| * referenced; if UPLO = 'L' the lower triangular part of the |
|
| * inverse is formed and the part of A above the diagonal is |
|
| * not referenced. |
|
| * |
|
| * LDA (input) INTEGER |
|
| * The leading dimension of the array A. LDA >= max(1,N). |
|
| * |
|
| * IPIV (input) INTEGER array, dimension (N) |
|
| * Details of the interchanges and the NNB structure of D |
|
| * as determined by ZHETRF. |
|
| * |
|
| * WORK (workspace) COMPLEX*16 array, dimension (N+NNB+1,NNB+3) |
|
| * |
|
| * NB (input) INTEGER |
|
| * Block size |
|
| * |
|
| * INFO (output) INTEGER |
|
| * = 0: successful exit |
|
| * < 0: if INFO = -i, the i-th argument had an illegal value |
|
| * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its |
|
| * inverse could not be computed. |
|
| * |
|
| * ===================================================================== |
* ===================================================================== |
| * |
* |
| * .. Parameters .. |
* .. Parameters .. |
| REAL ONE |
DOUBLE PRECISION ONE |
| COMPLEX*16 CONE, ZERO |
COMPLEX*16 CONE, ZERO |
| PARAMETER ( ONE = 1.0D+0, |
PARAMETER ( ONE = 1.0D+0, |
| $ CONE = ( 1.0D+0, 0.0D+0 ), |
$ CONE = ( 1.0D+0, 0.0D+0 ), |
|
Line 145
|
Line 212
|
| INFO = 0 |
INFO = 0 |
| * |
* |
| * Splitting Workspace |
* Splitting Workspace |
| * U01 is a block (N,NB+1) |
* U01 is a block (N,NB+1) |
| * The first element of U01 is in WORK(1,1) |
* The first element of U01 is in WORK(1,1) |
| * U11 is a block (NB+1,NB+1) |
* U11 is a block (NB+1,NB+1) |
| * The first element of U11 is in WORK(N+1,1) |
* The first element of U11 is in WORK(N+1,1) |
|
Line 161
|
Line 228
|
| CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) |
CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) |
| * |
* |
| * inv(D) and inv(D)*inv(U) |
* inv(D) and inv(D)*inv(U) |
| * |
* |
| K=1 |
K=1 |
| DO WHILE ( K .LE. N ) |
DO WHILE ( K .LE. N ) |
| IF( IPIV( K ).GT.0 ) THEN |
IF( IPIV( K ).GT.0 ) THEN |
|
Line 172
|
Line 239
|
| ELSE |
ELSE |
| * 2 x 2 diagonal NNB |
* 2 x 2 diagonal NNB |
| T = ABS ( WORK(K+1,1) ) |
T = ABS ( WORK(K+1,1) ) |
| AK = REAL ( A( K, K ) ) / T |
AK = DBLE ( A( K, K ) ) / T |
| AKP1 = REAL ( A( K+1, K+1 ) ) / T |
AKP1 = DBLE ( A( K+1, K+1 ) ) / T |
| AKKP1 = WORK(K+1,1) / T |
AKKP1 = WORK(K+1,1) / T |
| D = T*( AK*AKP1-ONE ) |
D = T*( AK*AKP1-ONE ) |
| WORK(K,INVD) = AKP1 / D |
WORK(K,INVD) = AKP1 / D |
| WORK(K+1,INVD+1) = AK / D |
WORK(K+1,INVD+1) = AK / D |
| WORK(K,INVD+1) = -AKKP1 / D |
WORK(K,INVD+1) = -AKKP1 / D |
| WORK(K+1,INVD) = DCONJG (WORK(K,INVD+1) ) |
WORK(K+1,INVD) = DCONJG (WORK(K,INVD+1) ) |
| K=K+2 |
K=K+2 |
| END IF |
END IF |
|
Line 195
|
Line 262
|
| NNB=CUT |
NNB=CUT |
| ELSE |
ELSE |
| COUNT = 0 |
COUNT = 0 |
| * count negative elements, |
* count negative elements, |
| DO I=CUT+1-NNB,CUT |
DO I=CUT+1-NNB,CUT |
| IF (IPIV(I) .LT. 0) COUNT=COUNT+1 |
IF (IPIV(I) .LT. 0) COUNT=COUNT+1 |
| END DO |
END DO |
|
Line 205
|
Line 272
|
| |
|
| CUT=CUT-NNB |
CUT=CUT-NNB |
| * |
* |
| * U01 Block |
* U01 Block |
| * |
* |
| DO I=1,CUT |
DO I=1,CUT |
| DO J=1,NNB |
DO J=1,NNB |
|
Line 268
|
Line 335
|
| I=I+2 |
I=I+2 |
| END IF |
END IF |
| END DO |
END DO |
| * |
* |
| * U11**H*invD1*U11->U11 |
* U11**H*invD1*U11->U11 |
| * |
* |
| CALL ZTRMM('L','U','C','U',NNB, NNB, |
CALL ZTRMM('L','U','C','U',NNB, NNB, |
|
Line 312
|
Line 379
|
| END DO |
END DO |
| * |
* |
| * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H |
* Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H |
| * |
* |
| I=1 |
I=1 |
| DO WHILE ( I .LE. N ) |
DO WHILE ( I .LE. N ) |
| IF( IPIV(I) .GT. 0 ) THEN |
IF( IPIV(I) .GT. 0 ) THEN |
|
Line 322
|
Line 389
|
| ELSE |
ELSE |
| IP=-IPIV(I) |
IP=-IPIV(I) |
| I=I+1 |
I=I+1 |
| IF ( (I-1) .LT. IP) |
IF ( (I-1) .LT. IP) |
| $ CALL ZHESWAPR( UPLO, N, A, LDA, I-1 ,IP ) |
$ CALL ZHESWAPR( UPLO, N, A, LDA, I-1 ,IP ) |
| IF ( (I-1) .GT. IP) |
IF ( (I-1) .GT. IP) |
| $ CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I-1 ) |
$ CALL ZHESWAPR( UPLO, N, A, LDA, IP ,I-1 ) |
| ENDIF |
ENDIF |
| I=I+1 |
I=I+1 |
|
Line 338
|
Line 405
|
| CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) |
CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) |
| * |
* |
| * inv(D) and inv(D)*inv(U) |
* inv(D) and inv(D)*inv(U) |
| * |
* |
| K=N |
K=N |
| DO WHILE ( K .GE. 1 ) |
DO WHILE ( K .GE. 1 ) |
| IF( IPIV( K ).GT.0 ) THEN |
IF( IPIV( K ).GT.0 ) THEN |
|
Line 349
|
Line 416
|
| ELSE |
ELSE |
| * 2 x 2 diagonal NNB |
* 2 x 2 diagonal NNB |
| T = ABS ( WORK(K-1,1) ) |
T = ABS ( WORK(K-1,1) ) |
| AK = REAL ( A( K-1, K-1 ) ) / T |
AK = DBLE ( A( K-1, K-1 ) ) / T |
| AKP1 = REAL ( A( K, K ) ) / T |
AKP1 = DBLE ( A( K, K ) ) / T |
| AKKP1 = WORK(K-1,1) / T |
AKKP1 = WORK(K-1,1) / T |
| D = T*( AK*AKP1-ONE ) |
D = T*( AK*AKP1-ONE ) |
| WORK(K-1,INVD) = AKP1 / D |
WORK(K-1,INVD) = AKP1 / D |
| WORK(K,INVD) = AK / D |
WORK(K,INVD) = AK / D |
| WORK(K,INVD+1) = -AKKP1 / D |
WORK(K,INVD+1) = -AKKP1 / D |
| WORK(K-1,INVD+1) = DCONJG (WORK(K,INVD+1) ) |
WORK(K-1,INVD+1) = DCONJG (WORK(K,INVD+1) ) |
| K=K-2 |
K=K-2 |
| END IF |
END IF |
|
Line 372
|
Line 439
|
| NNB=N-CUT |
NNB=N-CUT |
| ELSE |
ELSE |
| COUNT = 0 |
COUNT = 0 |
| * count negative elements, |
* count negative elements, |
| DO I=CUT+1,CUT+NNB |
DO I=CUT+1,CUT+NNB |
| IF (IPIV(I) .LT. 0) COUNT=COUNT+1 |
IF (IPIV(I) .LT. 0) COUNT=COUNT+1 |
| END DO |
END DO |
|
Line 439
|
Line 506
|
| I=I-2 |
I=I-2 |
| END IF |
END IF |
| END DO |
END DO |
| * |
* |
| * L11**H*invD1*L11->L11 |
* L11**H*invD1*L11->L11 |
| * |
* |
| CALL ZTRMM('L',UPLO,'C','U',NNB, NNB, |
CALL ZTRMM('L',UPLO,'C','U',NNB, NNB, |
|
Line 457
|
Line 524
|
| * |
* |
| CALL ZGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1) |
CALL ZGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1) |
| $ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) |
$ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) |
| |
|
| * |
* |
| * L11 = L11**H*invD1*L11 + U01**H*invD*U01 |
* L11 = L11**H*invD1*L11 + U01**H*invD*U01 |
| * |
* |
|
Line 495
|
Line 562
|
| END DO |
END DO |
| * |
* |
| * Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H |
* Apply PERMUTATIONS P and P**H: P * inv(U**H)*inv(D)*inv(U) *P**H |
| * |
* |
| I=N |
I=N |
| DO WHILE ( I .GE. 1 ) |
DO WHILE ( I .GE. 1 ) |
| IF( IPIV(I) .GT. 0 ) THEN |
IF( IPIV(I) .GT. 0 ) THEN |
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetri2x.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack