LAPACK: SRC/chptri.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CHPTRI( 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 *  CHPTRI computes the inverse of a complex Hermitian indefinite matrix
00021 *  A in packed storage using the factorization A = U*D*U**H or
00022 *  A = L*D*L**H computed by CHPTRF.
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**H;
00031 *          = 'L':  Lower triangular, form is A = L*D*L**H.
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 CHPTRF,
00039 *          stored as a packed triangular matrix.
00040 *
00041 *          On exit, if INFO = 0, the (Hermitian) 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 CHPTRF.
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       REAL               ONE
00064       COMPLEX            CONE, ZERO
00065       PARAMETER          ( ONE = 1.0E+0, CONE = ( 1.0E+0, 0.0E+0 ),
00066      $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
00067 *     ..
00068 *     .. Local Scalars ..
00069       LOGICAL            UPPER
00070       INTEGER            J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
00071       REAL               AK, AKP1, D, T
00072       COMPLEX            AKKP1, TEMP
00073 *     ..
00074 *     .. External Functions ..
00075       LOGICAL            LSAME
00076       COMPLEX            CDOTC
00077       EXTERNAL           LSAME, CDOTC
00078 *     ..
00079 *     .. External Subroutines ..
00080       EXTERNAL           CCOPY, CHPMV, CSWAP, XERBLA
00081 *     ..
00082 *     .. Intrinsic Functions ..
00083       INTRINSIC          ABS, CONJG, REAL
00084 *     ..
00085 *     .. Executable Statements ..
00086 *
00087 *     Test the input parameters.
00088 *
00089       INFO = 0
00090       UPPER = LSAME( UPLO, 'U' )
00091       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00092          INFO = -1
00093       ELSE IF( N.LT.0 ) THEN
00094          INFO = -2
00095       END IF
00096       IF( INFO.NE.0 ) THEN
00097          CALL XERBLA( 'CHPTRI', -INFO )
00098          RETURN
00099       END IF
00100 *
00101 *     Quick return if possible
00102 *
00103       IF( N.EQ.0 )
00104      $   RETURN
00105 *
00106 *     Check that the diagonal matrix D is nonsingular.
00107 *
00108       IF( UPPER ) THEN
00109 *
00110 *        Upper triangular storage: examine D from bottom to top
00111 *
00112          KP = N*( N+1 ) / 2
00113          DO 10 INFO = N, 1, -1
00114             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
00115      $         RETURN
00116             KP = KP - INFO
00117    10    CONTINUE
00118       ELSE
00119 *
00120 *        Lower triangular storage: examine D from top to bottom.
00121 *
00122          KP = 1
00123          DO 20 INFO = 1, N
00124             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
00125      $         RETURN
00126             KP = KP + N - INFO + 1
00127    20    CONTINUE
00128       END IF
00129       INFO = 0
00130 *
00131       IF( UPPER ) THEN
00132 *
00133 *        Compute inv(A) from the factorization A = U*D*U'.
00134 *
00135 *        K is the main loop index, increasing from 1 to N in steps of
00136 *        1 or 2, depending on the size of the diagonal blocks.
00137 *
00138          K = 1
00139          KC = 1
00140    30    CONTINUE
00141 *
00142 *        If K > N, exit from loop.
00143 *
00144          IF( K.GT.N )
00145      $      GO TO 50
00146 *
00147          KCNEXT = KC + K
00148          IF( IPIV( K ).GT.0 ) THEN
00149 *
00150 *           1 x 1 diagonal block
00151 *
00152 *           Invert the diagonal block.
00153 *
00154             AP( KC+K-1 ) = ONE / REAL( AP( KC+K-1 ) )
00155 *
00156 *           Compute column K of the inverse.
00157 *
00158             IF( K.GT.1 ) THEN
00159                CALL CCOPY( K-1, AP( KC ), 1, WORK, 1 )
00160                CALL CHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
00161      $                     AP( KC ), 1 )
00162                AP( KC+K-1 ) = AP( KC+K-1 ) -
00163      $                        REAL( CDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
00164             END IF
00165             KSTEP = 1
00166          ELSE
00167 *
00168 *           2 x 2 diagonal block
00169 *
00170 *           Invert the diagonal block.
00171 *
00172             T = ABS( AP( KCNEXT+K-1 ) )
00173             AK = REAL( AP( KC+K-1 ) ) / T
00174             AKP1 = REAL( AP( KCNEXT+K ) ) / T
00175             AKKP1 = AP( KCNEXT+K-1 ) / T
00176             D = T*( AK*AKP1-ONE )
00177             AP( KC+K-1 ) = AKP1 / D
00178             AP( KCNEXT+K ) = AK / D
00179             AP( KCNEXT+K-1 ) = -AKKP1 / D
00180 *
00181 *           Compute columns K and K+1 of the inverse.
00182 *
00183             IF( K.GT.1 ) THEN
00184                CALL CCOPY( K-1, AP( KC ), 1, WORK, 1 )
00185                CALL CHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
00186      $                     AP( KC ), 1 )
00187                AP( KC+K-1 ) = AP( KC+K-1 ) -
00188      $                        REAL( CDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
00189                AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
00190      $                            CDOTC( K-1, AP( KC ), 1, AP( KCNEXT ),
00191      $                            1 )
00192                CALL CCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
00193                CALL CHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
00194      $                     AP( KCNEXT ), 1 )
00195                AP( KCNEXT+K ) = AP( KCNEXT+K ) -
00196      $                          REAL( CDOTC( K-1, WORK, 1, AP( KCNEXT ),
00197      $                          1 ) )
00198             END IF
00199             KSTEP = 2
00200             KCNEXT = KCNEXT + K + 1
00201          END IF
00202 *
00203          KP = ABS( IPIV( K ) )
00204          IF( KP.NE.K ) THEN
00205 *
00206 *           Interchange rows and columns K and KP in the leading
00207 *           submatrix A(1:k+1,1:k+1)
00208 *
00209             KPC = ( KP-1 )*KP / 2 + 1
00210             CALL CSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
00211             KX = KPC + KP - 1
00212             DO 40 J = KP + 1, K - 1
00213                KX = KX + J - 1
00214                TEMP = CONJG( AP( KC+J-1 ) )
00215                AP( KC+J-1 ) = CONJG( AP( KX ) )
00216                AP( KX ) = TEMP
00217    40       CONTINUE
00218             AP( KC+KP-1 ) = CONJG( AP( KC+KP-1 ) )
00219             TEMP = AP( KC+K-1 )
00220             AP( KC+K-1 ) = AP( KPC+KP-1 )
00221             AP( KPC+KP-1 ) = TEMP
00222             IF( KSTEP.EQ.2 ) THEN
00223                TEMP = AP( KC+K+K-1 )
00224                AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
00225                AP( KC+K+KP-1 ) = TEMP
00226             END IF
00227          END IF
00228 *
00229          K = K + KSTEP
00230          KC = KCNEXT
00231          GO TO 30
00232    50    CONTINUE
00233 *
00234       ELSE
00235 *
00236 *        Compute inv(A) from the factorization A = L*D*L'.
00237 *
00238 *        K is the main loop index, increasing from 1 to N in steps of
00239 *        1 or 2, depending on the size of the diagonal blocks.
00240 *
00241          NPP = N*( N+1 ) / 2
00242          K = N
00243          KC = NPP
00244    60    CONTINUE
00245 *
00246 *        If K < 1, exit from loop.
00247 *
00248          IF( K.LT.1 )
00249      $      GO TO 80
00250 *
00251          KCNEXT = KC - ( N-K+2 )
00252          IF( IPIV( K ).GT.0 ) THEN
00253 *
00254 *           1 x 1 diagonal block
00255 *
00256 *           Invert the diagonal block.
00257 *
00258             AP( KC ) = ONE / REAL( AP( KC ) )
00259 *
00260 *           Compute column K of the inverse.
00261 *
00262             IF( K.LT.N ) THEN
00263                CALL CCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
00264                CALL CHPMV( UPLO, N-K, -CONE, AP( KC+N-K+1 ), WORK, 1,
00265      $                     ZERO, AP( KC+1 ), 1 )
00266                AP( KC ) = AP( KC ) - REAL( CDOTC( N-K, WORK, 1,
     $                    AP( KC+1 ), 1 ) )
00267             END IF
00268             KSTEP = 1
00269          ELSE
00270 *
00271 *           2 x 2 diagonal block
00272 *
00273 *           Invert the diagonal block.
00274 *
00275             T = ABS( AP( KCNEXT+1 ) )
00276             AK = REAL( AP( KCNEXT ) ) / T
00277             AKP1 = REAL( AP( KC ) ) / T
00278             AKKP1 = AP( KCNEXT+1 ) / T
00279             D = T*( AK*AKP1-ONE )
00280             AP( KCNEXT ) = AKP1 / D
00281             AP( KC ) = AK / D
00282             AP( KCNEXT+1 ) = -AKKP1 / D
00283 *
00284 *           Compute columns K-1 and K of the inverse.
00285 *
00286             IF( K.LT.N ) THEN
00287                CALL CCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
00288                CALL CHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
00289      $                     1, ZERO, AP( KC+1 ), 1 )
00290                AP( KC ) = AP( KC ) - REAL( CDOTC( N-K, WORK, 1,
     $                    AP( KC+1 ), 1 ) )
