LAPACK 3.3.0

dsytri2x.f

Go to the documentation of this file.
00001       SUBROUTINE DSYTRI2X( 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 PRECISION   A( LDA, * ), WORK( N+NB+1,* )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  DSYTRI2X computes the inverse of a real symmetric indefinite matrix
00023 *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
00024 *  DSYTRF.
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 PRECISION 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 DSYTRF.
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 DSYTRF.
00055 *
00056 *  WORK    (workspace) DOUBLE PRECISION 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 PRECISION   ONE, ZERO
00071       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00072 *     ..
00073 *     .. Local Scalars ..
00074       LOGICAL            UPPER
00075       INTEGER            I, IINFO, IP, K, CUT, NNB
00076       INTEGER            COUNT
00077       INTEGER            J, U11, INVD
00078 
00079       DOUBLE PRECISION   AK, AKKP1, AKP1, D, T
00080       DOUBLE PRECISION   U01_I_J, U01_IP1_J
00081       DOUBLE PRECISION   U11_I_J, U11_IP1_J
00082 *     ..
00083 *     .. External Functions ..
00084       LOGICAL            LSAME
00085       EXTERNAL           LSAME
00086 *     ..
00087 *     .. External Subroutines ..
00088       EXTERNAL           DSYCONV, XERBLA, DTRTRI
00089       EXTERNAL           DGEMM, DTRMM, DSYSWAPR
00090 *     ..
00091 *     .. Intrinsic Functions ..
00092       INTRINSIC          MAX
00093 *     ..
00094 *     .. Executable Statements ..
00095 *
00096 *     Test the input parameters.
00097 *
00098       INFO = 0
00099       UPPER = LSAME( UPLO, 'U' )
00100       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00101          INFO = -1
00102       ELSE IF( N.LT.0 ) THEN
00103          INFO = -2
00104       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00105          INFO = -4
00106       END IF
00107 *
00108 *     Quick return if possible
00109 *
00110 *
00111       IF( INFO.NE.0 ) THEN
00112          CALL XERBLA( 'DSYTRI2X', -INFO )
00113          RETURN
00114       END IF
00115       IF( N.EQ.0 )
00116      $   RETURN
00117 *
00118 *     Convert A
00119 *     Workspace got Non-diag elements of D
00120 *
00121       CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
00122 *
00123 *     Check that the diagonal matrix D is nonsingular.
00124 *
00125       IF( UPPER ) THEN
00126 *
00127 *        Upper triangular storage: examine D from bottom to top
00128 *
00129          DO INFO = N, 1, -1
00130             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00131      $         RETURN
00132          END DO
00133       ELSE
00134 *
00135 *        Lower triangular storage: examine D from top to bottom.
00136 *
00137          DO INFO = 1, N
00138             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
00139      $         RETURN
00140          END DO
00141       END IF
00142       INFO = 0
00143 *
00144 *  Splitting Workspace
00145 *     U01 is a block (N,NB+1) 
00146 *     The first element of U01 is in WORK(1,1)
00147 *     U11 is a block (NB+1,NB+1)
00148 *     The first element of U11 is in WORK(N+1,1)
00149       U11 = N
00150 *     INVD is a block (N,2)
00151 *     The first element of INVD is in WORK(1,INVD)
00152       INVD = NB+2
00153 
00154       IF( UPPER ) THEN
00155 *
00156 *        invA = P * inv(U')*inv(D)*inv(U)*P'.
00157 *
00158         CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
00159 *
00160 *       inv(D) and inv(D)*inv(U)
00161 * 
00162         K=1
00163         DO WHILE ( K .LE. N )
00164          IF( IPIV( K ).GT.0 ) THEN
00165 *           1 x 1 diagonal NNB
00166              WORK(K,INVD) = ONE /  A( K, K )
00167              WORK(K,INVD+1) = 0
00168             K=K+1
00169          ELSE
00170 *           2 x 2 diagonal NNB
00171              T = WORK(K+1,1)
00172              AK = A( K, K ) / T
00173              AKP1 = A( K+1, K+1 ) / T
00174              AKKP1 = WORK(K+1,1)  / T
00175              D = T*( AK*AKP1-ONE )
00176              WORK(K,INVD) = AKP1 / D
00177              WORK(K+1,INVD+1) = AK / D
00178              WORK(K,INVD+1) = -AKKP1 / D  
00179              WORK(K+1,INVD) = -AKKP1 / D  
00180             K=K+2
00181          END IF
00182         END DO
00183 *
00184 *       inv(U') = (inv(U))'
00185 *
00186 *       inv(U')*inv(D)*inv(U)
00187 *
00188         CUT=N
00189         DO WHILE (CUT .GT. 0)
00190            NNB=NB
00191            IF (CUT .LE. NNB) THEN
00192               NNB=CUT
00193            ELSE
00194               COUNT = 0
00195 *             count negative elements, 
00196               DO I=CUT+1-NNB,CUT
00197                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
00198               END DO
00199 *             need a even number for a clear cut
00200               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
00201            END IF
00202 
00203            CUT=CUT-NNB
00204 *
00205 *          U01 Block 
00206 *
00207            DO I=1,CUT
00208              DO J=1,NNB
00209               WORK(I,J)=A(I,CUT+J)
00210              END DO
00211            END DO
00212 *
00213 *          U11 Block
00214 *
00215            DO I=1,NNB
00216              WORK(U11+I,I)=ONE
00217              DO J=1,I-1
00218                 WORK(U11+I,J)=ZERO
00219              END DO
00220              DO J=I+1,NNB
00221                 WORK(U11+I,J)=A(CUT+I,CUT+J)
00222              END DO
00223            END DO
00224 *
00225 *          invD*U01
00226 *
00227            I=1
00228            DO WHILE (I .LE. CUT)
00229              IF (IPIV(I) > 0) THEN
00230                 DO J=1,NNB
00231                     WORK(I,J)=WORK(I,INVD)*WORK(I,J)
00232                 END DO
00233                 I=I+1
00234              ELSE
00235                 DO J=1,NNB
00236                    U01_I_J = WORK(I,J)
00237                    U01_IP1_J = WORK(I+1,J)
00238                    WORK(I,J)=WORK(I,INVD)*U01_I_J+
00239      $                      WORK(I,INVD+1)*U01_IP1_J
00240                    WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
00241      $                      WORK(I+1,INVD+1)*U01_IP1_J
00242                 END DO
00243                 I=I+2
00244              END IF
00245            END DO
00246 *
00247 *        invD1*U11
00248 *
00249            I=1
00250            DO WHILE (I .LE. NNB)
00251              IF (IPIV(CUT+I) > 0) THEN
00252                 DO J=I,NNB
00253                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
00254                 END DO
00255                 I=I+1
00256              ELSE
00257                 DO J=I,NNB
00258                    U11_I_J = WORK(U11+I,J)
00259                    U11_IP1_J = WORK(U11+I+1,J)
00260                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
00261      $                      WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
00262                 WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
00263      $                      WORK(CUT+I+1,INVD+1)*U11_IP1_J
00264                 END DO
00265                 I=I+2
00266              END IF
00267            END DO
00268 *    
00269 *       U11T*invD1*U11->U11
00270 *
00271         CALL DTRMM('L','U','T','U',NNB, NNB,
00272      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00273 *
00274 *          U01'invD*U01->A(CUT+I,CUT+J)
00275 *
00276          CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
00277      $              WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
00278 *
00279 *        U11 =  U11T*invD1*U11 + U01'invD*U01
00280 *
00281          DO I=1,NNB
00282             DO J=I,NNB
00283               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00284             END DO
00285          END DO
00286 *
00287 *        U01 =  U00T*invD0*U01
00288 *
00289          CALL DTRMM('L',UPLO,'T','U',CUT, NNB,
00290      $             ONE,A,LDA,WORK,N+NB+1)
00291 
00292 *
00293 *        Update U01
00294 *
00295          DO I=1,CUT
00296            DO J=1,NNB
00297             A(I,CUT+J)=WORK(I,J)
00298            END DO
00299          END DO
00300 *
00301 *      Next Block
00302 *
00303        END DO
00304 *
00305 *        Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
00306 *  
00307             I=1
00308             DO WHILE ( I .LE. N )
00309                IF( IPIV(I) .GT. 0 ) THEN
00310                   IP=IPIV(I)
00311                  IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP )
00312                  IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I )
00313                ELSE
00314                  IP=-IPIV(I)
00315                  I=I+1
00316                  IF ( (I-1) .LT. IP) 
00317      $                  CALL DSYSWAPR( UPLO, N, A, I-1 ,IP )
00318                  IF ( (I-1) .GT. IP) 
00319      $                  CALL DSYSWAPR( UPLO, N, A, IP ,I-1 )
00320               ENDIF
00321                I=I+1
00322             END DO
00323       ELSE
00324 *
00325 *        LOWER...
00326 *
00327 *        invA = P * inv(U')*inv(D)*inv(U)*P'.
00328 *
00329          CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO )
00330 *
00331 *       inv(D) and inv(D)*inv(U)
00332 * 
00333         K=N
00334         DO WHILE ( K .GE. 1 )
00335          IF( IPIV( K ).GT.0 ) THEN
00336 *           1 x 1 diagonal NNB
00337              WORK(K,INVD) = ONE /  A( K, K )
00338              WORK(K,INVD+1) = 0
00339             K=K-1
00340          ELSE
00341 *           2 x 2 diagonal NNB
00342              T = WORK(K-1,1)
00343              AK = A( K-1, K-1 ) / T
00344              AKP1 = A( K, K ) / T
00345              AKKP1 = WORK(K-1,1) / T
00346              D = T*( AK*AKP1-ONE )
00347              WORK(K-1,INVD) = AKP1 / D
00348              WORK(K,INVD) = AK / D
00349              WORK(K,INVD+1) = -AKKP1 / D  
00350              WORK(K-1,INVD+1) = -AKKP1 / D  
00351             K=K-2
00352          END IF
00353         END DO
00354 *
00355 *       inv(U') = (inv(U))'
00356 *
00357 *       inv(U')*inv(D)*inv(U)
00358 *
00359         CUT=0
00360         DO WHILE (CUT .LT. N)
00361            NNB=NB
00362            IF (CUT + NNB .GT. N) THEN
00363               NNB=N-CUT
00364            ELSE
00365               COUNT = 0
00366 *             count negative elements, 
00367               DO I=CUT+1,CUT+NNB
00368                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
00369               END DO
00370 *             need a even number for a clear cut
00371               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
00372            END IF
00373 *     L21 Block
00374            DO I=1,N-CUT-NNB
00375              DO J=1,NNB
00376               WORK(I,J)=A(CUT+NNB+I,CUT+J)
00377              END DO
00378            END DO
00379 *     L11 Block
00380            DO I=1,NNB
00381              WORK(U11+I,I)=ONE
00382              DO J=I+1,NNB
00383                 WORK(U11+I,J)=ZERO
00384              END DO
00385              DO J=1,I-1
00386                 WORK(U11+I,J)=A(CUT+I,CUT+J)
00387              END DO
00388            END DO
00389 *
00390 *          invD*L21
00391 *
00392            I=N-CUT-NNB
00393            DO WHILE (I .GE. 1)
00394              IF (IPIV(CUT+NNB+I) > 0) THEN
00395                 DO J=1,NNB
00396                     WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
00397                 END DO
00398                 I=I-1
00399              ELSE
00400                 DO J=1,NNB
00401                    U01_I_J = WORK(I,J)
00402                    U01_IP1_J = WORK(I-1,J)
00403                    WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
00404      $                        WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
00405                    WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
00406      $                        WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
00407                 END DO
00408                 I=I-2
00409              END IF
00410            END DO
00411 *
00412 *        invD1*L11
00413 *
00414            I=NNB
00415            DO WHILE (I .GE. 1)
00416              IF (IPIV(CUT+I) > 0) THEN
00417                 DO J=1,NNB
00418                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
00419                 END DO
00420                 I=I-1
00421              ELSE
00422                 DO J=1,NNB
00423                    U11_I_J = WORK(U11+I,J)
00424                    U11_IP1_J = WORK(U11+I-1,J)
00425                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
00426      $                      WORK(CUT+I,INVD+1)*U11_IP1_J
00427                 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
00428      $                      WORK(CUT+I-1,INVD)*U11_IP1_J
00429                 END DO
00430                 I=I-2
00431              END IF
00432            END DO
00433 *    
00434 *       L11T*invD1*L11->L11
00435 *
00436         CALL DTRMM('L',UPLO,'T','U',NNB, NNB,
00437      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00438 
00439         IF ( (CUT+NNB) .LT. N ) THEN
00440 *
00441 *          L21T*invD2*L21->A(CUT+I,CUT+J)
00442 *
00443          CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
00444      $             ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
00445        
00446 *
00447 *        L11 =  L11T*invD1*L11 + U01'invD*U01
00448 *
00449          DO I=1,NNB
00450             DO J=1,I
00451               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00452             END DO
00453          END DO
00454 *
00455 *        L01 =  L22T*invD2*L21
00456 *
00457          CALL DTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
00458      $             ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
00459 *
00460 *      Update L21
00461 *
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 
00468        ELSE
00469 *
00470 *        L11 =  L11T*invD1*L11
00471 *
00472          DO I=1,NNB
00473             DO J=1,I
00474               A(CUT+I,CUT+J)=WORK(U11+I,J)
00475             END DO
00476          END DO
00477        END IF
00478 *
00479 *      Next Block
00480 *
00481            CUT=CUT+NNB
00482        END DO
00483 *
00484 *        Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
00485 * 
00486             I=N
00487             DO WHILE ( I .GE. 1 )
00488                IF( IPIV(I) .GT. 0 ) THEN
00489                   IP=IPIV(I)
00490                  IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP  )
00491                  IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I )
00492                ELSE
00493                  IP=-IPIV(I)
00494                  IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP )
00495                  IF ( I .GT. IP) CALL DSYSWAPR(  UPLO, N, A, IP ,I )
00496                  I=I-1
00497                ENDIF
00498                I=I-1
00499             END DO
00500       END IF
00501 *
00502       RETURN
00503 *
00504 *     End of DSYTRI2X
00505 *
00506       END
00507 
 All Files Functions