version 1.8, 2014/01/27 09:28:43
|
version 1.14, 2023/08/07 08:39:39
|
Line 2
|
Line 2
|
* |
* |
* =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION =========== |
* |
* |
* Online html documentation available at |
* Online html documentation available at |
* http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
*> \htmlonly |
*> \htmlonly |
*> Download ZSYTRI2X + dependencies |
*> Download ZSYTRI2X + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytri2x.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytri2x.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytri2x.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytri2x.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytri2x.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytri2x.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
* |
* |
* SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) |
* SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* CHARACTER UPLO |
* CHARACTER UPLO |
* INTEGER INFO, LDA, N, NB |
* INTEGER INFO, LDA, N, NB |
Line 28
|
Line 28
|
* INTEGER IPIV( * ) |
* INTEGER IPIV( * ) |
* COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* ) |
* COMPLEX*16 A( LDA, * ), WORK( N+NB+1,* ) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
Line 87
|
Line 87
|
*> |
*> |
*> \param[out] WORK |
*> \param[out] WORK |
*> \verbatim |
*> \verbatim |
*> WORK is COMPLEX*16 array, dimension (N+NNB+1,NNB+3) |
*> WORK is COMPLEX*16 array, dimension (N+NB+1,NB+3) |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] NB |
*> \param[in] NB |
Line 108
|
Line 108
|
* Authors: |
* Authors: |
* ======== |
* ======== |
* |
* |
*> \author Univ. of Tennessee |
*> \author Univ. of Tennessee |
*> \author Univ. of California Berkeley |
*> \author Univ. of California Berkeley |
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
|
*> \date November 2011 |
|
* |
* |
*> \ingroup complex16SYcomputational |
*> \ingroup complex16SYcomputational |
* |
* |
* ===================================================================== |
* ===================================================================== |
SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) |
SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) |
* |
* |
* -- LAPACK computational routine (version 3.4.0) -- |
* -- 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..-- |
* November 2011 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER UPLO |
CHARACTER UPLO |
Line 213
|
Line 210
|
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 229
|
Line 226
|
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 |
* 1 x 1 diagonal NNB |
* 1 x 1 diagonal NNB |
WORK(K,INVD) = 1/ A( K, K ) |
WORK(K,INVD) = ONE / A( K, K ) |
WORK(K,INVD+1) = 0 |
WORK(K,INVD+1) = 0 |
K=K+1 |
K=K+1 |
ELSE |
ELSE |
Line 246
|
Line 243
|
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) = -AKKP1 / D |
WORK(K+1,INVD) = -AKKP1 / D |
K=K+2 |
K=K+2 |
END IF |
END IF |
END DO |
END DO |
Line 263
|
Line 260
|
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 273
|
Line 270
|
|
|
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 336
|
Line 333
|
I=I+2 |
I=I+2 |
END IF |
END IF |
END DO |
END DO |
* |
* |
* U11**T*invD1*U11->U11 |
* U11**T*invD1*U11->U11 |
* |
* |
CALL ZTRMM('L','U','T','U',NNB, NNB, |
CALL ZTRMM('L','U','T','U',NNB, NNB, |
Line 346
|
Line 343
|
DO J=I,NNB |
DO J=I,NNB |
A(CUT+I,CUT+J)=WORK(U11+I,J) |
A(CUT+I,CUT+J)=WORK(U11+I,J) |
END DO |
END DO |
END DO |
END DO |
* |
* |
* U01**T*invD*U01->A(CUT+I,CUT+J) |
* U01**T*invD*U01->A(CUT+I,CUT+J) |
* |
* |
Line 380
|
Line 377
|
END DO |
END DO |
* |
* |
* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T |
* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T |
* |
* |
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 390
|
Line 387
|
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 ZSYSWAPR( UPLO, N, A, LDA, I-1 ,IP ) |
$ CALL ZSYSWAPR( UPLO, N, A, LDA, I-1 ,IP ) |
IF ( (I-1) .GT. IP) |
IF ( (I-1) .GT. IP) |
$ CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I-1 ) |
$ CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I-1 ) |
ENDIF |
ENDIF |
I=I+1 |
I=I+1 |
Line 406
|
Line 403
|
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 |
* 1 x 1 diagonal NNB |
* 1 x 1 diagonal NNB |
WORK(K,INVD) = 1/ A( K, K ) |
WORK(K,INVD) = ONE / A( K, K ) |
WORK(K,INVD+1) = 0 |
WORK(K,INVD+1) = 0 |
K=K-1 |
K=K-1 |
ELSE |
ELSE |
Line 423
|
Line 420
|
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) = -AKKP1 / D |
WORK(K-1,INVD+1) = -AKKP1 / D |
K=K-2 |
K=K-2 |
END IF |
END IF |
END DO |
END DO |
Line 440
|
Line 437
|
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 507
|
Line 504
|
I=I-2 |
I=I-2 |
END IF |
END IF |
END DO |
END DO |
* |
* |
* L11**T*invD1*L11->L11 |
* L11**T*invD1*L11->L11 |
* |
* |
CALL ZTRMM('L',UPLO,'T','U',NNB, NNB, |
CALL ZTRMM('L',UPLO,'T','U',NNB, NNB, |
Line 526
|
Line 523
|
* |
* |
CALL ZGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1) |
CALL ZGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,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**T*invD1*L11 + U01**T*invD*U01 |
* L11 = L11**T*invD1*L11 + U01**T*invD*U01 |
* |
* |
Line 564
|
Line 561
|
END DO |
END DO |
* |
* |
* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T |
* Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T |
* |
* |
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 |