Annotation of rpl/lapack/lapack/zsytri_rook.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZSYTRI_ROOK
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZSYTRI_ROOK + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsytri_rook.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsytri_rook.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsytri_rook.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZSYTRI_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: *> ZSYTRI_ROOK computes the inverse of a complex symmetric
! 39: *> matrix A using the factorization A = U*D*U**T or A = L*D*L**T
! 40: *> computed by ZSYTRF_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**T;
! 52: *> = 'L': Lower triangular, form is A = L*D*L**T.
! 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 ZSYTRF_ROOK.
! 66: *>
! 67: *> On exit, if INFO = 0, the (symmetric) 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 ZSYTRF_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 2011
! 111: *
! 112: *> \ingroup complex16SYcomputational
! 113: *
! 114: *> \par Contributors:
! 115: * ==================
! 116: *>
! 117: *> \verbatim
! 118: *>
! 119: *> November 2011, 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: *>
! 127: *> \endverbatim
! 128: *
! 129: * =====================================================================
! 130: SUBROUTINE ZSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
! 131: *
! 132: * -- LAPACK computational routine (version 3.4.0) --
! 133: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 134: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 135: * November 2011
! 136: *
! 137: * .. Scalar Arguments ..
! 138: CHARACTER UPLO
! 139: INTEGER INFO, LDA, N
! 140: * ..
! 141: * .. Array Arguments ..
! 142: INTEGER IPIV( * )
! 143: COMPLEX*16 A( LDA, * ), WORK( * )
! 144: * ..
! 145: *
! 146: * =====================================================================
! 147: *
! 148: * .. Parameters ..
! 149: COMPLEX*16 CONE, CZERO
! 150: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
! 151: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
! 152: * ..
! 153: * .. Local Scalars ..
! 154: LOGICAL UPPER
! 155: INTEGER K, KP, KSTEP
! 156: COMPLEX*16 AK, AKKP1, AKP1, D, T, TEMP
! 157: * ..
! 158: * .. External Functions ..
! 159: LOGICAL LSAME
! 160: COMPLEX*16 ZDOTU
! 161: EXTERNAL LSAME, ZDOTU
! 162: * ..
! 163: * .. External Subroutines ..
! 164: EXTERNAL ZCOPY, ZSWAP, ZSYMV, XERBLA
! 165: * ..
! 166: * .. Intrinsic Functions ..
! 167: INTRINSIC MAX
! 168: * ..
! 169: * .. Executable Statements ..
! 170: *
! 171: * Test the input parameters.
! 172: *
! 173: INFO = 0
! 174: UPPER = LSAME( UPLO, 'U' )
! 175: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 176: INFO = -1
! 177: ELSE IF( N.LT.0 ) THEN
! 178: INFO = -2
! 179: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 180: INFO = -4
! 181: END IF
! 182: IF( INFO.NE.0 ) THEN
! 183: CALL XERBLA( 'ZSYTRI_ROOK', -INFO )
! 184: RETURN
! 185: END IF
! 186: *
! 187: * Quick return if possible
! 188: *
! 189: IF( N.EQ.0 )
! 190: $ RETURN
! 191: *
! 192: * Check that the diagonal matrix D is nonsingular.
! 193: *
! 194: IF( UPPER ) THEN
! 195: *
! 196: * Upper triangular storage: examine D from bottom to top
! 197: *
! 198: DO 10 INFO = N, 1, -1
! 199: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
! 200: $ RETURN
! 201: 10 CONTINUE
! 202: ELSE
! 203: *
! 204: * Lower triangular storage: examine D from top to bottom.
! 205: *
! 206: DO 20 INFO = 1, N
! 207: IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
! 208: $ RETURN
! 209: 20 CONTINUE
! 210: END IF
! 211: INFO = 0
! 212: *
! 213: IF( UPPER ) THEN
! 214: *
! 215: * Compute inv(A) from the factorization A = U*D*U**T.
! 216: *
! 217: * K is the main loop index, increasing from 1 to N in steps of
! 218: * 1 or 2, depending on the size of the diagonal blocks.
! 219: *
! 220: K = 1
! 221: 30 CONTINUE
! 222: *
! 223: * If K > N, exit from loop.
! 224: *
! 225: IF( K.GT.N )
! 226: $ GO TO 40
! 227: *
! 228: IF( IPIV( K ).GT.0 ) THEN
! 229: *
! 230: * 1 x 1 diagonal block
! 231: *
! 232: * Invert the diagonal block.
! 233: *
! 234: A( K, K ) = CONE / A( K, K )
! 235: *
! 236: * Compute column K of the inverse.
! 237: *
! 238: IF( K.GT.1 ) THEN
! 239: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
! 240: CALL ZSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
! 241: $ A( 1, K ), 1 )
! 242: A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
! 243: $ 1 )
! 244: END IF
! 245: KSTEP = 1
! 246: ELSE
! 247: *
! 248: * 2 x 2 diagonal block
! 249: *
! 250: * Invert the diagonal block.
! 251: *
! 252: T = A( K, K+1 )
! 253: AK = A( K, K ) / T
! 254: AKP1 = A( K+1, K+1 ) / T
! 255: AKKP1 = A( K, K+1 ) / T
! 256: D = T*( AK*AKP1-CONE )
! 257: A( K, K ) = AKP1 / D
! 258: A( K+1, K+1 ) = AK / D
! 259: A( K, K+1 ) = -AKKP1 / D
! 260: *
! 261: * Compute columns K and K+1 of the inverse.
! 262: *
! 263: IF( K.GT.1 ) THEN
! 264: CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
! 265: CALL ZSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
! 266: $ A( 1, K ), 1 )
! 267: A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
! 268: $ 1 )
! 269: A( K, K+1 ) = A( K, K+1 ) -
! 270: $ ZDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
! 271: CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
! 272: CALL ZSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
! 273: $ A( 1, K+1 ), 1 )
! 274: A( K+1, K+1 ) = A( K+1, K+1 ) -
! 275: $ ZDOTU( K-1, WORK, 1, A( 1, K+1 ), 1 )
! 276: END IF
! 277: KSTEP = 2
! 278: END IF
! 279: *
! 280: IF( KSTEP.EQ.1 ) THEN
! 281: *
! 282: * Interchange rows and columns K and IPIV(K) in the leading
! 283: * submatrix A(1:k+1,1:k+1)
! 284: *
! 285: KP = IPIV( K )
! 286: IF( KP.NE.K ) THEN
! 287: IF( KP.GT.1 )
! 288: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
! 289: CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
! 290: TEMP = A( K, K )
! 291: A( K, K ) = A( KP, KP )
! 292: A( KP, KP ) = TEMP
! 293: END IF
! 294: ELSE
! 295: *
! 296: * Interchange rows and columns K and K+1 with -IPIV(K) and
! 297: * -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
! 298: *
! 299: KP = -IPIV( K )
! 300: IF( KP.NE.K ) THEN
! 301: IF( KP.GT.1 )
! 302: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
! 303: CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
! 304: *
! 305: TEMP = A( K, K )
! 306: A( K, K ) = A( KP, KP )
! 307: A( KP, KP ) = TEMP
! 308: TEMP = A( K, K+1 )
! 309: A( K, K+1 ) = A( KP, K+1 )
! 310: A( KP, K+1 ) = TEMP
! 311: END IF
! 312: *
! 313: K = K + 1
! 314: KP = -IPIV( K )
! 315: IF( KP.NE.K ) THEN
! 316: IF( KP.GT.1 )
! 317: $ CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
! 318: CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
! 319: TEMP = A( K, K )
! 320: A( K, K ) = A( KP, KP )
! 321: A( KP, KP ) = TEMP
! 322: END IF
! 323: END IF
! 324: *
! 325: K = K + 1
! 326: GO TO 30
! 327: 40 CONTINUE
! 328: *
! 329: ELSE
! 330: *
! 331: * Compute inv(A) from the factorization A = L*D*L**T.
! 332: *
! 333: * K is the main loop index, increasing from 1 to N in steps of
! 334: * 1 or 2, depending on the size of the diagonal blocks.
! 335: *
! 336: K = N
! 337: 50 CONTINUE
! 338: *
! 339: * If K < 1, exit from loop.
! 340: *
! 341: IF( K.LT.1 )
! 342: $ GO TO 60
! 343: *
! 344: IF( IPIV( K ).GT.0 ) THEN
! 345: *
! 346: * 1 x 1 diagonal block
! 347: *
! 348: * Invert the diagonal block.
! 349: *
! 350: A( K, K ) = CONE / A( K, K )
! 351: *
! 352: * Compute column K of the inverse.
! 353: *
! 354: IF( K.LT.N ) THEN
! 355: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
! 356: CALL ZSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
! 357: $ CZERO, A( K+1, K ), 1 )
! 358: A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
! 359: $ 1 )
! 360: END IF
! 361: KSTEP = 1
! 362: ELSE
! 363: *
! 364: * 2 x 2 diagonal block
! 365: *
! 366: * Invert the diagonal block.
! 367: *
! 368: T = A( K, K-1 )
! 369: AK = A( K-1, K-1 ) / T
! 370: AKP1 = A( K, K ) / T
! 371: AKKP1 = A( K, K-1 ) / T
! 372: D = T*( AK*AKP1-CONE )
! 373: A( K-1, K-1 ) = AKP1 / D
! 374: A( K, K ) = AK / D
! 375: A( K, K-1 ) = -AKKP1 / D
! 376: *
! 377: * Compute columns K-1 and K of the inverse.
! 378: *
! 379: IF( K.LT.N ) THEN
! 380: CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
! 381: CALL ZSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
! 382: $ CZERO, A( K+1, K ), 1 )
! 383: A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
! 384: $ 1 )
! 385: A( K, K-1 ) = A( K, K-1 ) -
! 386: $ ZDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
! 387: $ 1 )
! 388: CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
! 389: CALL ZSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
! 390: $ CZERO, A( K+1, K-1 ), 1 )
! 391: A( K-1, K-1 ) = A( K-1, K-1 ) -
! 392: $ ZDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 )
! 393: END IF
! 394: KSTEP = 2
! 395: END IF
! 396: *
! 397: IF( KSTEP.EQ.1 ) THEN
! 398: *
! 399: * Interchange rows and columns K and IPIV(K) in the trailing
! 400: * submatrix A(k-1:n,k-1:n)
! 401: *
! 402: KP = IPIV( K )
! 403: IF( KP.NE.K ) THEN
! 404: IF( KP.LT.N )
! 405: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
! 406: CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
! 407: TEMP = A( K, K )
! 408: A( K, K ) = A( KP, KP )
! 409: A( KP, KP ) = TEMP
! 410: END IF
! 411: ELSE
! 412: *
! 413: * Interchange rows and columns K and K-1 with -IPIV(K) and
! 414: * -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
! 415: *
! 416: KP = -IPIV( K )
! 417: IF( KP.NE.K ) THEN
! 418: IF( KP.LT.N )
! 419: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
! 420: CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
! 421: *
! 422: TEMP = A( K, K )
! 423: A( K, K ) = A( KP, KP )
! 424: A( KP, KP ) = TEMP
! 425: TEMP = A( K, K-1 )
! 426: A( K, K-1 ) = A( KP, K-1 )
! 427: A( KP, K-1 ) = TEMP
! 428: END IF
! 429: *
! 430: K = K - 1
! 431: KP = -IPIV( K )
! 432: IF( KP.NE.K ) THEN
! 433: IF( KP.LT.N )
! 434: $ CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
! 435: CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
! 436: TEMP = A( K, K )
! 437: A( K, K ) = A( KP, KP )
! 438: A( KP, KP ) = TEMP
! 439: END IF
! 440: END IF
! 441: *
! 442: K = K - 1
! 443: GO TO 50
! 444: 60 CONTINUE
! 445: END IF
! 446: *
! 447: RETURN
! 448: *
! 449: * End of ZSYTRI_ROOK
! 450: *
! 451: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zsytri_rook.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack