SRC/zhptri.f

   1 *> \brief \b ZHPTRI
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHPTRI + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhptri.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhptri.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhptri.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHPTRI( UPLO, N, AP, IPIV, WORK, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       INTEGER            IPIV( * )
  29 *       COMPLEX*16         AP( * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> ZHPTRI computes the inverse of a complex Hermitian indefinite matrix
  39 *> A in packed storage using the factorization A = U*D*U**H or
  40 *> A = L*D*L**H computed by ZHPTRF.
  41 *> \endverbatim
  42 *
  43 *  Arguments:
  44 *  ==========
  45 *
  46 *> \param[in] UPLO
  47 *> \verbatim
  48 *>          UPLO is CHARACTER*1
  49 *>          Specifies whether the details of the factorization are stored
  50 *>          as an upper or lower triangular matrix.
  51 *>          = 'U':  Upper triangular, form is A = U*D*U**H;
  52 *>          = 'L':  Lower triangular, form is A = L*D*L**H.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] N
  56 *> \verbatim
  57 *>          N is INTEGER
  58 *>          The order of the matrix A.  N >= 0.
  59 *> \endverbatim
  60 *>
  61 *> \param[in,out] AP
  62 *> \verbatim
  63 *>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
  64 *>          On entry, the block diagonal matrix D and the multipliers
  65 *>          used to obtain the factor U or L as computed by ZHPTRF,
  66 *>          stored as a packed triangular matrix.
  67 *>
  68 *>          On exit, if INFO = 0, the (Hermitian) inverse of the original
  69 *>          matrix, stored as a packed triangular matrix. The j-th column
  70 *>          of inv(A) is stored in the array AP as follows:
  71 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
  72 *>          if UPLO = 'L',
  73 *>             AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] IPIV
  77 *> \verbatim
  78 *>          IPIV is INTEGER array, dimension (N)
  79 *>          Details of the interchanges and the block structure of D
  80 *>          as determined by ZHPTRF.
  81 *> \endverbatim
  82 *>
  83 *> \param[out] WORK
  84 *> \verbatim
  85 *>          WORK is COMPLEX*16 array, dimension (N)
  86 *> \endverbatim
  87 *>
  88 *> \param[out] INFO
  89 *> \verbatim
  90 *>          INFO is INTEGER
  91 *>          = 0: successful exit
  92 *>          < 0: if INFO = -i, the i-th argument had an illegal value
  93 *>          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
  94 *>               inverse could not be computed.
  95 *> \endverbatim
  96 *
  97 *  Authors:
  98 *  ========
  99 *
 100 *> \author Univ. of Tennessee
 101 *> \author Univ. of California Berkeley
 102 *> \author Univ. of Colorado Denver
 103 *> \author NAG Ltd.
 104 *
 105 *> \date November 2011
 106 *
 107 *> \ingroup complex16OTHERcomputational
 108 *
 109 *  =====================================================================
 110       SUBROUTINE ZHPTRI( UPLO, N, AP, IPIV, WORK, INFO )
 111 *
 112 *  -- LAPACK computational routine (version 3.4.0) --
 113 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 114 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 115 *     November 2011
 116 *
 117 *     .. Scalar Arguments ..
 118       CHARACTER          UPLO
 119       INTEGER            INFO, N
 120 *     ..
 121 *     .. Array Arguments ..
 122       INTEGER            IPIV( * )
 123       COMPLEX*16         AP( * ), WORK( * )
 124 *     ..
 125 *
 126 *  =====================================================================
 127 *
 128 *     .. Parameters ..
 129       DOUBLE PRECISION   ONE
 130       COMPLEX*16         CONE, ZERO
 131       PARAMETER          ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
 132      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
 133 *     ..
 134 *     .. Local Scalars ..
 135       LOGICAL            UPPER
 136       INTEGER            J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
 137       DOUBLE PRECISION   AK, AKP1, D, T
 138       COMPLEX*16         AKKP1, TEMP
 139 *     ..
 140 *     .. External Functions ..
 141       LOGICAL            LSAME
 142       COMPLEX*16         ZDOTC
 143       EXTERNAL           LSAME, ZDOTC
 144 *     ..
 145 *     .. External Subroutines ..
 146       EXTERNAL           XERBLA, ZCOPY, ZHPMV, ZSWAP
 147 *     ..
 148 *     .. Intrinsic Functions ..
 149       INTRINSIC          ABS, DBLE, DCONJG
 150 *     ..
 151 *     .. Executable Statements ..
 152 *
 153 *     Test the input parameters.
 154 *
 155       INFO = 0
 156       UPPER = LSAME( UPLO, 'U' )
 157       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 158          INFO = -1
 159       ELSE IF( N.LT.0 ) THEN
 160          INFO = -2
 161       END IF
 162       IF( INFO.NE.0 ) THEN
 163          CALL XERBLA( 'ZHPTRI', -INFO )
 164          RETURN
 165       END IF
 166 *
 167 *     Quick return if possible
 168 *
 169       IF( N.EQ.0 )
 170      $   RETURN
 171 *
 172 *     Check that the diagonal matrix D is nonsingular.
 173 *
 174       IF( UPPER ) THEN
 175 *
 176 *        Upper triangular storage: examine D from bottom to top
 177 *
 178          KP = N*( N+1 ) / 2
 179          DO 10 INFO = N, 1, -1
 180             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
 181      $         RETURN
 182             KP = KP - INFO
 183    10    CONTINUE
 184       ELSE
 185 *
 186 *        Lower triangular storage: examine D from top to bottom.
 187 *
 188          KP = 1
 189          DO 20 INFO = 1, N
 190             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
 191      $         RETURN
 192             KP = KP + N - INFO + 1
 193    20    CONTINUE
 194       END IF
 195       INFO = 0
 196 *
 197       IF( UPPER ) THEN
 198 *
 199 *        Compute inv(A) from the factorization A = U*D*U**H.
 200 *
 201 *        K is the main loop index, increasing from 1 to N in steps of
 202 *        1 or 2, depending on the size of the diagonal blocks.
 203 *
 204          K = 1
 205          KC = 1
 206    30    CONTINUE
 207 *
 208 *        If K > N, exit from loop.
 209 *
 210          IF( K.GT.N )
 211      $      GO TO 50
 212 *
 213          KCNEXT = KC + K
 214          IF( IPIV( K ).GT.0 ) THEN
 215 *
 216 *           1 x 1 diagonal block
 217 *
 218 *           Invert the diagonal block.
 219 *
 220             AP( KC+K-1 ) = ONE / DBLE( AP( KC+K-1 ) )
 221 *
 222 *           Compute column K of the inverse.
 223 *
 224             IF( K.GT.1 ) THEN
 225                CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
 226                CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
 227      $                     AP( KC ), 1 )
 228                AP( KC+K-1 ) = AP( KC+K-1 ) -
 229      $                        DBLE( ZDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
 230             END IF
 231             KSTEP = 1
 232          ELSE
 233 *
 234 *           2 x 2 diagonal block
 235 *
 236 *           Invert the diagonal block.
 237 *
 238             T = ABS( AP( KCNEXT+K-1 ) )
 239             AK = DBLE( AP( KC+K-1 ) ) / T
 240             AKP1 = DBLE( AP( KCNEXT+K ) ) / T
 241             AKKP1 = AP( KCNEXT+K-1 ) / T
 242             D = T*( AK*AKP1-ONE )
 243             AP( KC+K-1 ) = AKP1 / D
 244             AP( KCNEXT+K ) = AK / D
 245             AP( KCNEXT+K-1 ) = -AKKP1 / D
 246 *
 247 *           Compute columns K and K+1 of the inverse.
 248 *
 249             IF( K.GT.1 ) THEN
 250                CALL ZCOPY( K-1, AP( KC ), 1, WORK, 1 )
 251                CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
 252      $                     AP( KC ), 1 )
 253                AP( KC+K-1 ) = AP( KC+K-1 ) -
 254      $                        DBLE( ZDOTC( K-1, WORK, 1, AP( KC ), 1 ) )
 255                AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
 256      $                            ZDOTC( K-1, AP( KC ), 1, AP( KCNEXT ),
 257      $                            1 )
 258                CALL ZCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
 259                CALL ZHPMV( UPLO, K-1, -CONE, AP, WORK, 1, ZERO,
 260      $                     AP( KCNEXT ), 1 )
 261                AP( KCNEXT+K ) = AP( KCNEXT+K ) -
 262      $                          DBLE( ZDOTC( K-1, WORK, 1, AP( KCNEXT ),
 263      $                          1 ) )
 264             END IF
 265             KSTEP = 2
 266             KCNEXT = KCNEXT + K + 1
 267          END IF
 268 *
 269          KP = ABS( IPIV( K ) )
 270          IF( KP.NE.K ) THEN
 271 *
 272 *           Interchange rows and columns K and KP in the leading
 273 *           submatrix A(1:k+1,1:k+1)
 274 *
 275             KPC = ( KP-1 )*KP / 2 + 1
 276             CALL ZSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
 277             KX = KPC + KP - 1
 278             DO 40 J = KP + 1, K - 1
 279                KX = KX + J - 1
 280                TEMP = DCONJG( AP( KC+J-1 ) )
 281                AP( KC+J-1 ) = DCONJG( AP( KX ) )
 282                AP( KX ) = TEMP
 283    40       CONTINUE
 284             AP( KC+KP-1 ) = DCONJG( AP( KC+KP-1 ) )
 285             TEMP = AP( KC+K-1 )
 286             AP( KC+K-1 ) = AP( KPC+KP-1 )
 287             AP( KPC+KP-1 ) = TEMP
 288             IF( KSTEP.EQ.2 ) THEN
 289                TEMP = AP( KC+K+K-1 )
 290                AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
 291                AP( KC+K+KP-1 ) = TEMP
 292             END IF
 293          END IF
 294 *
 295          K = K + KSTEP
 296          KC = KCNEXT
 297          GO TO 30
 298    50    CONTINUE
 299 *
 300       ELSE
 301 *
 302 *        Compute inv(A) from the factorization A = L*D*L**H.
 303 *
 304 *        K is the main loop index, increasing from 1 to N in steps of
 305 *        1 or 2, depending on the size of the diagonal blocks.
 306 *
 307          NPP = N*( N+1 ) / 2
 308          K = N
 309          KC = NPP
 310    60    CONTINUE
 311 *
 312 *        If K < 1, exit from loop.
 313 *
 314          IF( K.LT.1 )
 315      $      GO TO 80
 316 *
 317          KCNEXT = KC - ( N-K+2 )
 318          IF( IPIV( K ).GT.0 ) THEN
 319 *
 320 *           1 x 1 diagonal block
 321 *
 322 *           Invert the diagonal block.
 323 *
 324             AP( KC ) = ONE / DBLE( AP( KC ) )
 325 *
 326 *           Compute column K of the inverse.
 327 *
 328             IF( K.LT.N ) THEN
 329                CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
 330                CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+N-K+1 ), WORK, 1,
 331      $                     ZERO, AP( KC+1 ), 1 )
 332                AP( KC ) = AP( KC ) - DBLE( ZDOTC( N-K, WORK, 1,
 333      $                    AP( KC+1 ), 1 ) )
 334             END IF
 335             KSTEP = 1
 336          ELSE
 337 *
 338 *           2 x 2 diagonal block
 339 *
 340 *           Invert the diagonal block.
 341 *
 342             T = ABS( AP( KCNEXT+1 ) )
 343             AK = DBLE( AP( KCNEXT ) ) / T
 344             AKP1 = DBLE( AP( KC ) ) / T
 345             AKKP1 = AP( KCNEXT+1 ) / T
 346             D = T*( AK*AKP1-ONE )
 347             AP( KCNEXT ) = AKP1 / D
 348             AP( KC ) = AK / D
 349             AP( KCNEXT+1 ) = -AKKP1 / D
 350 *
 351 *           Compute columns K-1 and K of the inverse.
 352 *
 353             IF( K.LT.N ) THEN
 354                CALL ZCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
 355                CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
 356      $                     1, ZERO, AP( KC+1 ), 1 )
 357                AP( KC ) = AP( KC ) - DBLE( ZDOTC( N-K, WORK, 1,
 358      $                    AP( KC+1 ), 1 ) )
 359                AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
 360      $                          ZDOTC( N-K, AP( KC+1 ), 1,
 361      $                          AP( KCNEXT+2 ), 1 )
 362                CALL ZCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
 363                CALL ZHPMV( UPLO, N-K, -CONE, AP( KC+( N-K+1 ) ), WORK,
 364      $                     1, ZERO, AP( KCNEXT+2 ), 1 )
 365                AP( KCNEXT ) = AP( KCNEXT ) -
 366      $                        DBLE( ZDOTC( N-K, WORK, 1, AP( KCNEXT+2 ),
 367      $                        1 ) )
 368             END IF
 369             KSTEP = 2
 370             KCNEXT = KCNEXT - ( N-K+3 )
 371          END IF
 372 *
 373          KP = ABS( IPIV( K ) )
 374          IF( KP.NE.K ) THEN
 375 *
 376 *           Interchange rows and columns K and KP in the trailing
 377 *           submatrix A(k-1:n,k-1:n)
 378 *
 379             KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
 380             IF( KP.LT.N )
 381      $         CALL ZSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
 382             KX = KC + KP - K
 383             DO 70 J = K + 1, KP - 1
 384                KX = KX + N - J + 1
 385                TEMP = DCONJG( AP( KC+J-K ) )
 386                AP( KC+J-K ) = DCONJG( AP( KX ) )
 387                AP( KX ) = TEMP
 388    70       CONTINUE
 389             AP( KC+KP-K ) = DCONJG( AP( KC+KP-K ) )
 390             TEMP = AP( KC )
 391             AP( KC ) = AP( KPC )
 392             AP( KPC ) = TEMP
 393             IF( KSTEP.EQ.2 ) THEN
 394                TEMP = AP( KC-N+K-1 )
 395                AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
 396                AP( KC-N+KP-1 ) = TEMP
 397             END IF
 398          END IF
 399 *
 400          K = K - KSTEP
 401          KC = KCNEXT
 402          GO TO 60
 403    80    CONTINUE
 404       END IF
 405 *
 406       RETURN
 407 *
 408 *     End of ZHPTRI
 409 *
 410       END