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