version 1.2, 2010/12/21 13:53:56
|
version 1.3, 2011/07/22 07:38:20
|
Line 1
|
Line 1
|
SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) |
SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) |
* |
* |
* -- LAPACK routine (version 3.3.0) -- |
* -- LAPACK routine (version 3.3.1) -- |
* -- 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 2010 |
* -- April 2011 -- |
* |
* |
* -- Written by Julie Langou of the Univ. of TN -- |
* -- Written by Julie Langou of the Univ. of TN -- |
* |
* |
Line 154
|
Line 154
|
|
|
IF( UPPER ) THEN |
IF( UPPER ) THEN |
* |
* |
* invA = P * inv(U')*inv(D)*inv(U)*P'. |
* invA = P * inv(U**T)*inv(D)*inv(U)*P**T. |
* |
* |
CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) |
CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) |
* |
* |
Line 182
|
Line 182
|
END IF |
END IF |
END DO |
END DO |
* |
* |
* inv(U') = (inv(U))' |
* inv(U**T) = (inv(U))**T |
* |
* |
* inv(U')*inv(D)*inv(U) |
* inv(U**T)*inv(D)*inv(U) |
* |
* |
CUT=N |
CUT=N |
DO WHILE (CUT .GT. 0) |
DO WHILE (CUT .GT. 0) |
Line 267
|
Line 267
|
END IF |
END IF |
END DO |
END DO |
* |
* |
* U11T*invD1*U11->U11 |
* U11**T*invD1*U11->U11 |
* |
* |
CALL ZTRMM('L','U','T','U',NNB, NNB, |
CALL ZTRMM('L','U','T','U',NNB, NNB, |
$ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) |
$ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) |
* |
* |
* U01'invD*U01->A(CUT+I,CUT+J) |
DO I=1,NNB |
|
DO J=I,NNB |
|
A(CUT+I,CUT+J)=WORK(U11+I,J) |
|
END DO |
|
END DO |
|
* |
|
* U01**T*invD*U01->A(CUT+I,CUT+J) |
* |
* |
CALL ZGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA, |
CALL ZGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA, |
$ WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA) |
$ WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) |
* |
* |
* U11 = U11T*invD1*U11 + U01'invD*U01 |
* U11 = U11**T*invD1*U11 + U01**T*invD*U01 |
* |
* |
DO I=1,NNB |
DO I=1,NNB |
DO J=I,NNB |
DO J=I,NNB |
Line 285
|
Line 291
|
END DO |
END DO |
END DO |
END DO |
* |
* |
* U01 = U00T*invD0*U01 |
* U01 = U00**T*invD0*U01 |
* |
* |
CALL ZTRMM('L',UPLO,'T','U',CUT, NNB, |
CALL ZTRMM('L',UPLO,'T','U',CUT, NNB, |
$ ONE,A,LDA,WORK,N+NB+1) |
$ ONE,A,LDA,WORK,N+NB+1) |
Line 303
|
Line 309
|
* |
* |
END DO |
END DO |
* |
* |
* Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P' |
* 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 |
IP=IPIV(I) |
IP=IPIV(I) |
IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP ) |
IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP ) |
IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I ) |
IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I ) |
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, 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, IP ,I-1 ) |
$ CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I-1 ) |
ENDIF |
ENDIF |
I=I+1 |
I=I+1 |
END DO |
END DO |
Line 325
|
Line 331
|
* |
* |
* LOWER... |
* LOWER... |
* |
* |
* invA = P * inv(U')*inv(D)*inv(U)*P'. |
* invA = P * inv(U**T)*inv(D)*inv(U)*P**T. |
* |
* |
CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) |
CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO ) |
* |
* |
Line 353
|
Line 359
|
END IF |
END IF |
END DO |
END DO |
* |
* |
* inv(U') = (inv(U))' |
* inv(U**T) = (inv(U))**T |
* |
* |
* inv(U')*inv(D)*inv(U) |
* inv(U**T)*inv(D)*inv(U) |
* |
* |
CUT=0 |
CUT=0 |
DO WHILE (CUT .LT. N) |
DO WHILE (CUT .LT. N) |
Line 432
|
Line 438
|
END IF |
END IF |
END DO |
END DO |
* |
* |
* L11T*invD1*L11->L11 |
* L11**T*invD1*L11->L11 |
* |
* |
CALL ZTRMM('L',UPLO,'T','U',NNB, NNB, |
CALL ZTRMM('L',UPLO,'T','U',NNB, NNB, |
$ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) |
$ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) |
|
* |
|
DO I=1,NNB |
|
DO J=1,I |
|
A(CUT+I,CUT+J)=WORK(U11+I,J) |
|
END DO |
|
END DO |
|
* |
|
|
IF ( (CUT+NNB) .LT. N ) THEN |
IF ( (CUT+NNB) .LT. N ) THEN |
* |
* |
* L21T*invD2*L21->A(CUT+I,CUT+J) |
* L21**T*invD2*L21->A(CUT+I,CUT+J) |
* |
* |
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, A(CUT+1,CUT+1), LDA) |
$ ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1) |
|
|
* |
* |
* L11 = L11T*invD1*L11 + U01'invD*U01 |
* L11 = L11**T*invD1*L11 + U01**T*invD*U01 |
* |
* |
DO I=1,NNB |
DO I=1,NNB |
DO J=1,I |
DO J=1,I |
Line 453
|
Line 466
|
END DO |
END DO |
END DO |
END DO |
* |
* |
* U01 = L22T*invD2*L21 |
* U01 = L22**T*invD2*L21 |
* |
* |
CALL ZTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB, |
CALL ZTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB, |
$ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1) |
$ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1) |
Line 466
|
Line 479
|
END DO |
END DO |
ELSE |
ELSE |
* |
* |
* L11 = L11T*invD1*L11 |
* L11 = L11**T*invD1*L11 |
* |
* |
DO I=1,NNB |
DO I=1,NNB |
DO J=1,I |
DO J=1,I |
Line 480
|
Line 493
|
CUT=CUT+NNB |
CUT=CUT+NNB |
END DO |
END DO |
* |
* |
* Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P' |
* 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 |
IP=IPIV(I) |
IP=IPIV(I) |
IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP ) |
IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP ) |
IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I ) |
IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I ) |
ELSE |
ELSE |
IP=-IPIV(I) |
IP=-IPIV(I) |
IF ( I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP ) |
IF ( I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP ) |
IF ( I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I ) |
IF ( I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I ) |
I=I-1 |
I=I-1 |
ENDIF |
ENDIF |
I=I-1 |
I=I-1 |