LAPACK: SRC/zhptri.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZHPTRI( 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*16         AP( * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZHPTRI 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 ZHPTRF.
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*16 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 ZHPTRF,
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 ZHPTRF.
00051 *
00052 *  WORK    (workspace) COMPLEX*16 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       DOUBLE PRECISION   ONE
00064       COMPLEX*16         CONE, ZERO
00065       PARAMETER          ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
00066      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
00067 *     ..
00068 *     .. Local Scalars ..
00069       LOGICAL            UPPER
00070       INTEGER            J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
00071       DOUBLE PRECISION   AK, AKP1, D, T
00072       COMPLEX*16         AKKP1, TEMP
00073 *     ..
00074 *     .. External Functions ..
00075       LOGICAL            LSAME
00076       COMPLEX*16         ZDOTC
00077       EXTERNAL           LSAME, ZDOTC
00078 *     ..
00079 *     .. External Subroutines ..
00080       EXTERNAL           XERBLA, ZCOPY, ZHPMV, ZSWAP
00081 *     ..
00082 *     .. Intrinsic Functions ..
00083       INTRINSIC          ABS, DBLE, DCONJG
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( 'ZHPTRI', -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 / DBLE( AP( KC+K-1 ) )
00155 *
00156 *           Compute column K of the inverse.
00157 *
00158             IF( K.GT.1 ) THEN
00159                CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
00160                CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
00161      $                     AP( KC ), 1 )
00162                AP( KC+K-1 ) = AP( KC+K-1 ) -
00163      $                        DBLE( ZDOTC( 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 = DBLE( AP( KC+K-1 ) ) / T
00174             AKP1 = DBLE( 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 ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
00185                CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
00186      $                     AP( KC ), 1 )
00187                AP( KC+K-1 ) = AP( KC+K-1 ) -
00188      $                        DBLE( ZDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
00189                AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
00190      $                            ZDOTC( K-1, AP( KC ), 1, AP( KCNEXT ),
00191      $                            1 )
00192                CALL ZCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
00193                CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
00194      $                     AP( KCNEXT ), 1 )
00195                AP( KCNEXT+K ) = AP( KCNEXT+K ) -
00196      $                          DBLE( ZDOTC( 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 ZSWAP( 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 = DCONJG( AP( KC+J-1 ) )
00215                AP( KC+J-1 ) = DCONJG( AP( KX ) )
00216                AP( KX ) = TEMP
00217    40       CONTINUE
00218             AP( KC+KP-1 ) = DCONJG( 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 / DBLE( AP( KC ) )
00259 *
00260 *           Compute column K of the inverse.
00261 *
00262             IF( K.LT.N ) THEN
00263                CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
00264                CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+N-K+1 ), WORK, 1,
00265      $                     ZERO, AP( KC+1 ), 1 )
00266                AP( KC ) = AP( KC ) - DBLE( ZDOTC( N-K, WORK, 1,
00267      $                    AP( KC+1 ), 1 ) )
00268             END IF
00269             KSTEP = 1
00270          ELSE
00271 *
00272 *           2 x 2 diagonal block
00273 *
00274 *           Invert the diagonal block.
00275 *
00276             T = ABS( AP( KCNEXT+1 ) )
00277             AK = DBLE( AP( KCNEXT ) ) / T
00278             AKP1 = DBLE( AP( KC ) ) / T
00279             AKKP1 = AP( KCNEXT+1 ) / T
00280             D = T*( AK*AKP1-ONE )
00281             AP( KCNEXT ) = AKP1 / D
00282             AP( KC ) = AK / D
00283             AP( KCNEXT+1 ) = -AKKP1 / D
00284 *
00285 *           Compute columns K-1 and K of the inverse.
00286 *
00287             IF( K.LT.N ) THEN
00288                CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
00289                CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
00290      $                     1, ZERO, AP( KC+1 ), 1 )
00291                AP( KC ) = AP( KC ) - DBLE( ZDOTC( N-K, WORK, 1,
00292      $                    AP( KC+1 ), 1 ) )
00293                AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
00294      $                          ZDOTC( N-K, AP( KC+1 ), 1,
00295      $                          AP( KCNEXT+2 ), 1 )
00296                CALL ZCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
00297                CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
00298      $                     1, ZERO, AP( KCNEXT+2 ), 1 )
00299                AP( KCNEXT ) = AP( KCNEXT ) -
00300      $                        DBLE( ZDOTC( N-K, WORK, 1, AP( KCNEXT+2 ),
00301      $                        1 ) )
00302             END IF
00303             KSTEP = 2
00304             KCNEXT = KCNEXT - ( N-K+3 )
00305          END IF
00306 *
00307          KP = ABS( IPIV( K ) )
00308          IF( KP.NE.K ) THEN
00309 *
00310 *           Interchange rows and columns K and KP in the trailing
00311 *           submatrix A(k-1:n,k-1:n)
00312 *
00313             KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
00314             IF( KP.LT.N )
00315      $         CALL ZSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
00316             KX = KC + KP - K
00317             DO 70 J = K + 1, KP - 1
00318                KX = KX + N - J + 1
00319                TEMP = DCONJG( AP( KC+J-K ) )
00320                AP( KC+J-K ) = DCONJG( AP( KX ) )
00321                AP( KX ) = TEMP
00322    70       CONTINUE
00323             AP( KC+KP-K ) = DCONJG( AP( KC+KP-K ) )
00324             TEMP = AP( KC )
00325             AP( KC ) = AP( KPC )
00326             AP( KPC ) = TEMP
00327             IF( KSTEP.EQ.2 ) THEN
00328                TEMP = AP( KC-N+K-1 )
00329                AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
00330                AP( KC-N+KP-1 ) = TEMP
00331             END IF
00332          END IF
00333 *
00334          K = K - KSTEP
00335          KC = KCNEXT
00336          GO TO 60
00337    80    CONTINUE
00338       END IF
00339 *
00340       RETURN
00341 *
00342 *     End of ZHPTRI
00343 *
00344       END