LAPACK 3.3.1
Linear Algebra PACKage

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.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 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**T)*inv(D)*inv(U)*P**T.
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**T) = (inv(U))**T
00185 *
00186 *       inv(U**T)*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 *       U11**T*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          DO I=1,NNB
00275             DO J=I,NNB
00276               A(CUT+I,CUT+J)=WORK(U11+I,J)
00277             END DO
00278          END DO         
00279 *
00280 *          U01**T*invD*U01->A(CUT+I,CUT+J)
00281 *
00282          CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA,
00283      $              WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
00284         
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 DTRMM('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 DSYSWAPR( UPLO, N, A, LDA, I ,IP )
00319                  IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
00320                ELSE
00321                  IP=-IPIV(I)
00322                  I=I+1
00323                  IF ( (I-1) .LT. IP) 
00324      $                  CALL DSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
00325                  IF ( (I-1) .GT. IP) 
00326      $                  CALL DSYSWAPR( 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 DTRTRI( 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) = ONE /  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 .GT. 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 DTRMM('L',UPLO,'T','U',NNB, NNB,
00444      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00445 
00446 *
00447          DO I=1,NNB
00448             DO J=1,I
00449               A(CUT+I,CUT+J)=WORK(U11+I,J)
00450             END DO
00451          END DO
00452 *
00453         IF ( (CUT+NNB) .LT. N ) THEN
00454 *
00455 *          L21**T*invD2*L21->A(CUT+I,CUT+J)
00456 *
00457          CALL DGEMM('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 *        L01 =  L22**T*invD2*L21
00470 *
00471          CALL DTRMM('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 *
00476          DO I=1,N-CUT-NNB
00477            DO J=1,NNB
00478               A(CUT+NNB+I,CUT+J)=WORK(I,J)
00479            END DO
00480          END DO
00481 
00482        ELSE
00483 *
00484 *        L11 =  L11**T*invD1*L11
00485 *
00486          DO I=1,NNB
00487             DO J=1,I
00488               A(CUT+I,CUT+J)=WORK(U11+I,J)
00489             END DO
00490          END DO
00491        END IF
00492 *
00493 *      Next Block
00494 *
00495            CUT=CUT+NNB
00496        END DO
00497 *
00498 *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
00499 * 
00500             I=N
00501             DO WHILE ( I .GE. 1 )
00502                IF( IPIV(I) .GT. 0 ) THEN
00503                   IP=IPIV(I)
00504                  IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP  )
00505                  IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP ,I )
00506                ELSE
00507                  IP=-IPIV(I)
00508                  IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, I ,IP )
00509                  IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, LDA, IP, I )
00510                  I=I-1
00511                ENDIF
00512                I=I-1
00513             END DO
00514       END IF
00515 *
00516       RETURN
00517 *
00518 *     End of DSYTRI2X
00519 *
00520       END
00521 
 All Files Functions