LAPACK: SRC/csptri.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CSPTRI( UPLO, N, AP, IPIV, WORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
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 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IPIV( * )
00014       COMPLEX            AP( * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  CSPTRI computes the inverse of a complex symmetric indefinite matrix
00021 *  A in packed storage using the factorization A = U*D*U**T or
00022 *  A = L*D*L**T computed by CSPTRF.
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  UPLO    (input) CHARACTER*1
00028 *          Specifies whether the details of the factorization are stored
00029 *          as an upper or lower triangular matrix.
00030 *          = 'U':  Upper triangular, form is A = U*D*U**T;
00031 *          = 'L':  Lower triangular, form is A = L*D*L**T.
00032 *
00033 *  N       (input) INTEGER
00034 *          The order of the matrix A.  N >= 0.
00035 *
00036 *  AP      (input/output) COMPLEX array, dimension (N*(N+1)/2)
00037 *          On entry, the block diagonal matrix D and the multipliers
00038 *          used to obtain the factor U or L as computed by CSPTRF,
00039 *          stored as a packed triangular matrix.
00040 *
00041 *          On exit, if INFO = 0, the (symmetric) inverse of the original
00042 *          matrix, stored as a packed triangular matrix. The j-th column
00043 *          of inv(A) is stored in the array AP as follows:
00044 *          if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
00045 *          if UPLO = 'L',
00046 *             AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
00047 *
00048 *  IPIV    (input) INTEGER array, dimension (N)
00049 *          Details of the interchanges and the block structure of D
00050 *          as determined by CSPTRF.
00051 *
00052 *  WORK    (workspace) COMPLEX array, dimension (N)
00053 *
00054 *  INFO    (output) INTEGER
00055 *          = 0: successful exit
00056 *          < 0: if INFO = -i, the i-th argument had an illegal value
00057 *          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
00058 *               inverse could not be computed.
00059 *
00060 *  =====================================================================
00061 *
00062 *     .. Parameters ..
00063       COMPLEX            ONE, ZERO
00064       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
00065      $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
00066 *     ..
00067 *     .. Local Scalars ..
00068       LOGICAL            UPPER
00069       INTEGER            J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
00070       COMPLEX            AK, AKKP1, AKP1, D, T, TEMP
00071 *     ..
00072 *     .. External Functions ..
00073       LOGICAL            LSAME
00074       COMPLEX            CDOTU
00075       EXTERNAL           LSAME, CDOTU
00076 *     ..
00077 *     .. External Subroutines ..
00078       EXTERNAL           CCOPY, CSPMV, CSWAP, XERBLA
00079 *     ..
00080 *     .. Intrinsic Functions ..
00081       INTRINSIC          ABS
00082 *     ..
00083 *     .. Executable Statements ..
00084 *
00085 *     Test the input parameters.
00086 *
00087       INFO = 0
00088       UPPER = LSAME( UPLO, 'U' )
00089       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00090          INFO = -1
00091       ELSE IF( N.LT.0 ) THEN
00092          INFO = -2
00093       END IF
00094       IF( INFO.NE.0 ) THEN
00095          CALL XERBLA( 'CSPTRI', -INFO )
00096          RETURN
00097       END IF
00098 *
00099 *     Quick return if possible
00100 *
00101       IF( N.EQ.0 )
00102      $   RETURN
00103 *
00104 *     Check that the diagonal matrix D is nonsingular.
00105 *
00106       IF( UPPER ) THEN
00107 *
00108 *        Upper triangular storage: examine D from bottom to top
00109 *
00110          KP = N*( N+1 ) / 2
00111          DO 10 INFO = N, 1, -1
00112             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
00113      $         RETURN
00114             KP = KP - INFO
00115    10    CONTINUE
00116       ELSE
00117 *
00118 *        Lower triangular storage: examine D from top to bottom.
00119 *
00120          KP = 1
00121          DO 20 INFO = 1, N
00122             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
00123      $         RETURN
00124             KP = KP + N - INFO + 1
00125    20    CONTINUE
00126       END IF
00127       INFO = 0
00128 *
00129       IF( UPPER ) THEN
00130 *
00131 *        Compute inv(A) from the factorization A = U*D*U'.
00132 *
00133 *        K is the main loop index, increasing from 1 to N in steps of
00134 *        1 or 2, depending on the size of the diagonal blocks.
00135 *
00136          K = 1
00137          KC = 1
00138    30    CONTINUE
00139 *
00140 *        If K > N, exit from loop.
00141 *
00142          IF( K.GT.N )
00143      $      GO TO 50
00144 *
00145          KCNEXT = KC + K
00146          IF( IPIV( K ).GT.0 ) THEN
00147 *
00148 *           1 x 1 diagonal block
00149 *
00150 *           Invert the diagonal block.
00151 *
00152             AP( KC+K-1 ) = ONE / AP( KC+K-1 )
00153 *
00154 *           Compute column K of the inverse.
00155 *
00156             IF( K.GT.1 ) THEN
00157                CALL CCOPY( K-1, AP( KC ), 1, WORK, 1 )
00158                CALL CSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
00159      $                     1 )
00160                AP( KC+K-1 ) = AP( KC+K-1 ) -
00161      $                        CDOTU( K-1, WORK, 1, AP( KC ), 1 )
00162             END IF
00163             KSTEP = 1
00164          ELSE
00165 *
00166 *           2 x 2 diagonal block
00167 *
00168 *           Invert the diagonal block.
00169 *
00170             T = AP( KCNEXT+K-1 )
00171             AK = AP( KC+K-1 ) / T
00172             AKP1 = AP( KCNEXT+K ) / T
00173             AKKP1 = AP( KCNEXT+K-1 ) / T
00174             D = T*( AK*AKP1-ONE )
00175             AP( KC+K-1 ) = AKP1 / D
00176             AP( KCNEXT+K ) = AK / D
00177             AP( KCNEXT+K-1 ) = -AKKP1 / D
00178 *
00179 *           Compute columns K and K+1 of the inverse.
00180 *
00181             IF( K.GT.1 ) THEN
00182                CALL CCOPY( K-1, AP( KC ), 1, WORK, 1 )
00183                CALL CSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
00184      $                     1 )
00185                AP( KC+K-1 ) = AP( KC+K-1 ) -
00186      $                        CDOTU( K-1, WORK, 1, AP( KC ), 1 )
00187                AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
00188      $                            CDOTU( K-1, AP( KC ), 1, AP( KCNEXT ),
00189      $                            1 )
00190                CALL CCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
00191                CALL CSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
00192      $                     AP( KCNEXT ), 1 )
00193                AP( KCNEXT+K ) = AP( KCNEXT+K ) -
00194      $                          CDOTU( K-1, WORK, 1, AP( KCNEXT ), 1 )
00195             END IF
00196             KSTEP = 2
00197             KCNEXT = KCNEXT + K + 1
00198          END IF
00199 *
00200          KP = ABS( IPIV( K ) )
00201          IF( KP.NE.K ) THEN
00202 *
00203 *           Interchange rows and columns K and KP in the leading
00204 *           submatrix A(1:k+1,1:k+1)
00205 *
00206             KPC = ( KP-1 )*KP / 2 + 1
00207             CALL CSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
00208             KX = KPC + KP - 1
00209             DO 40 J = KP + 1, K - 1
00210                KX = KX + J - 1
00211                TEMP = AP( KC+J-1 )
00212                AP( KC+J-1 ) = AP( KX )
00213                AP( KX ) = TEMP
00214    40       CONTINUE
00215             TEMP = AP( KC+K-1 )
00216             AP( KC+K-1 ) = AP( KPC+KP-1 )
00217             AP( KPC+KP-1 ) = TEMP
00218             IF( KSTEP.EQ.2 ) THEN
00219                TEMP = AP( KC+K+K-1 )
00220                AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
00221                AP( KC+K+KP-1 ) = TEMP
00222             END IF
00223          END IF
00224 *
00225          K = K + KSTEP
00226          KC = KCNEXT
00227          GO TO 30
00228    50    CONTINUE
00229 *
00230       ELSE
00231 *
00232 *        Compute inv(A) from the factorization A = L*D*L'.
00233 *
00234 *        K is the main loop index, increasing from 1 to N in steps of
00235 *        1 or 2, depending on the size of the diagonal blocks.
00236 *
00237          NPP = N*( N+1 ) / 2
00238          K = N
00239          KC = NPP
00240    60    CONTINUE
00241 *
00242 *        If K < 1, exit from loop.
00243 *
00244          IF( K.LT.1 )
00245      $      GO TO 80
00246 *
00247          KCNEXT = KC - ( N-K+2 )
00248          IF( IPIV( K ).GT.0 ) THEN
00249 *
00250 *           1 x 1 diagonal block
00251 *
00252 *           Invert the diagonal block.
00253 *
00254             AP( KC ) = ONE / AP( KC )
00255 *
00256 *           Compute column K of the inverse.
00257 *
00258             IF( K.LT.N ) THEN
00259                CALL CCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
00260                CALL CSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
00261      $                     ZERO, AP( KC+1 ), 1 )
00262                AP( KC ) = AP( KC ) - CDOTU( N-K, WORK, 1, AP( KC+1 ),
00263      $                    1 )
00264             END IF
00265             KSTEP = 1
00266          ELSE
00267 *
00268 *           2 x 2 diagonal block
00269 *
00270 *           Invert the diagonal block.
00271 *
00272             T = AP( KCNEXT+1 )
00273             AK = AP( KCNEXT ) / T
00274             AKP1 = AP( KC ) / T
00275             AKKP1 = AP( KCNEXT+1 ) / T
00276             D = T*( AK*AKP1-ONE )
00277             AP( KCNEXT ) = AKP1 / D
00278             AP( KC ) = AK / D
00279             AP( KCNEXT+1 ) = -AKKP1 / D
00280 *
00281 *           Compute columns K-1 and K of the inverse.
00282 *
00283             IF( K.LT.N ) THEN
00284                CALL CCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
00285                CALL CSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
00286      $                     ZERO, AP( KC+1 ), 1 )
00287                AP( KC ) = AP( KC ) - CDOTU( N-K, WORK, 1, AP( KC+1 ),
00288      $                    1 )
00289                AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
00290      $                          CDOTU( N-K, AP( KC+1 ), 1,
00291      $                          AP( KCNEXT+2 ), 1 )
00292                CALL CCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
00293                CALL CSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
00294      $                     ZERO, AP( KCNEXT+2 ), 1 )
00295                AP( KCNEXT ) = AP( KCNEXT ) -
00296      $                        CDOTU( N-K, WORK, 1, AP( KCNEXT+2 ), 1 )
00297             END IF
00298             KSTEP = 2
00299             KCNEXT = KCNEXT - ( N-K+3 )
00300          END IF
00301 *
00302          KP = ABS( IPIV( K ) )
00303          IF( KP.NE.K ) THEN
00304 *
00305 *           Interchange rows and columns K and KP in the trailing
00306 *           submatrix A(k-1:n,k-1:n)
00307 *
00308             KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
00309             IF( KP.LT.N )
00310      $         CALL CSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
00311             KX = KC + KP - K
00312             DO 70 J = K + 1, KP - 1
00313                KX = KX + N - J + 1
00314                TEMP = AP( KC+J-K )
00315                AP( KC+J-K ) = AP( KX )
00316                AP( KX ) = TEMP
00317    70       CONTINUE
00318             TEMP = AP( KC )
00319             AP( KC ) = AP( KPC )
00320             AP( KPC ) = TEMP
00321             IF( KSTEP.EQ.2 ) THEN
00322                TEMP = AP( KC-N+K-1 )
00323                AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
00324                AP( KC-N+KP-1 ) = TEMP
00325             END IF
00326          END IF
00327 *
00328          K = K - KSTEP
00329          KC = KCNEXT
00330          GO TO 60
00331    80    CONTINUE
00332       END IF
00333 *
00334       RETURN
00335 *
00336 *     End of CSPTRI
00337 *
00338       END