LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dsytri2x.f
Go to the documentation of this file.
00001 *> \brief \b DSYTRI2X
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DSYTRI2X + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytri2x.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytri2x.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytri2x.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DSYTRI2X( 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 *       DOUBLE PRECISION   A( LDA, * ), WORK( N+NB+1,* )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> DSYTRI2X computes the inverse of a real symmetric indefinite matrix
00039 *> A using the factorization A = U*D*U**T or A = L*D*L**T computed by
00040 *> DSYTRF.
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 DOUBLE PRECISION 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 DSYTRF.
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 DSYTRF.
00086 *> \endverbatim
00087 *>
00088 *> \param[out] WORK
00089 *> \verbatim
00090 *>          WORK is DOUBLE PRECISION 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 doubleSYcomputational
00119 *
00120 *  =====================================================================
00121       SUBROUTINE DSYTRI2X( 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       DOUBLE PRECISION   A( LDA, * ), WORK( N+NB+1,* )
00135 *     ..
00136 *
00137 *  =====================================================================
00138 *
00139 *     .. Parameters ..
00140       DOUBLE PRECISION   ONE, ZERO
00141       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00142 *     ..
00143 *     .. Local Scalars ..
00144       LOGICAL            UPPER
00145       INTEGER            I, IINFO, IP, K, CUT, NNB
00146       INTEGER            COUNT
00147       INTEGER            J, U11, INVD
00148 
00149       DOUBLE PRECISION   AK, AKKP1, AKP1, D, T
00150       DOUBLE PRECISION   U01_I_J, U01_IP1_J
00151       DOUBLE PRECISION   U11_I_J, U11_IP1_J
00152 *     ..
00153 *     .. External Functions ..
00154       LOGICAL            LSAME
00155       EXTERNAL           LSAME
00156 *     ..
00157 *     .. External Subroutines ..
00158       EXTERNAL           DSYCONV, XERBLA, DTRTRI
00159       EXTERNAL           DGEMM, DTRMM, DSYSWAPR
00160 *     ..
00161 *     .. Intrinsic Functions ..
00162       INTRINSIC          MAX
00163 *     ..
00164 *     .. Executable Statements ..
00165 *
00166 *     Test the input parameters.
00167 *
00168       INFO = 0
00169       UPPER = LSAME( UPLO, 'U' )
00170       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00171          INFO = -1
00172       ELSE IF( N.LT.0 ) THEN
00173          INFO = -2
00174       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00175          INFO = -4
00176       END IF
00177 *
00178 *     Quick return if possible
00179 *
00180 *
00181       IF( INFO.NE.0 ) THEN
00182          CALL XERBLA( 'DSYTRI2X', -INFO )
00183          RETURN
00184       END IF
00185       IF( N.EQ.0 )
00186      $   RETURN
00187 *
00188 *     Convert A
00189 *     Workspace got Non-diag elements of D
00190 *
00191       CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
00192 *
00193 *     Check that the diagonal matrix D is nonsingular.
00194 *
00195       IF( UPPER ) THEN
00196 *
00197 *        Upper triangular storage: examine D from bottom to top
00198 *
00199          DO INFO = N, 1, -1
00200             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00201      $         RETURN
00202          END DO
00203       ELSE
00204 *
00205 *        Lower triangular storage: examine D from top to bottom.
00206 *
00207          DO INFO = 1, N
00208             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00209      $         RETURN
00210          END DO
00211       END IF
00212       INFO = 0
00213 *
00214 *  Splitting Workspace
00215 *     U01 is a block (N,NB+1) 
00216 *     The first element of U01 is in WORK(1,1)
00217 *     U11 is a block (NB+1,NB+1)
00218 *     The first element of U11 is in WORK(N+1,1)
00219       U11 = N
00220 *     INVD is a block (N,2)
00221 *     The first element of INVD is in WORK(1,INVD)
00222       INVD = NB+2
00223 
00224       IF( UPPER ) THEN
00225 *
00226 *        invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
00227 *
00228         CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
00229 *
00230 *       inv(D) and inv(D)*inv(U)
00231 * 
00232         K=1
00233         DO WHILE ( K .LE. N )
00234          IF( IPIV( K ).GT.0 ) THEN
00235 *           1 x 1 diagonal NNB
00236              WORK(K,INVD) = ONE /  A( K, K )
00237              WORK(K,INVD+1) = 0
00238             K=K+1
00239          ELSE
00240 *           2 x 2 diagonal NNB
00241              T = WORK(K+1,1)
00242              AK = A( K, K ) / T
00243              AKP1 = A( K+1, K+1 ) / T
00244              AKKP1 = WORK(K+1,1)  / T
00245              D = T*( AK*AKP1-ONE )
00246              WORK(K,INVD) = AKP1 / D
00247              WORK(K+1,INVD+1) = AK / D
00248              WORK(K,INVD+1) = -AKKP1 / D  
00249              WORK(K+1,INVD) = -AKKP1 / D  
00250             K=K+2
00251          END IF
00252         END DO
00253 *
00254 *       inv(U**T) = (inv(U))**T
00255 *
00256 *       inv(U**T)*inv(D)*inv(U)
00257 *
00258         CUT=N
00259         DO WHILE (CUT .GT. 0)
00260            NNB=NB
00261            IF (CUT .LE. NNB) THEN
00262               NNB=CUT
00263            ELSE
00264               COUNT = 0
00265 *             count negative elements, 
00266               DO I=CUT+1-NNB,CUT
00267                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
00268               END DO
00269 *             need a even number for a clear cut
00270               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
00271            END IF
00272 
00273            CUT=CUT-NNB
00274 *
00275 *          U01 Block 
00276 *
00277            DO I=1,CUT
00278              DO J=1,NNB
00279               WORK(I,J)=A(I,CUT+J)
00280              END DO
00281            END DO
00282 *
00283 *          U11 Block
00284 *
00285            DO I=1,NNB
00286              WORK(U11+I,I)=ONE
00287              DO J=1,I-1
00288                 WORK(U11+I,J)=ZERO
00289              END DO
00290              DO J=I+1,NNB
00291                 WORK(U11+I,J)=A(CUT+I,CUT+J)
00292              END DO
00293            END DO
00294 *
00295 *          invD*U01
00296 *
00297            I=1
00298            DO WHILE (I .LE. CUT)
00299              IF (IPIV(I) > 0) THEN
00300                 DO J=1,NNB
00301                     WORK(I,J)=WORK(I,INVD)*WORK(I,J)
00302                 END DO
00303                 I=I+1
00304              ELSE
00305                 DO J=1,NNB
00306                    U01_I_J = WORK(I,J)
00307                    U01_IP1_J = WORK(I+1,J)
00308                    WORK(I,J)=WORK(I,INVD)*U01_I_J+
00309      $                      WORK(I,INVD+1)*U01_IP1_J
00310                    WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
00311      $                      WORK(I+1,INVD+1)*U01_IP1_J
00312                 END DO
00313                 I=I+2
00314              END IF
00315            END DO
00316 *
00317 *        invD1*U11
00318 *
00319            I=1
00320            DO WHILE (I .LE. NNB)
00321              IF (IPIV(CUT+I) > 0) THEN
00322                 DO J=I,NNB
00323                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
00324                 END DO
00325                 I=I+1
00326              ELSE
00327                 DO J=I,NNB
00328                    U11_I_J = WORK(U11+I,J)
00329                    U11_IP1_J = WORK(U11+I+1,J)
00330                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
00331      $                      WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
00332                 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
00333      $                      WORK(CUT+I+1,INVD+1)*U11_IP1_J
00334                 END DO
00335                 I=I+2
00336              END IF
00337            END DO
00338 *    
00339 *       U11**T*invD1*U11->U11
00340 *
00341         CALL DTRMM('L','U','T','U',NNB, NNB,
00342      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00343 *
00344          DO I=1,NNB
00345             DO J=I,NNB
00346               A(CUT+I,CUT+J)=WORK(U11+I,J)
00347             END DO
00348          END DO         
00349 *
00350 *          U01**T*invD*U01->A(CUT+I,CUT+J)
00351 *
00352          CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
00353      $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
00354         
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 DTRMM('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 DSYSWAPR( UPLO, N, A, LDA, I ,IP )
00389                  IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
00390                ELSE
00391                  IP=-IPIV(I)
00392                  I=I+1
00393                  IF ( (I-1) .LT. IP) 
00394      $                  CALL DSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
00395                  IF ( (I-1) .GT. IP) 
00396      $                  CALL DSYSWAPR( 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 DTRTRI( 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) = ONE /  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 .GT. 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 DTRMM('L',UPLO,'T','U',NNB, NNB,
00514      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00515 
00516 *
00517          DO I=1,NNB
00518             DO J=1,I
00519               A(CUT+I,CUT+J)=WORK(U11+I,J)
00520             END DO
00521          END DO
00522 *
00523         IF ( (CUT+NNB) .LT. N ) THEN
00524 *
00525 *          L21**T*invD2*L21->A(CUT+I,CUT+J)
00526 *
00527          CALL DGEMM('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 *        L01 =  L22**T*invD2*L21
00540 *
00541          CALL DTRMM('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 *
00546          DO I=1,N-CUT-NNB
00547            DO J=1,NNB
00548               A(CUT+NNB+I,CUT+J)=WORK(I,J)
00549            END DO
00550          END DO
00551 
00552        ELSE
00553 *
00554 *        L11 =  L11**T*invD1*L11
00555 *
00556          DO I=1,NNB
00557             DO J=1,I
00558               A(CUT+I,CUT+J)=WORK(U11+I,J)
00559             END DO
00560          END DO
00561        END IF
00562 *
00563 *      Next Block
00564 *
00565            CUT=CUT+NNB
00566        END DO
00567 *
00568 *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
00569 * 
00570             I=N
00571             DO WHILE ( I .GE. 1 )
00572                IF( IPIV(I) .GT. 0 ) THEN
00573                   IP=IPIV(I)
00574                  IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP  )
00575                  IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
00576                ELSE
00577                  IP=-IPIV(I)
00578                  IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
00579                  IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP, I )
00580                  I=I-1
00581                ENDIF
00582                I=I-1
00583             END DO
00584       END IF
00585 *
00586       RETURN
00587 *
00588 *     End of DSYTRI2X
00589 *
00590       END
00591 
 All Files Functions