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