LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ssytri2x.f
Go to the documentation of this file.
00001 *> \brief \b SSYTRI2X
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SSYTRI2X + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytri2x.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytri2x.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytri2x.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SSYTRI2X( 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 *       REAL               A( LDA, * ), WORK( N+NB+1,* )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> SSYTRI2X 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 *> SSYTRF.
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 REAL 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 SSYTRF.
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 SSYTRF.
00086 *> \endverbatim
00087 *>
00088 *> \param[out] WORK
00089 *> \verbatim
00090 *>          WORK is REAL 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 realSYcomputational
00119 *
00120 *  =====================================================================
00121       SUBROUTINE SSYTRI2X( 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       REAL               A( LDA, * ), WORK( N+NB+1,* )
00135 *     ..
00136 *
00137 *  =====================================================================
00138 *
00139 *     .. Parameters ..
00140       REAL               ONE, ZERO
00141       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+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       REAL               AK, AKKP1, AKP1, D, T
00150       REAL               U01_I_J, U01_IP1_J
00151       REAL               U11_I_J, U11_IP1_J
00152 *     ..
00153 *     .. External Functions ..
00154       LOGICAL            LSAME
00155       EXTERNAL           LSAME
00156 *     ..
00157 *     .. External Subroutines ..
00158       EXTERNAL           SSYCONV, XERBLA, STRTRI
00159       EXTERNAL           SGEMM, STRMM, SSYSWAPR
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( 'SSYTRI2X', -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 SSYCONV( 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 STRTRI( 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 STRMM('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 SGEMM('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 *        U11 =  U11**T*invD1*U11 + U01**T*invD*U01
00356 *
00357          DO I=1,NNB
00358             DO J=I,NNB
00359               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00360             END DO
00361          END DO
00362 *
00363 *        U01 =  U00**T*invD0*U01
00364 *
00365          CALL STRMM('L',UPLO,'T','U',CUT, NNB,
00366      $             ONE,A,LDA,WORK,N+NB+1)
00367 
00368 *
00369 *        Update U01
00370 *
00371          DO I=1,CUT
00372            DO J=1,NNB
00373             A(I,CUT+J)=WORK(I,J)
00374            END DO
00375          END DO
00376 *
00377 *      Next Block
00378 *
00379        END DO
00380 *
00381 *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
00382 *  
00383             I=1
00384             DO WHILE ( I .LE. N )
00385                IF( IPIV(I) .GT. 0 ) THEN
00386                   IP=IPIV(I)
00387                  IF (I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP )
00388                  IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I )
00389                ELSE
00390                  IP=-IPIV(I)
00391                  I=I+1
00392                  IF ( (I-1) .LT. IP) 
00393      $                  CALL SSYSWAPR( UPLO, N, A, LDA, I-1 ,IP )
00394                  IF ( (I-1) .GT. IP) 
00395      $                  CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I-1 )
00396               ENDIF
00397                I=I+1
00398             END DO
00399       ELSE
00400 *
00401 *        LOWER...
00402 *
00403 *        invA = P * inv(U**T)*inv(D)*inv(U)*P**T.
00404 *
00405          CALL STRTRI( UPLO, 'U', N, A, LDA, INFO )
00406 *
00407 *       inv(D) and inv(D)*inv(U)
00408 * 
00409         K=N
00410         DO WHILE ( K .GE. 1 )
00411          IF( IPIV( K ).GT.0 ) THEN
00412 *           1 x 1 diagonal NNB
00413              WORK(K,INVD) = ONE /  A( K, K )
00414              WORK(K,INVD+1) = 0
00415             K=K-1
00416          ELSE
00417 *           2 x 2 diagonal NNB
00418              T = WORK(K-1,1)
00419              AK = A( K-1, K-1 ) / T
00420              AKP1 = A( K, K ) / T
00421              AKKP1 = WORK(K-1,1) / T
00422              D = T*( AK*AKP1-ONE )
00423              WORK(K-1,INVD) = AKP1 / D
00424              WORK(K,INVD) = AK / D
00425              WORK(K,INVD+1) = -AKKP1 / D  
00426              WORK(K-1,INVD+1) = -AKKP1 / D  
00427             K=K-2
00428          END IF
00429         END DO
00430 *
00431 *       inv(U**T) = (inv(U))**T
00432 *
00433 *       inv(U**T)*inv(D)*inv(U)
00434 *
00435         CUT=0
00436         DO WHILE (CUT .LT. N)
00437            NNB=NB
00438            IF (CUT + NNB .GT. N) THEN
00439               NNB=N-CUT
00440            ELSE
00441               COUNT = 0
00442 *             count negative elements, 
00443               DO I=CUT+1,CUT+NNB
00444                   IF (IPIV(I) .LT. 0) COUNT=COUNT+1
00445               END DO
00446 *             need a even number for a clear cut
00447               IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
00448            END IF
00449 *     L21 Block
00450            DO I=1,N-CUT-NNB
00451              DO J=1,NNB
00452               WORK(I,J)=A(CUT+NNB+I,CUT+J)
00453              END DO
00454            END DO
00455 *     L11 Block
00456            DO I=1,NNB
00457              WORK(U11+I,I)=ONE
00458              DO J=I+1,NNB
00459                 WORK(U11+I,J)=ZERO
00460              END DO
00461              DO J=1,I-1
00462                 WORK(U11+I,J)=A(CUT+I,CUT+J)
00463              END DO
00464            END DO
00465 *
00466 *          invD*L21
00467 *
00468            I=N-CUT-NNB
00469            DO WHILE (I .GE. 1)
00470              IF (IPIV(CUT+NNB+I) > 0) THEN
00471                 DO J=1,NNB
00472                     WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
00473                 END DO
00474                 I=I-1
00475              ELSE
00476                 DO J=1,NNB
00477                    U01_I_J = WORK(I,J)
00478                    U01_IP1_J = WORK(I-1,J)
00479                    WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
00480      $                        WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
00481                    WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
00482      $                        WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
00483                 END DO
00484                 I=I-2
00485              END IF
00486            END DO
00487 *
00488 *        invD1*L11
00489 *
00490            I=NNB
00491            DO WHILE (I .GE. 1)
00492              IF (IPIV(CUT+I) > 0) THEN
00493                 DO J=1,NNB
00494                     WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
00495                 END DO
00496                 I=I-1
00497              ELSE
00498                 DO J=1,NNB
00499                    U11_I_J = WORK(U11+I,J)
00500                    U11_IP1_J = WORK(U11+I-1,J)
00501                 WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
00502      $                      WORK(CUT+I,INVD+1)*U11_IP1_J
00503                 WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
00504      $                      WORK(CUT+I-1,INVD)*U11_IP1_J
00505                 END DO
00506                 I=I-2
00507              END IF
00508            END DO
00509 *    
00510 *       L11**T*invD1*L11->L11
00511 *
00512         CALL STRMM('L',UPLO,'T','U',NNB, NNB,
00513      $             ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
00514 
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         IF ( (CUT+NNB) .LT. N ) THEN
00523 *
00524 *          L21**T*invD2*L21->A(CUT+I,CUT+J)
00525 *
00526          CALL SGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1)
00527      $             ,LDA,WORK,N+NB+1, ZERO, WORK(U11+1,1), N+NB+1)
00528        
00529 *
00530 *        L11 =  L11**T*invD1*L11 + U01**T*invD*U01
00531 *
00532          DO I=1,NNB
00533             DO J=1,I
00534               A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
00535             END DO
00536          END DO
00537 *
00538 *        L01 =  L22**T*invD2*L21
00539 *
00540          CALL STRMM('L',UPLO,'T','U', N-NNB-CUT, NNB,
00541      $             ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
00542 *
00543 *      Update L21
00544 *
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 
00551        ELSE
00552 *
00553 *        L11 =  L11**T*invD1*L11
00554 *
00555          DO I=1,NNB
00556             DO J=1,I
00557               A(CUT+I,CUT+J)=WORK(U11+I,J)
00558             END DO
00559          END DO
00560        END IF
00561 *
00562 *      Next Block
00563 *
00564            CUT=CUT+NNB
00565        END DO
00566 *
00567 *        Apply PERMUTATIONS P and P**T: P * inv(U**T)*inv(D)*inv(U) *P**T
00568 * 
00569             I=N
00570             DO WHILE ( I .GE. 1 )
00571                IF( IPIV(I) .GT. 0 ) THEN
00572                   IP=IPIV(I)
00573                  IF (I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP  )
00574                  IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I )
00575                ELSE
00576                  IP=-IPIV(I)
00577                  IF ( I .LT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, I ,IP )
00578                  IF ( I .GT. IP) CALL SSYSWAPR( UPLO, N, A, LDA, IP ,I )
00579                  I=I-1
00580                ENDIF
00581                I=I-1
00582             END DO
00583       END IF
00584 *
00585       RETURN
00586 *
00587 *     End of SSYTRI2X
00588 *
00589       END
00590 
 All Files Functions