Annotation of rpl/lapack/lapack/zhetri_rook.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZHETRI_ROOK computes the inverse of HE matrix using the factorization obtained with the bounded Bunch-Kaufman ("rook") diagonal pivoting method.
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZHETRI_ROOK + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetri_rook.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetri_rook.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetri_rook.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZHETRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * CHARACTER UPLO
! 25: * INTEGER INFO, LDA, N
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * INTEGER IPIV( * )
! 29: * COMPLEX*16 A( LDA, * ), WORK( * )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> ZHETRI_ROOK computes the inverse of a complex Hermitian indefinite matrix
! 39: *> A using the factorization A = U*D*U**H or A = L*D*L**H computed by
! 40: *> ZHETRF_ROOK.
! 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] A
! 62: *> \verbatim
! 63: *> A is COMPLEX*16 array, dimension (LDA,N)
! 64: *> On entry, the block diagonal matrix D and the multipliers
! 65: *> used to obtain the factor U or L as computed by ZHETRF_ROOK.
! 66: *>
! 67: *> On exit, if INFO = 0, the (Hermitian) inverse of the original
! 68: *> matrix. If UPLO = 'U', the upper triangular part of the
! 69: *> inverse is formed and the part of A below the diagonal is not
! 70: *> referenced; if UPLO = 'L' the lower triangular part of the
! 71: *> inverse is formed and the part of A above the diagonal is
! 72: *> not referenced.
! 73: *> \endverbatim
! 74: *>
! 75: *> \param[in] LDA
! 76: *> \verbatim
! 77: *> LDA is INTEGER
! 78: *> The leading dimension of the array A. LDA >= max(1,N).
! 79: *> \endverbatim
! 80: *>
! 81: *> \param[in] IPIV
! 82: *> \verbatim
! 83: *> IPIV is INTEGER array, dimension (N)
! 84: *> Details of the interchanges and the block structure of D
! 85: *> as determined by ZHETRF_ROOK.
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[out] WORK
! 89: *> \verbatim
! 90: *> WORK is COMPLEX*16 array, dimension (N)
! 91: *> \endverbatim
! 92: *>
! 93: *> \param[out] INFO
! 94: *> \verbatim
! 95: *> INFO is INTEGER
! 96: *> = 0: successful exit
! 97: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 98: *> > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
! 99: *> inverse could not be computed.
! 100: *> \endverbatim
! 101: *
! 102: * Authors:
! 103: * ========
! 104: *
! 105: *> \author Univ. of Tennessee
! 106: *> \author Univ. of California Berkeley
! 107: *> \author Univ. of Colorado Denver
! 108: *> \author NAG Ltd.
! 109: *
! 110: *> \date November 2013
! 111: *
! 112: *> \ingroup complex16HEcomputational
! 113: *
! 114: *> \par Contributors:
! 115: * ==================
! 116: *>
! 117: *> \verbatim
! 118: *>
! 119: *> November 2013, Igor Kozachenko,
! 120: *> Computer Science Division,
! 121: *> University of California, Berkeley
! 122: *>
! 123: *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
! 124: *> School of Mathematics,
! 125: *> University of Manchester
! 126: *> \endverbatim
! 127: *
! 128: * =====================================================================
! 129: SUBROUTINE ZHETRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
! 130: *
! 131: * -- LAPACK computational routine (version 3.5.0) --
! 132: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 133: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 134: * November 2013
! 135: *
! 136: * .. Scalar Arguments ..
! 137: CHARACTER UPLO
! 138: INTEGER INFO, LDA, N
! 139: * ..
! 140: * .. Array Arguments ..
! 141: INTEGER IPIV( * )
! 142: COMPLEX*16 A( LDA, * ), WORK( * )
! 143: * ..
! 144: *
! 145: * =====================================================================
! 146: *
! 147: * .. Parameters ..
! 148: DOUBLE PRECISION ONE
! 149: COMPLEX*16 CONE, CZERO
! 150: PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ),
! 151: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
! 152: * ..
! 153: * .. Local Scalars ..
! 154: LOGICAL UPPER
! 155: INTEGER J, K, KP, KSTEP
! 156: DOUBLE PRECISION AK, AKP1, D, T
! 157: COMPLEX*16 AKKP1, TEMP
! 158: * ..
! 159: * .. External Functions ..
! 160: LOGICAL LSAME
! 161: COMPLEX*16 ZDOTC
! 162: EXTERNAL LSAME, ZDOTC
! 163: * ..
! 164: * .. External Subroutines ..
! 165: EXTERNAL ZCOPY, ZHEMV, ZSWAP, XERBLA
! 166: * ..
! 167: * .. Intrinsic Functions ..
! 168: INTRINSIC ABS, DCONJG, MAX, DBLE
! 169: * ..
! 170: * .. Executable Statements ..
! 171: *
! 172: * Test the input parameters.
! 173: *
! 174: INFO = 0
! 175: UPPER = LSAME( UPLO, 'U' )
! 176: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 177: INFO = -1
! 178: ELSE IF( N.LT.0 ) THEN
! 179: INFO = -2
! 180: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 181: INFO = -4
! 182: END IF
! 183: IF( INFO.NE.0 ) THEN
! 184: CALL XERBLA( 'ZHETRI_ROOK', -INFO )
! 185: RETURN
! 186: END IF
! 187: *
! 188: * Quick return if possible
! 189: *
! 190: IF( N.EQ.0 )
! 191: $ RETURN
! 192: *
! 193: * Check that the diagonal matrix D is nonsingular.
! 194: *
! 195: IF( UPPER ) THEN
! 196: *
! 197: * Upper triangular storage: examine D from bottom to top
! 198: *
! 199: DO 10 INFO = N, 1, -1
! 200: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
! 201: $ RETURN
! 202: 10 CONTINUE
! 203: ELSE
! 204: *
! 205: * Lower triangular storage: examine D from top to bottom.
! 206: *
! 207: DO 20 INFO = 1, N
! 208: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
! 209: $ RETURN
! 210: 20 CONTINUE
! 211: END IF
! 212: INFO = 0
! 213: *
! 214: IF( UPPER ) THEN
! 215: *
! 216: * Compute inv(A) from the factorization A = U*D*U**H.
! 217: *
! 218: * K is the main loop index, increasing from 1 to N in steps of
! 219: * 1 or 2, depending on the size of the diagonal blocks.
! 220: *
! 221: K = 1
! 222: 30 CONTINUE
! 223: *
! 224: * If K > N, exit from loop.
! 225: *
! 226: IF( K.GT.N )
! 227: $ GO TO 70
! 228: *
! 229: IF( IPIV( K ).GT.0 ) THEN
! 230: *
! 231: * 1 x 1 diagonal block
! 232: *
! 233: * Invert the diagonal block.
! 234: *
! 235: A( K, K ) = ONE / DBLE( A( K, K ) )
! 236: *
! 237: * Compute column K of the inverse.
! 238: *
! 239: IF( K.GT.1 ) THEN
! 240: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
! 241: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
! 242: $ A( 1, K ), 1 )
! 243: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
! 244: $ K ), 1 ) )
! 245: END IF
! 246: KSTEP = 1
! 247: ELSE
! 248: *
! 249: * 2 x 2 diagonal block
! 250: *
! 251: * Invert the diagonal block.
! 252: *
! 253: T = ABS( A( K, K+1 ) )
! 254: AK = DBLE( A( K, K ) ) / T
! 255: AKP1 = DBLE( A( K+1, K+1 ) ) / T
! 256: AKKP1 = A( K, K+1 ) / T
! 257: D = T*( AK*AKP1-ONE )
! 258: A( K, K ) = AKP1 / D
! 259: A( K+1, K+1 ) = AK / D
! 260: A( K, K+1 ) = -AKKP1 / D
! 261: *
! 262: * Compute columns K and K+1 of the inverse.
! 263: *
! 264: IF( K.GT.1 ) THEN
! 265: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
! 266: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
! 267: $ A( 1, K ), 1 )
! 268: A( K, K ) = A( K, K ) - DBLE( ZDOTC( K-1, WORK, 1, A( 1,
! 269: $ K ), 1 ) )
! 270: A( K, K+1 ) = A( K, K+1 ) -
! 271: $ ZDOTC( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
! 272: CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
! 273: CALL ZHEMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
! 274: $ A( 1, K+1 ), 1 )
! 275: A( K+1, K+1 ) = A( K+1, K+1 ) -
! 276: $ DBLE( ZDOTC( K-1, WORK, 1, A( 1, K+1 ),
! 277: $ 1 ) )
! 278: END IF
! 279: KSTEP = 2
! 280: END IF
! 281: *
! 282: IF( KSTEP.EQ.1 ) THEN
! 283: *
! 284: * Interchange rows and columns K and IPIV(K) in the leading
! 285: * submatrix A(1:k,1:k)
! 286: *
! 287: KP = IPIV( K )
! 288: IF( KP.NE.K ) THEN
! 289: *
! 290: IF( KP.GT.1 )
! 291: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
! 292: *
! 293: DO 40 J = KP + 1, K - 1
! 294: TEMP = DCONJG( A( J, K ) )
! 295: A( J, K ) = DCONJG( A( KP, J ) )
! 296: A( KP, J ) = TEMP
! 297: 40 CONTINUE
! 298: *
! 299: A( KP, K ) = DCONJG( A( KP, K ) )
! 300: *
! 301: TEMP = A( K, K )
! 302: A( K, K ) = A( KP, KP )
! 303: A( KP, KP ) = TEMP
! 304: END IF
! 305: ELSE
! 306: *
! 307: * Interchange rows and columns K and K+1 with -IPIV(K) and
! 308: * -IPIV(K+1) in the leading submatrix A(k+1:n,k+1:n)
! 309: *
! 310: * (1) Interchange rows and columns K and -IPIV(K)
! 311: *
! 312: KP = -IPIV( K )
! 313: IF( KP.NE.K ) THEN
! 314: *
! 315: IF( KP.GT.1 )
! 316: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
! 317: *
! 318: DO 50 J = KP + 1, K - 1
! 319: TEMP = DCONJG( A( J, K ) )
! 320: A( J, K ) = DCONJG( A( KP, J ) )
! 321: A( KP, J ) = TEMP
! 322: 50 CONTINUE
! 323: *
! 324: A( KP, K ) = DCONJG( A( KP, K ) )
! 325: *
! 326: TEMP = A( K, K )
! 327: A( K, K ) = A( KP, KP )
! 328: A( KP, KP ) = TEMP
! 329: *
! 330: TEMP = A( K, K+1 )
! 331: A( K, K+1 ) = A( KP, K+1 )
! 332: A( KP, K+1 ) = TEMP
! 333: END IF
! 334: *
! 335: * (2) Interchange rows and columns K+1 and -IPIV(K+1)
! 336: *
! 337: K = K + 1
! 338: KP = -IPIV( K )
! 339: IF( KP.NE.K ) THEN
! 340: *
! 341: IF( KP.GT.1 )
! 342: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
! 343: *
! 344: DO 60 J = KP + 1, K - 1
! 345: TEMP = DCONJG( A( J, K ) )
! 346: A( J, K ) = DCONJG( A( KP, J ) )
! 347: A( KP, J ) = TEMP
! 348: 60 CONTINUE
! 349: *
! 350: A( KP, K ) = DCONJG( A( KP, K ) )
! 351: *
! 352: TEMP = A( K, K )
! 353: A( K, K ) = A( KP, KP )
! 354: A( KP, KP ) = TEMP
! 355: END IF
! 356: END IF
! 357: *
! 358: K = K + 1
! 359: GO TO 30
! 360: 70 CONTINUE
! 361: *
! 362: ELSE
! 363: *
! 364: * Compute inv(A) from the factorization A = L*D*L**H.
! 365: *
! 366: * K is the main loop index, decreasing from N to 1 in steps of
! 367: * 1 or 2, depending on the size of the diagonal blocks.
! 368: *
! 369: K = N
! 370: 80 CONTINUE
! 371: *
! 372: * If K < 1, exit from loop.
! 373: *
! 374: IF( K.LT.1 )
! 375: $ GO TO 120
! 376: *
! 377: IF( IPIV( K ).GT.0 ) THEN
! 378: *
! 379: * 1 x 1 diagonal block
! 380: *
! 381: * Invert the diagonal block.
! 382: *
! 383: A( K, K ) = ONE / DBLE( A( K, K ) )
! 384: *
! 385: * Compute column K of the inverse.
! 386: *
! 387: IF( K.LT.N ) THEN
! 388: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
! 389: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
! 390: $ 1, CZERO, A( K+1, K ), 1 )
! 391: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
! 392: $ A( K+1, K ), 1 ) )
! 393: END IF
! 394: KSTEP = 1
! 395: ELSE
! 396: *
! 397: * 2 x 2 diagonal block
! 398: *
! 399: * Invert the diagonal block.
! 400: *
! 401: T = ABS( A( K, K-1 ) )
! 402: AK = DBLE( A( K-1, K-1 ) ) / T
! 403: AKP1 = DBLE( A( K, K ) ) / T
! 404: AKKP1 = A( K, K-1 ) / T
! 405: D = T*( AK*AKP1-ONE )
! 406: A( K-1, K-1 ) = AKP1 / D
! 407: A( K, K ) = AK / D
! 408: A( K, K-1 ) = -AKKP1 / D
! 409: *
! 410: * Compute columns K-1 and K of the inverse.
! 411: *
! 412: IF( K.LT.N ) THEN
! 413: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
! 414: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
! 415: $ 1, CZERO, A( K+1, K ), 1 )
! 416: A( K, K ) = A( K, K ) - DBLE( ZDOTC( N-K, WORK, 1,
! 417: $ A( K+1, K ), 1 ) )
! 418: A( K, K-1 ) = A( K, K-1 ) -
! 419: $ ZDOTC( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
! 420: $ 1 )
! 421: CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
! 422: CALL ZHEMV( UPLO, N-K, -CONE, A( K+1, K+1 ), LDA, WORK,
! 423: $ 1, CZERO, A( K+1, K-1 ), 1 )
! 424: A( K-1, K-1 ) = A( K-1, K-1 ) -
! 425: $ DBLE( ZDOTC( N-K, WORK, 1, A( K+1, K-1 ),
! 426: $ 1 ) )
! 427: END IF
! 428: KSTEP = 2
! 429: END IF
! 430: *
! 431: IF( KSTEP.EQ.1 ) THEN
! 432: *
! 433: * Interchange rows and columns K and IPIV(K) in the trailing
! 434: * submatrix A(k:n,k:n)
! 435: *
! 436: KP = IPIV( K )
! 437: IF( KP.NE.K ) THEN
! 438: *
! 439: IF( KP.LT.N )
! 440: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
! 441: *
! 442: DO 90 J = K + 1, KP - 1
! 443: TEMP = DCONJG( A( J, K ) )
! 444: A( J, K ) = DCONJG( A( KP, J ) )
! 445: A( KP, J ) = TEMP
! 446: 90 CONTINUE
! 447: *
! 448: A( KP, K ) = DCONJG( A( KP, K ) )
! 449: *
! 450: TEMP = A( K, K )
! 451: A( K, K ) = A( KP, KP )
! 452: A( KP, KP ) = TEMP
! 453: END IF
! 454: ELSE
! 455: *
! 456: * Interchange rows and columns K and K-1 with -IPIV(K) and
! 457: * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
! 458: *
! 459: * (1) Interchange rows and columns K and -IPIV(K)
! 460: *
! 461: KP = -IPIV( K )
! 462: IF( KP.NE.K ) THEN
! 463: *
! 464: IF( KP.LT.N )
! 465: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
! 466: *
! 467: DO 100 J = K + 1, KP - 1
! 468: TEMP = DCONJG( A( J, K ) )
! 469: A( J, K ) = DCONJG( A( KP, J ) )
! 470: A( KP, J ) = TEMP
! 471: 100 CONTINUE
! 472: *
! 473: A( KP, K ) = DCONJG( A( KP, K ) )
! 474: *
! 475: TEMP = A( K, K )
! 476: A( K, K ) = A( KP, KP )
! 477: A( KP, KP ) = TEMP
! 478: *
! 479: TEMP = A( K, K-1 )
! 480: A( K, K-1 ) = A( KP, K-1 )
! 481: A( KP, K-1 ) = TEMP
! 482: END IF
! 483: *
! 484: * (2) Interchange rows and columns K-1 and -IPIV(K-1)
! 485: *
! 486: K = K - 1
! 487: KP = -IPIV( K )
! 488: IF( KP.NE.K ) THEN
! 489: *
! 490: IF( KP.LT.N )
! 491: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
! 492: *
! 493: DO 110 J = K + 1, KP - 1
! 494: TEMP = DCONJG( A( J, K ) )
! 495: A( J, K ) = DCONJG( A( KP, J ) )
! 496: A( KP, J ) = TEMP
! 497: 110 CONTINUE
! 498: *
! 499: A( KP, K ) = DCONJG( A( KP, K ) )
! 500: *
! 501: TEMP = A( K, K )
! 502: A( K, K ) = A( KP, KP )
! 503: A( KP, KP ) = TEMP
! 504: END IF
! 505: END IF
! 506: *
! 507: K = K - 1
! 508: GO TO 80
! 509: 120 CONTINUE
! 510: END IF
! 511: *
! 512: RETURN
! 513: *
! 514: * End of ZHETRI_ROOK
! 515: *
! 516: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetri_rook.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack