Annotation of rpl/lapack/lapack/dsytrs_rook.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DSYTRS_ROOK
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DSYTRS_ROOK + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytrs_rook.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytrs_rook.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytrs_rook.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DSYTRS_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: * DOUBLE PRECISION A( LDA, * ), B( LDB, * )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> DSYTRS_ROOK solves a system of linear equations A*X = B with
! 39: *> a real symmetric matrix A using the factorization A = U*D*U**T or
! 40: *> A = L*D*L**T computed by DSYTRF_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] 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 DOUBLE PRECISION 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 DSYTRF_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 DSYTRF_ROOK.
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[in,out] B
! 89: *> \verbatim
! 90: *> B is DOUBLE PRECISION 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 April 2012
! 117: *
! 118: *> \ingroup doubleSYcomputational
! 119: *
! 120: *> \par Contributors:
! 121: * ==================
! 122: *>
! 123: *> \verbatim
! 124: *>
! 125: *> April 2012, 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 DSYTRS_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
! 137: $ INFO )
! 138: *
! 139: * -- LAPACK computational routine (version 3.4.1) --
! 140: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 141: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 142: * April 2012
! 143: *
! 144: * .. Scalar Arguments ..
! 145: CHARACTER UPLO
! 146: INTEGER INFO, LDA, LDB, N, NRHS
! 147: * ..
! 148: * .. Array Arguments ..
! 149: INTEGER IPIV( * )
! 150: DOUBLE PRECISION A( LDA, * ), B( LDB, * )
! 151: * ..
! 152: *
! 153: * =====================================================================
! 154: *
! 155: * .. Parameters ..
! 156: DOUBLE PRECISION ONE
! 157: PARAMETER ( ONE = 1.0D+0 )
! 158: * ..
! 159: * .. Local Scalars ..
! 160: LOGICAL UPPER
! 161: INTEGER J, K, KP
! 162: DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
! 163: * ..
! 164: * .. External Functions ..
! 165: LOGICAL LSAME
! 166: EXTERNAL LSAME
! 167: * ..
! 168: * .. External Subroutines ..
! 169: EXTERNAL DGEMV, DGER, DSCAL, DSWAP, XERBLA
! 170: * ..
! 171: * .. Intrinsic Functions ..
! 172: INTRINSIC MAX
! 173: * ..
! 174: * .. Executable Statements ..
! 175: *
! 176: INFO = 0
! 177: UPPER = LSAME( UPLO, 'U' )
! 178: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 179: INFO = -1
! 180: ELSE IF( N.LT.0 ) THEN
! 181: INFO = -2
! 182: ELSE IF( NRHS.LT.0 ) THEN
! 183: INFO = -3
! 184: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 185: INFO = -5
! 186: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 187: INFO = -8
! 188: END IF
! 189: IF( INFO.NE.0 ) THEN
! 190: CALL XERBLA( 'DSYTRS_ROOK', -INFO )
! 191: RETURN
! 192: END IF
! 193: *
! 194: * Quick return if possible
! 195: *
! 196: IF( N.EQ.0 .OR. NRHS.EQ.0 )
! 197: $ RETURN
! 198: *
! 199: IF( UPPER ) THEN
! 200: *
! 201: * Solve A*X = B, where A = U*D*U**T.
! 202: *
! 203: * First solve U*D*X = B, overwriting B with X.
! 204: *
! 205: * K is the main loop index, decreasing from N to 1 in steps of
! 206: * 1 or 2, depending on the size of the diagonal blocks.
! 207: *
! 208: K = N
! 209: 10 CONTINUE
! 210: *
! 211: * If K < 1, exit from loop.
! 212: *
! 213: IF( K.LT.1 )
! 214: $ GO TO 30
! 215: *
! 216: IF( IPIV( K ).GT.0 ) THEN
! 217: *
! 218: * 1 x 1 diagonal block
! 219: *
! 220: * Interchange rows K and IPIV(K).
! 221: *
! 222: KP = IPIV( K )
! 223: IF( KP.NE.K )
! 224: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
! 225: *
! 226: * Multiply by inv(U(K)), where U(K) is the transformation
! 227: * stored in column K of A.
! 228: *
! 229: CALL DGER( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB,
! 230: $ B( 1, 1 ), LDB )
! 231: *
! 232: * Multiply by the inverse of the diagonal block.
! 233: *
! 234: CALL DSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB )
! 235: K = K - 1
! 236: ELSE
! 237: *
! 238: * 2 x 2 diagonal block
! 239: *
! 240: * Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1)
! 241: *
! 242: KP = -IPIV( K )
! 243: IF( KP.NE.K )
! 244: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
! 245: *
! 246: KP = -IPIV( K-1 )
! 247: IF( KP.NE.K-1 )
! 248: $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
! 249: *
! 250: * Multiply by inv(U(K)), where U(K) is the transformation
! 251: * stored in columns K-1 and K of A.
! 252: *
! 253: IF( K.GT.2 ) THEN
! 254: CALL DGER( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ),
! 255: $ LDB, B( 1, 1 ), LDB )
! 256: CALL DGER( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ),
! 257: $ LDB, B( 1, 1 ), LDB )
! 258: END IF
! 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 ) / AKM1K
! 265: DENOM = AKM1*AK - ONE
! 266: DO 20 J = 1, NRHS
! 267: BKM1 = B( K-1, J ) / AKM1K
! 268: BK = B( K, J ) / 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**T *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**T(K)), where U(K) is the transformation
! 296: * stored in column K of A.
! 297: *
! 298: IF( K.GT.1 )
! 299: $ CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B,
! 300: $ LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
! 301: *
! 302: * Interchange rows K and IPIV(K).
! 303: *
! 304: KP = IPIV( K )
! 305: IF( KP.NE.K )
! 306: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
! 307: K = K + 1
! 308: ELSE
! 309: *
! 310: * 2 x 2 diagonal block
! 311: *
! 312: * Multiply by inv(U**T(K+1)), where U(K+1) is the transformation
! 313: * stored in columns K and K+1 of A.
! 314: *
! 315: IF( K.GT.1 ) THEN
! 316: CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B,
! 317: $ LDB, A( 1, K ), 1, ONE, B( K, 1 ), LDB )
! 318: CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B,
! 319: $ LDB, A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB )
! 320: END IF
! 321: *
! 322: * Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1).
! 323: *
! 324: KP = -IPIV( K )
! 325: IF( KP.NE.K )
! 326: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
! 327: *
! 328: KP = -IPIV( K+1 )
! 329: IF( KP.NE.K+1 )
! 330: $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
! 331: *
! 332: K = K + 2
! 333: END IF
! 334: *
! 335: GO TO 40
! 336: 50 CONTINUE
! 337: *
! 338: ELSE
! 339: *
! 340: * Solve A*X = B, where A = L*D*L**T.
! 341: *
! 342: * First solve L*D*X = B, overwriting B with X.
! 343: *
! 344: * K is the main loop index, increasing from 1 to N in steps of
! 345: * 1 or 2, depending on the size of the diagonal blocks.
! 346: *
! 347: K = 1
! 348: 60 CONTINUE
! 349: *
! 350: * If K > N, exit from loop.
! 351: *
! 352: IF( K.GT.N )
! 353: $ GO TO 80
! 354: *
! 355: IF( IPIV( K ).GT.0 ) THEN
! 356: *
! 357: * 1 x 1 diagonal block
! 358: *
! 359: * Interchange rows K and IPIV(K).
! 360: *
! 361: KP = IPIV( K )
! 362: IF( KP.NE.K )
! 363: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
! 364: *
! 365: * Multiply by inv(L(K)), where L(K) is the transformation
! 366: * stored in column K of A.
! 367: *
! 368: IF( K.LT.N )
! 369: $ CALL DGER( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ),
! 370: $ LDB, B( K+1, 1 ), LDB )
! 371: *
! 372: * Multiply by the inverse of the diagonal block.
! 373: *
! 374: CALL DSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB )
! 375: K = K + 1
! 376: ELSE
! 377: *
! 378: * 2 x 2 diagonal block
! 379: *
! 380: * Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1)
! 381: *
! 382: KP = -IPIV( K )
! 383: IF( KP.NE.K )
! 384: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
! 385: *
! 386: KP = -IPIV( K+1 )
! 387: IF( KP.NE.K+1 )
! 388: $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
! 389: *
! 390: * Multiply by inv(L(K)), where L(K) is the transformation
! 391: * stored in columns K and K+1 of A.
! 392: *
! 393: IF( K.LT.N-1 ) THEN
! 394: CALL DGER( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ),
! 395: $ LDB, B( K+2, 1 ), LDB )
! 396: CALL DGER( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1,
! 397: $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
! 398: END IF
! 399: *
! 400: * Multiply by the inverse of the diagonal block.
! 401: *
! 402: AKM1K = A( K+1, K )
! 403: AKM1 = A( K, K ) / AKM1K
! 404: AK = A( K+1, K+1 ) / AKM1K
! 405: DENOM = AKM1*AK - ONE
! 406: DO 70 J = 1, NRHS
! 407: BKM1 = B( K, J ) / AKM1K
! 408: BK = B( K+1, J ) / AKM1K
! 409: B( K, J ) = ( AK*BKM1-BK ) / DENOM
! 410: B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
! 411: 70 CONTINUE
! 412: K = K + 2
! 413: END IF
! 414: *
! 415: GO TO 60
! 416: 80 CONTINUE
! 417: *
! 418: * Next solve L**T *X = B, overwriting B with X.
! 419: *
! 420: * K is the main loop index, decreasing from N to 1 in steps of
! 421: * 1 or 2, depending on the size of the diagonal blocks.
! 422: *
! 423: K = N
! 424: 90 CONTINUE
! 425: *
! 426: * If K < 1, exit from loop.
! 427: *
! 428: IF( K.LT.1 )
! 429: $ GO TO 100
! 430: *
! 431: IF( IPIV( K ).GT.0 ) THEN
! 432: *
! 433: * 1 x 1 diagonal block
! 434: *
! 435: * Multiply by inv(L**T(K)), where L(K) is the transformation
! 436: * stored in column K of A.
! 437: *
! 438: IF( K.LT.N )
! 439: $ CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
! 440: $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
! 441: *
! 442: * Interchange rows K and IPIV(K).
! 443: *
! 444: KP = IPIV( K )
! 445: IF( KP.NE.K )
! 446: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
! 447: K = K - 1
! 448: ELSE
! 449: *
! 450: * 2 x 2 diagonal block
! 451: *
! 452: * Multiply by inv(L**T(K-1)), where L(K-1) is the transformation
! 453: * stored in columns K-1 and K of A.
! 454: *
! 455: IF( K.LT.N ) THEN
! 456: CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
! 457: $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
! 458: CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
! 459: $ LDB, A( K+1, K-1 ), 1, ONE, B( K-1, 1 ),
! 460: $ LDB )
! 461: END IF
! 462: *
! 463: * Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1)
! 464: *
! 465: KP = -IPIV( K )
! 466: IF( KP.NE.K )
! 467: $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
! 468: *
! 469: KP = -IPIV( K-1 )
! 470: IF( KP.NE.K-1 )
! 471: $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
! 472: *
! 473: K = K - 2
! 474: END IF
! 475: *
! 476: GO TO 90
! 477: 100 CONTINUE
! 478: END IF
! 479: *
! 480: RETURN
! 481: *
! 482: * End of DSYTRS_ROOK
! 483: *
! 484: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytrs_rook.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack