LAPACK 3.3.0

zsytri2x.f

Go to the documentation of this file.
00001       SUBROUTINE ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.0) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2010
00007 *
00008 *  -- Written by Julie Langou of the Univ. of TN    --
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          UPLO
00012       INTEGER            INFO, LDA, N, NB
00013 *     ..
00014 *     .. Array Arguments ..
00015       INTEGER            IPIV( * )
00016       DOUBLE COMPLEX     A( LDA, * ), WORK( N+NB+1,* )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  ZSYTRI2X computes the inverse of a complex symmetric indefinite matrix
00023 *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
00024 *  ZSYTRF.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  UPLO    (input) CHARACTER*1
00030 *          Specifies whether the details of the factorization are stored
00031 *          as an upper or lower triangular matrix.
00032 *          = 'U':  Upper triangular, form is A = U*D*U**T;
00033 *          = 'L':  Lower triangular, form is A = L*D*L**T.
00034 *
00035 *  N       (input) INTEGER
00036 *          The order of the matrix A.  N >= 0.
00037 *
00038 *  A       (input/output) DOUBLE COMPLEX array, dimension (LDA,N)
00039 *          On entry, the NNB diagonal matrix D and the multipliers
00040 *          used to obtain the factor U or L as computed by ZSYTRF.
00041 *
00042 *          On exit, if INFO = 0, the (symmetric) inverse of the original
00043 *          matrix.  If UPLO = 'U', the upper triangular part of the
00044 *          inverse is formed and the part of A below the diagonal is not
00045 *          referenced; if UPLO = 'L' the lower triangular part of the
00046 *          inverse is formed and the part of A above the diagonal is
00047 *          not referenced.
00048 *
00049 *  LDA     (input) INTEGER
00050 *          The leading dimension of the array A.  LDA >= max(1,N).
00051 *
00052 *  IPIV    (input) INTEGER array, dimension (N)
00053 *          Details of the interchanges and the NNB structure of D
00054 *          as determined by ZSYTRF.
00055 *
00056 *  WORK    (workspace) DOUBLE COMPLEX array, dimension (N+NNB+1,NNB+3)
00057 *
00058 *  NB      (input) INTEGER
00059 *          Block size
00060 *
00061 *  INFO    (output) INTEGER
00062 *          = 0: successful exit
00063 *          < 0: if INFO = -i, the i-th argument had an illegal value
00064 *          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
00065 *               inverse could not be computed.
00066 *
00067 *  =====================================================================
00068 *
00069 *     .. Parameters ..
00070       DOUBLE COMPLEX     ONE, ZERO
00071       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
00072      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
00073 *     ..
00074 *     .. Local Scalars ..
00075       LOGICAL            UPPER
00076       INTEGER            I, IINFO, IP, K, CUT, NNB
00077       INTEGER            COUNT
00078       INTEGER            J, U11, INVD
00079 
00080       DOUBLE COMPLEX     AK, AKKP1, AKP1, D, T
00081       DOUBLE COMPLEX     U01_I_J, U01_IP1_J
00082       DOUBLE COMPLEX     U11_I_J, U11_IP1_J
00083 *     ..
00084 *     .. External Functions ..
00085       LOGICAL            LSAME
00086       EXTERNAL           LSAME
00087 *     ..
00088 *     .. External Subroutines ..
00089       EXTERNAL           ZSYCONV, XERBLA, ZTRTRI
00090       EXTERNAL           ZGEMM, ZTRMM, ZSYSWAPR
00091 *     ..
00092 *     .. Intrinsic Functions ..
00093       INTRINSIC          MAX
00094 *     ..
00095 *     .. Executable Statements ..
00096 *
00097 *     Test the input parameters.
00098 *
00099       INFO = 0
00100       UPPER = LSAME( UPLO, 'U' )
00101       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00102          INFO = -1
00103       ELSE IF( N.LT.0 ) THEN
00104          INFO = -2
00105       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00106          INFO = -4
00107       END IF
00108 *
00109 *     Quick return if possible
00110 *
00111 *
00112       IF( INFO.NE.0 ) THEN
00113          CALL XERBLA( 'ZSYTRI2X', -INFO )
00114          RETURN
00115       END IF
00116       IF( N.EQ.0 )
00117      $   RETURN
00118 *
00119 *     Convert A
00120 *     Workspace got Non-diag elements of D
00121 *
00122       CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
00123 *
00124 *     Check that the diagonal matrix D is nonsingular.
00125 *
00126       IF( UPPER ) THEN
00127 *
00128 *        Upper triangular storage: examine D from bottom to top
00129 *
00130          DO INFO = N, 1, -1
00131             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00132      $         RETURN
00133          END DO
00134       ELSE
00135 *
00136 *        Lower triangular storage: examine D from top to bottom.
00137 *
00138          DO INFO = 1, N
00139             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00140      $         RETURN
00141          END DO
00142       END IF
00143       INFO = 0
00144 *
00145 *  Splitting Workspace
00146 *     U01 is a block (N,NB+1) 
00147 *     The first element of U01 is in WORK(1,1)
00148 *     U11 is a block (NB+1,NB+1)
00149 *     The first element of U11 is in WORK(N+1,1)
00150       U11 = N
00151 *     INVD is a block (N,2)
00152 *     The first element of INVD is in WORK(1,INVD)
00153       INVD = NB+2
00154 
00155       IF( UPPER ) THEN
00156 *
00157 *        invA = P * inv(U')*inv(D)*inv(U)*P'.
00158 *
00159         CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
00160 *
00161 *       inv(D) and inv(D)*inv(U)
00162 * 
00163         K=1
00164         DO WHILE ( K .LE. N )
00165          IF( IPIV( K ).GT.0 ) THEN
00166 *           1 x 1 diagonal NNB
00167              WORK(K,INVD) = 1/  A( K, K )
00168              WORK(K,INVD+1) = 0
00169             K=K+1
00170          ELSE
00171 *           2 x 2 diagonal NNB
00172              T = WORK(K+1,1)
00173              AK = A( K, K ) / T
00174              AKP1 = A( K+1, K+1 ) / T
00175              AKKP1 = WORK(K+1,1)  / T
00176              D = T*( AK*AKP1-ONE )
00177              WORK(K,INVD) = AKP1 / D
00178              WORK(K+1,INVD+1) = AK / D
00179              WORK(K,INVD+1) = -AKKP1 / D  
00180              WORK(K+1,INVD) = -AKKP1 / D  
00181             K=K+2
00182          END IF
00183         END DO
00184 *
00185 *       inv(U') = (inv(U))'
00186 *
00187 *       inv(U')*inv(D)*inv(U)
00188 *
00189         CUT=N
00190         DO WHILE (CUT .GT. 0)
00191            NNB=NB
00192            IF (CUT .LE. NNB) THEN
00193               NNB=CUT
00194            ELSE
00195               COUNT = 0
00196 *             count negative elements, 
00197               DO I=CUT+1-NNB,CUT
00198                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
00199               END DO
00200 *             need a even number for a clear cut
00201               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
00202            END IF
00203 
00204            CUT=CUT-NNB
00205 *
00206 *          U01 Block 
00207 *
00208            DO I=1,CUT
00209              DO J=1,NNB
00210               WORK(I,J)=A(I,CUT+J)
00211              END DO
00212            END DO
00213 *
00214 *          U11 Block
00215 *
00216            DO I=1,NNB
00217              WORK(U11+I,I)=ONE
00218              DO J=1,I-1
00219                 WORK(U11+I,J)=ZERO
00220              END DO
00221              DO J=I+1,NNB
00222                 WORK(U11+I,J)=A(CUT+I,CUT+J)
00223              END DO
00224            END DO
00225 *
00226 *          invD*U01
00227 *
00228            I=1
00229            DO WHILE (I .LE. CUT)
00230              IF (IPIV(I) > 0) THEN
00231                 DO J=1,NNB
00232                     WORK(I,J)=WORK(I,INVD)*WORK(I,J)
00233                 END DO
00234                 I=I+1
00235              ELSE
00236                 DO J=1,NNB
00237                    U01_I_J = WORK(I,J)
00238                    U01_IP1_J = WORK(I+1,J)
00239                    WORK(I,J)=WORK(I,INVD)*U01_I_J+
00240      $                      WORK(I,INVD+1)*U01_IP1_J
00241                    WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
00242      $                      WORK(I+1,INVD+1)*U01_IP1_J
00243                 END DO
00244                 I=I+2
00245              END IF
00246            END DO
00247 *
00248 *        invD1*U11
00249 *
00250            I=1
00251            DO WHILE (I .LE. NNB)
00252              IF (IPIV(CUT+I) > 0) THEN
00253                 DO J=I,NNB
00254                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
00255                 END DO
00256                 I=I+1
00257              ELSE
00258                 DO J=I,NNB
00259                    U11_I_J = WORK(U11+I,J)
00260                    U11_IP1_J = WORK(U11+I+1,J)
00261                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
00262      $                      WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
00263                 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
00264      $                      WORK(CUT+I+1,INVD+1)*U11_IP1_J
00265                 END DO
00266                 I=I+2
00267              END IF
00268            END DO
00269 *    
00270 *       U11T*invD1*U11->U11
00271 *
00272         CALL ZTRMM('L','U','T','U',NNB, NNB,
00273      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00274 *
00275 *          U01'invD*U01->A(CUT+I,CUT+J)
00276 *
00277          CALL ZGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
00278      $              WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
00279 *
00280 *        U11 =  U11T*invD1*U11 + U01'invD*U01
00281 *
00282          DO I=1,NNB
00283             DO J=I,NNB
00284               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00285             END DO
00286          END DO
00287 *
00288 *        U01 =  U00T*invD0*U01
00289 *
00290          CALL ZTRMM('L',UPLO,'T','U',CUT, NNB,
00291      $             ONE,A,LDA,WORK,N+NB+1)
00292 
00293 *
00294 *        Update U01
00295 *
00296          DO I=1,CUT
00297            DO J=1,NNB
00298             A(I,CUT+J)=WORK(I,J)
00299            END DO
00300          END DO
00301 *
00302 *      Next Block
00303 *
00304        END DO
00305 *
00306 *        Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
00307 *  
00308             I=1
00309             DO WHILE ( I .LE. N )
00310                IF( IPIV(I) .GT. 0 ) THEN
00311                   IP=IPIV(I)
00312                  IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP )
00313                  IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I )
00314                ELSE
00315                  IP=-IPIV(I)
00316                  I=I+1
00317                  IF ( (I-1) .LT. IP) 
00318      $                  CALL ZSYSWAPR( UPLO, N, A, I-1 ,IP )
00319                  IF ( (I-1) .GT. IP) 
00320      $                  CALL ZSYSWAPR( UPLO, N, A, IP ,I-1 )
00321               ENDIF
00322                I=I+1
00323             END DO
00324       ELSE
00325 *
00326 *        LOWER...
00327 *
00328 *        invA = P * inv(U')*inv(D)*inv(U)*P'.
00329 *
00330          CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
00331 *
00332 *       inv(D) and inv(D)*inv(U)
00333 * 
00334         K=N
00335         DO WHILE ( K .GE. 1 )
00336          IF( IPIV( K ).GT.0 ) THEN
00337 *           1 x 1 diagonal NNB
00338              WORK(K,INVD) = 1/  A( K, K )
00339              WORK(K,INVD+1) = 0
00340             K=K-1
00341          ELSE
00342 *           2 x 2 diagonal NNB
00343              T = WORK(K-1,1)
00344              AK = A( K-1, K-1 ) / T
00345              AKP1 = A( K, K ) / T
00346              AKKP1 = WORK(K-1,1) / T
00347              D = T*( AK*AKP1-ONE )
00348              WORK(K-1,INVD) = AKP1 / D
00349              WORK(K,INVD) = AK / D
00350              WORK(K,INVD+1) = -AKKP1 / D  
00351              WORK(K-1,INVD+1) = -AKKP1 / D  
00352             K=K-2
00353          END IF
00354         END DO
00355 *
00356 *       inv(U') = (inv(U))'
00357 *
00358 *       inv(U')*inv(D)*inv(U)
00359 *
00360         CUT=0
00361         DO WHILE (CUT .LT. N)
00362            NNB=NB
00363            IF (CUT + NNB .GE. N) THEN
00364               NNB=N-CUT
00365            ELSE
00366               COUNT = 0
00367 *             count negative elements, 
00368               DO I=CUT+1,CUT+NNB
00369                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
00370               END DO
00371 *             need a even number for a clear cut
00372               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
00373            END IF
00374 *      L21 Block
00375            DO I=1,N-CUT-NNB
00376              DO J=1,NNB
00377               WORK(I,J)=A(CUT+NNB+I,CUT+J)
00378              END DO
00379            END DO
00380 *     L11 Block
00381            DO I=1,NNB
00382              WORK(U11+I,I)=ONE
00383              DO J=I+1,NNB
00384                 WORK(U11+I,J)=ZERO
00385              END DO
00386              DO J=1,I-1
00387                 WORK(U11+I,J)=A(CUT+I,CUT+J)
00388              END DO
00389            END DO
00390 *
00391 *          invD*L21
00392 *
00393            I=N-CUT-NNB
00394            DO WHILE (I .GE. 1)
00395              IF (IPIV(CUT+NNB+I) > 0) THEN
00396                 DO J=1,NNB
00397                     WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
00398                 END DO
00399                 I=I-1
00400              ELSE
00401                 DO J=1,NNB
00402                    U01_I_J = WORK(I,J)
00403                    U01_IP1_J = WORK(I-1,J)
00404                    WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
00405      $                        WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
00406                    WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
00407      $                        WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
00408                 END DO
00409                 I=I-2
00410              END IF
00411            END DO
00412 *
00413 *        invD1*L11
00414 *
00415            I=NNB
00416            DO WHILE (I .GE. 1)
00417              IF (IPIV(CUT+I) > 0) THEN
00418                 DO J=1,NNB
00419                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
00420                 END DO
00421                 I=I-1
00422              ELSE
00423                 DO J=1,NNB
00424                    U11_I_J = WORK(U11+I,J)
00425                    U11_IP1_J = WORK(U11+I-1,J)
00426                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
00427      $                      WORK(CUT+I,INVD+1)*U11_IP1_J
00428                 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
00429      $                      WORK(CUT+I-1,INVD)*U11_IP1_J
00430                 END DO
00431                 I=I-2
00432              END IF
00433            END DO
00434 *    
00435 *       L11T*invD1*L11->L11
00436 *
00437         CALL ZTRMM('L',UPLO,'T','U',NNB, NNB,
00438      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00439 
00440         IF ( (CUT+NNB) .LT. N ) THEN
00441 *
00442 *          L21T*invD2*L21->A(CUT+I,CUT+J)
00443 *
00444          CALL ZGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
00445      $             ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
00446        
00447 *
00448 *        L11 =  L11T*invD1*L11 + U01'invD*U01
00449 *
00450          DO I=1,NNB
00451             DO J=1,I
00452               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00453             END DO
00454          END DO
00455 *
00456 *        U01 =  L22T*invD2*L21
00457 *
00458          CALL ZTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
00459      $             ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
00460 
00461 *      Update L21
00462          DO I=1,N-CUT-NNB
00463            DO J=1,NNB
00464               A(CUT+NNB+I,CUT+J)=WORK(I,J)
00465            END DO
00466          END DO
00467        ELSE
00468 *
00469 *        L11 =  L11T*invD1*L11
00470 *
00471          DO I=1,NNB
00472             DO J=1,I
00473               A(CUT+I,CUT+J)=WORK(U11+I,J)
00474             END DO
00475          END DO
00476        END IF
00477 *
00478 *      Next Block
00479 *
00480            CUT=CUT+NNB
00481        END DO
00482 *
00483 *        Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
00484 * 
00485             I=N
00486             DO WHILE ( I .GE. 1 )
00487                IF( IPIV(I) .GT. 0 ) THEN
00488                   IP=IPIV(I)
00489                  IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP  )
00490                  IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I )
00491                ELSE
00492                  IP=-IPIV(I)
00493                  IF ( I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP )
00494                  IF ( I .GT. IP) CALL ZSYSWAPR(  UPLO, N, A, IP ,I )
00495                  I=I-1
00496                ENDIF
00497                I=I-1
00498             END DO
00499       END IF
00500 *
00501       RETURN
00502 *
00503 *     End of ZSYTRI2X
00504 *
00505       END
00506 
 All Files Functions