00291                AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
00292      $                          CDOTC( N-K, AP( KC+1 ), 1,
00293      $                          AP( KCNEXT+2 ), 1 )
00294                CALL CCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
00295                CALL CHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
00296      $                     1, ZERO, AP( KCNEXT+2 ), 1 )
00297                AP( KCNEXT ) = AP( KCNEXT ) -
00298      $                        REAL( CDOTC( N-K, WORK, 1, AP( KCNEXT+2 ),
00299      $                        1 ) )
00300             END IF
00301             KSTEP = 2
00302             KCNEXT = KCNEXT - ( N-K+3 )
00303          END IF
00304 *
00305          KP = ABS( IPIV( K ) )
00306          IF( KP.NE.K ) THEN
00307 *
00308 *           Interchange rows and columns K and KP in the trailing
00309 *           submatrix A(k-1:n,k-1:n)
00310 *
00311             KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
00312             IF( KP.LT.N )
00313      $         CALL CSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
00314             KX = KC + KP - K
00315             DO 70 J = K + 1, KP - 1
00316                KX = KX + N - J + 1
00317                TEMP = CONJG( AP( KC+J-K ) )
00318                AP( KC+J-K ) = CONJG( AP( KX ) )
00319                AP( KX ) = TEMP
00320    70       CONTINUE
00321             AP( KC+KP-K ) = CONJG( AP( KC+KP-K ) )
00322             TEMP = AP( KC )
00323             AP( KC ) = AP( KPC )
00324             AP( KPC ) = TEMP
00325             IF( KSTEP.EQ.2 ) THEN
00326                TEMP = AP( KC-N+K-1 )
00327                AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
00328                AP( KC-N+KP-1 ) = TEMP
00329             END IF
00330          END IF
00331 *
00332          K = K - KSTEP
00333          KC = KCNEXT
00334          GO TO 60
00335    80    CONTINUE
00336       END IF
00337 *
00338       RETURN
00339 *
00340 *     End of CHPTRI
00341 *
00342       END
00343