LAPACK 3.3.1 Linear Algebra PACKage

ssytri2x.f

Go to the documentation of this file.
```00001       SUBROUTINE SSYTRI2X( 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       REAL               A( LDA, * ), WORK( N+NB+1,* )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  SSYTRI2X 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 *  SSYTRF.
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) REAL 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 SSYTRF.
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 SSYTRF.
00055 *
00056 *  WORK    (workspace) REAL 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       REAL               ONE, ZERO
00071       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+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       REAL               AK, AKKP1, AKP1, D, T
00080       REAL               U01_I_J, U01_IP1_J
00081       REAL               U11_I_J, U11_IP1_J
00082 *     ..
00083 *     .. External Functions ..
00084       LOGICAL            LSAME
00085       EXTERNAL           LSAME
00086 *     ..
00087 *     .. External Subroutines ..
00088       EXTERNAL           SSYCONV, XERBLA, STRTRI
00089       EXTERNAL           SGEMM, STRMM, SSYSWAPR
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( 'SSYTRI2X', -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 SSYCONV( 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 STRTRI( 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 STRMM('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 SGEMM('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 *        U11 =  U11**T*invD1*U11 + U01**T*invD*U01
00286 *
00287          DO I=1,NNB
00288             DO J=I,NNB
00289               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00290             END DO
00291          END DO
00292 *
00293 *        U01 =  U00**T*invD0*U01
00294 *
00295          CALL STRMM('L',UPLO,'T','U',CUT, NNB,
00296      \$             ONE,A,LDA,WORK,N+NB+1)
00297
00298 *
00299 *        Update U01
00300 *
00301          DO I=1,CUT
00302            DO J=1,NNB
00303             A(I,CUT+J)=WORK(I,J)
00304            END DO
00305          END DO
00306 *
00307 *      Next Block
00308 *
00309        END DO
00310 *
00311 *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
00312 *
00313             I=1
00314             DO WHILE ( I .LE. N )
00315                IF( IPIV(I) .GT. 0 ) THEN
00316                   IP=IPIV(I)
00317                  IF (I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP )
00318                  IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I )
00319                ELSE
00320                  IP=-IPIV(I)
00321                  I=I+1
00322                  IF ( (I-1) .LT. IP)
00323      \$                  CALL SSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
00324                  IF ( (I-1) .GT. IP)
00325      \$                  CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I-1 )
00326               ENDIF
00327                I=I+1
00328             END DO
00329       ELSE
00330 *
00331 *        LOWER...
00332 *
00333 *        invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
00334 *
00335          CALL STRTRI( UPLO, 'U', N, A, LDA, INFO )
00336 *
00337 *       inv(D) and inv(D)*inv(U)
00338 *
00339         K=N
00340         DO WHILE ( K .GE. 1 )
00341          IF( IPIV( K ).GT.0 ) THEN
00342 *           1 x 1 diagonal NNB
00343              WORK(K,INVD) = ONE /  A( K, K )
00344              WORK(K,INVD+1) = 0
00345             K=K-1
00346          ELSE
00347 *           2 x 2 diagonal NNB
00348              T = WORK(K-1,1)
00349              AK = A( K-1, K-1 ) / T
00350              AKP1 = A( K, K ) / T
00351              AKKP1 = WORK(K-1,1) / T
00352              D = T*( AK*AKP1-ONE )
00353              WORK(K-1,INVD) = AKP1 / D
00354              WORK(K,INVD) = AK / D
00355              WORK(K,INVD+1) = -AKKP1 / D
00356              WORK(K-1,INVD+1) = -AKKP1 / D
00357             K=K-2
00358          END IF
00359         END DO
00360 *
00361 *       inv(U**T) = (inv(U))**T
00362 *
00363 *       inv(U**T)*inv(D)*inv(U)
00364 *
00365         CUT=0
00366         DO WHILE (CUT .LT. N)
00367            NNB=NB
00368            IF (CUT + NNB .GT. N) THEN
00369               NNB=N-CUT
00370            ELSE
00371               COUNT = 0
00372 *             count negative elements,
00373               DO I=CUT+1,CUT+NNB
00374                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
00375               END DO
00376 *             need a even number for a clear cut
00377               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
00378            END IF
00379 *     L21 Block
00380            DO I=1,N-CUT-NNB
00381              DO J=1,NNB
00382               WORK(I,J)=A(CUT+NNB+I,CUT+J)
00383              END DO
00384            END DO
00385 *     L11 Block
00386            DO I=1,NNB
00387              WORK(U11+I,I)=ONE
00388              DO J=I+1,NNB
00389                 WORK(U11+I,J)=ZERO
00390              END DO
00391              DO J=1,I-1
00392                 WORK(U11+I,J)=A(CUT+I,CUT+J)
00393              END DO
00394            END DO
00395 *
00396 *          invD*L21
00397 *
00398            I=N-CUT-NNB
00399            DO WHILE (I .GE. 1)
00400              IF (IPIV(CUT+NNB+I) > 0) THEN
00401                 DO J=1,NNB
00402                     WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
00403                 END DO
00404                 I=I-1
00405              ELSE
00406                 DO J=1,NNB
00407                    U01_I_J = WORK(I,J)
00408                    U01_IP1_J = WORK(I-1,J)
00409                    WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
00410      \$                        WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
00411                    WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
00412      \$                        WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
00413                 END DO
00414                 I=I-2
00415              END IF
00416            END DO
00417 *
00418 *        invD1*L11
00419 *
00420            I=NNB
00421            DO WHILE (I .GE. 1)
00422              IF (IPIV(CUT+I) > 0) THEN
00423                 DO J=1,NNB
00424                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
00425                 END DO
00426                 I=I-1
00427              ELSE
00428                 DO J=1,NNB
00429                    U11_I_J = WORK(U11+I,J)
00430                    U11_IP1_J = WORK(U11+I-1,J)
00431                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
00432      \$                      WORK(CUT+I,INVD+1)*U11_IP1_J
00433                 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
00434      \$                      WORK(CUT+I-1,INVD)*U11_IP1_J
00435                 END DO
00436                 I=I-2
00437              END IF
00438            END DO
00439 *
00440 *       L11**T*invD1*L11->L11
00441 *
00442         CALL STRMM('L',UPLO,'T','U',NNB, NNB,
00443      \$             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00444
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         IF ( (CUT+NNB) .LT. N ) THEN
00453 *
00454 *          L21**T*invD2*L21->A(CUT+I,CUT+J)
00455 *
00456          CALL SGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
00457      \$             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
00458
00459 *
00460 *        L11 =  L11**T*invD1*L11 + U01**T*invD*U01
00461 *
00462          DO I=1,NNB
00463             DO J=1,I
00464               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00465             END DO
00466          END DO
00467 *
00468 *        L01 =  L22**T*invD2*L21
00469 *
00470          CALL STRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
00471      \$             ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
00472 *
00473 *      Update L21
00474 *
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
00481        ELSE
00482 *
00483 *        L11 =  L11**T*invD1*L11
00484 *
00485          DO I=1,NNB
00486             DO J=1,I
00487               A(CUT+I,CUT+J)=WORK(U11+I,J)
00488             END DO
00489          END DO
00490        END IF
00491 *
00492 *      Next Block
00493 *
00494            CUT=CUT+NNB
00495        END DO
00496 *
00497 *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
00498 *
00499             I=N
00500             DO WHILE ( I .GE. 1 )
00501                IF( IPIV(I) .GT. 0 ) THEN
00502                   IP=IPIV(I)
00503                  IF (I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP  )
00504                  IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I )
00505                ELSE
00506                  IP=-IPIV(I)
00507                  IF ( I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP )
00508                  IF ( I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I )
00509                  I=I-1
00510                ENDIF
00511                I=I-1
00512             END DO
00513       END IF
00514 *
00515       RETURN
00516 *
00517 *     End of SSYTRI2X
00518 *
00519       END
00520
```