LAPACK 3.3.1
Linear Algebra PACKage

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.1) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *  -- April 2011                                                      --
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**T)*inv(D)*inv(U)*P**T.
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**T) = (inv(U))**T
00186 *
00187 *       inv(U**T)*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 *       U11**T*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          DO I=1,NNB
00276             DO J=I,NNB
00277               A(CUT+I,CUT+J)=WORK(U11+I,J)
00278             END DO
00279          END DO         
00280 *
00281 *          U01**T*invD*U01->A(CUT+I,CUT+J)
00282 *
00283          CALL ZGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
00284      $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
00285 *
00286 *        U11 =  U11**T*invD1*U11 + U01**T*invD*U01
00287 *
00288          DO I=1,NNB
00289             DO J=I,NNB
00290               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00291             END DO
00292          END DO
00293 *
00294 *        U01 =  U00**T*invD0*U01
00295 *
00296          CALL ZTRMM('L',UPLO,'T','U',CUT, NNB,
00297      $             ONE,A,LDA,WORK,N+NB+1)
00298 
00299 *
00300 *        Update U01
00301 *
00302          DO I=1,CUT
00303            DO J=1,NNB
00304             A(I,CUT+J)=WORK(I,J)
00305            END DO
00306          END DO
00307 *
00308 *      Next Block
00309 *
00310        END DO
00311 *
00312 *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
00313 *  
00314             I=1
00315             DO WHILE ( I .LE. N )
00316                IF( IPIV(I) .GT. 0 ) THEN
00317                   IP=IPIV(I)
00318                  IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP )
00319                  IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
00320                ELSE
00321                  IP=-IPIV(I)
00322                  I=I+1
00323                  IF ( (I-1) .LT. IP) 
00324      $                  CALL ZSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
00325                  IF ( (I-1) .GT. IP) 
00326      $                  CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I-1 )
00327               ENDIF
00328                I=I+1
00329             END DO
00330       ELSE
00331 *
00332 *        LOWER...
00333 *
00334 *        invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
00335 *
00336          CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
00337 *
00338 *       inv(D) and inv(D)*inv(U)
00339 * 
00340         K=N
00341         DO WHILE ( K .GE. 1 )
00342          IF( IPIV( K ).GT.0 ) THEN
00343 *           1 x 1 diagonal NNB
00344              WORK(K,INVD) = 1/  A( K, K )
00345              WORK(K,INVD+1) = 0
00346             K=K-1
00347          ELSE
00348 *           2 x 2 diagonal NNB
00349              T = WORK(K-1,1)
00350              AK = A( K-1, K-1 ) / T
00351              AKP1 = A( K, K ) / T
00352              AKKP1 = WORK(K-1,1) / T
00353              D = T*( AK*AKP1-ONE )
00354              WORK(K-1,INVD) = AKP1 / D
00355              WORK(K,INVD) = AK / D
00356              WORK(K,INVD+1) = -AKKP1 / D  
00357              WORK(K-1,INVD+1) = -AKKP1 / D  
00358             K=K-2
00359          END IF
00360         END DO
00361 *
00362 *       inv(U**T) = (inv(U))**T
00363 *
00364 *       inv(U**T)*inv(D)*inv(U)
00365 *
00366         CUT=0
00367         DO WHILE (CUT .LT. N)
00368            NNB=NB
00369            IF (CUT + NNB .GE. N) THEN
00370               NNB=N-CUT
00371            ELSE
00372               COUNT = 0
00373 *             count negative elements, 
00374               DO I=CUT+1,CUT+NNB
00375                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
00376               END DO
00377 *             need a even number for a clear cut
00378               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
00379            END IF
00380 *      L21 Block
00381            DO I=1,N-CUT-NNB
00382              DO J=1,NNB
00383               WORK(I,J)=A(CUT+NNB+I,CUT+J)
00384              END DO
00385            END DO
00386 *     L11 Block
00387            DO I=1,NNB
00388              WORK(U11+I,I)=ONE
00389              DO J=I+1,NNB
00390                 WORK(U11+I,J)=ZERO
00391              END DO
00392              DO J=1,I-1
00393                 WORK(U11+I,J)=A(CUT+I,CUT+J)
00394              END DO
00395            END DO
00396 *
00397 *          invD*L21
00398 *
00399            I=N-CUT-NNB
00400            DO WHILE (I .GE. 1)
00401              IF (IPIV(CUT+NNB+I) > 0) THEN
00402                 DO J=1,NNB
00403                     WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
00404                 END DO
00405                 I=I-1
00406              ELSE
00407                 DO J=1,NNB
00408                    U01_I_J = WORK(I,J)
00409                    U01_IP1_J = WORK(I-1,J)
00410                    WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
00411      $                        WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
00412                    WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
00413      $                        WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
00414                 END DO
00415                 I=I-2
00416              END IF
00417            END DO
00418 *
00419 *        invD1*L11
00420 *
00421            I=NNB
00422            DO WHILE (I .GE. 1)
00423              IF (IPIV(CUT+I) > 0) THEN
00424                 DO J=1,NNB
00425                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
00426                 END DO
00427                 I=I-1
00428              ELSE
00429                 DO J=1,NNB
00430                    U11_I_J = WORK(U11+I,J)
00431                    U11_IP1_J = WORK(U11+I-1,J)
00432                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
00433      $                      WORK(CUT+I,INVD+1)*U11_IP1_J
00434                 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
00435      $                      WORK(CUT+I-1,INVD)*U11_IP1_J
00436                 END DO
00437                 I=I-2
00438              END IF
00439            END DO
00440 *    
00441 *       L11**T*invD1*L11->L11
00442 *
00443         CALL ZTRMM('L',UPLO,'T','U',NNB, NNB,
00444      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00445 *
00446          DO I=1,NNB
00447             DO J=1,I
00448               A(CUT+I,CUT+J)=WORK(U11+I,J)
00449             END DO
00450          END DO
00451 *
00452 
00453         IF ( (CUT+NNB) .LT. N ) THEN
00454 *
00455 *          L21**T*invD2*L21->A(CUT+I,CUT+J)
00456 *
00457          CALL ZGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
00458      $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
00459        
00460 *
00461 *        L11 =  L11**T*invD1*L11 + U01**T*invD*U01
00462 *
00463          DO I=1,NNB
00464             DO J=1,I
00465               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00466             END DO
00467          END DO
00468 *
00469 *        U01 =  L22**T*invD2*L21
00470 *
00471          CALL ZTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
00472      $             ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
00473 
00474 *      Update L21
00475          DO I=1,N-CUT-NNB
00476            DO J=1,NNB
00477               A(CUT+NNB+I,CUT+J)=WORK(I,J)
00478            END DO
00479          END DO
00480        ELSE
00481 *
00482 *        L11 =  L11**T*invD1*L11
00483 *
00484          DO I=1,NNB
00485             DO J=1,I
00486               A(CUT+I,CUT+J)=WORK(U11+I,J)
00487             END DO
00488          END DO
00489        END IF
00490 *
00491 *      Next Block
00492 *
00493            CUT=CUT+NNB
00494        END DO
00495 *
00496 *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
00497 * 
00498             I=N
00499             DO WHILE ( I .GE. 1 )
00500                IF( IPIV(I) .GT. 0 ) THEN
00501                   IP=IPIV(I)
00502                  IF (I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP  )
00503                  IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
00504                ELSE
00505                  IP=-IPIV(I)
00506                  IF ( I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, I ,IP )
00507                  IF ( I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, LDA, IP ,I )
00508                  I=I-1
00509                ENDIF
00510                I=I-1
00511             END DO
00512       END IF
00513 *
00514       RETURN
00515 *
00516 *     End of ZSYTRI2X
00517 *
00518       END
00519 
 All Files Functions