1*> \brief \b SLAVSY_ROOK 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE SLAVSY_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B, 12* LDB, INFO ) 13* 14* .. Scalar Arguments .. 15* CHARACTER DIAG, TRANS, UPLO 16* INTEGER INFO, LDA, LDB, N, NRHS 17* .. 18* .. Array Arguments .. 19* INTEGER IPIV( * ) 20* REAL A( LDA, * ), B( LDB, * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> SLAVSY_ROOK performs one of the matrix-vector operations 30*> x := A*x or x := A'*x, 31*> where x is an N element vector and A is one of the factors 32*> from the block U*D*U' or L*D*L' factorization computed by SSYTRF_ROOK. 33*> 34*> If TRANS = 'N', multiplies by U or U * D (or L or L * D) 35*> If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L') 36*> If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L') 37*> \endverbatim 38* 39* Arguments: 40* ========== 41* 42*> \param[in] UPLO 43*> \verbatim 44*> UPLO is CHARACTER*1 45*> Specifies whether the factor stored in A is upper or lower 46*> triangular. 47*> = 'U': Upper triangular 48*> = 'L': Lower triangular 49*> \endverbatim 50*> 51*> \param[in] TRANS 52*> \verbatim 53*> TRANS is CHARACTER*1 54*> Specifies the operation to be performed: 55*> = 'N': x := A*x 56*> = 'T': x := A'*x 57*> = 'C': x := A'*x 58*> \endverbatim 59*> 60*> \param[in] DIAG 61*> \verbatim 62*> DIAG is CHARACTER*1 63*> Specifies whether or not the diagonal blocks are unit 64*> matrices. If the diagonal blocks are assumed to be unit, 65*> then A = U or A = L, otherwise A = U*D or A = L*D. 66*> = 'U': Diagonal blocks are assumed to be unit matrices. 67*> = 'N': Diagonal blocks are assumed to be non-unit matrices. 68*> \endverbatim 69*> 70*> \param[in] N 71*> \verbatim 72*> N is INTEGER 73*> The number of rows and columns of the matrix A. N >= 0. 74*> \endverbatim 75*> 76*> \param[in] NRHS 77*> \verbatim 78*> NRHS is INTEGER 79*> The number of right hand sides, i.e., the number of vectors 80*> x to be multiplied by A. NRHS >= 0. 81*> \endverbatim 82*> 83*> \param[in] A 84*> \verbatim 85*> A is REAL array, dimension (LDA,N) 86*> The block diagonal matrix D and the multipliers used to 87*> obtain the factor U or L as computed by SSYTRF_ROOK. 88*> Stored as a 2-D triangular matrix. 89*> \endverbatim 90*> 91*> \param[in] LDA 92*> \verbatim 93*> LDA is INTEGER 94*> The leading dimension of the array A. LDA >= max(1,N). 95*> \endverbatim 96*> 97*> \param[in] IPIV 98*> \verbatim 99*> IPIV is INTEGER array, dimension (N) 100*> Details of the interchanges and the block structure of D, 101*> as determined by SSYTRF_ROOK. 102*> 103*> If UPLO = 'U': 104*> If IPIV(k) > 0, then rows and columns k and IPIV(k) 105*> were interchanged and D(k,k) is a 1-by-1 diagonal block. 106*> (If IPIV( k ) = k, no interchange was done). 107*> 108*> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and 109*> columns k and -IPIV(k) were interchanged and rows and 110*> columns k-1 and -IPIV(k-1) were inerchaged, 111*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. 112*> 113*> If UPLO = 'L': 114*> If IPIV(k) > 0, then rows and columns k and IPIV(k) 115*> were interchanged and D(k,k) is a 1-by-1 diagonal block. 116*> (If IPIV( k ) = k, no interchange was done). 117*> 118*> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and 119*> columns k and -IPIV(k) were interchanged and rows and 120*> columns k+1 and -IPIV(k+1) were inerchaged, 121*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 122*> \endverbatim 123*> 124*> \param[in,out] B 125*> \verbatim 126*> B is REAL array, dimension (LDB,NRHS) 127*> On entry, B contains NRHS vectors of length N. 128*> On exit, B is overwritten with the product A * B. 129*> \endverbatim 130*> 131*> \param[in] LDB 132*> \verbatim 133*> LDB is INTEGER 134*> The leading dimension of the array B. LDB >= max(1,N). 135*> \endverbatim 136*> 137*> \param[out] INFO 138*> \verbatim 139*> INFO is INTEGER 140*> = 0: successful exit 141*> < 0: if INFO = -k, the k-th argument had an illegal value 142*> \endverbatim 143* 144* Authors: 145* ======== 146* 147*> \author Univ. of Tennessee 148*> \author Univ. of California Berkeley 149*> \author Univ. of Colorado Denver 150*> \author NAG Ltd. 151* 152*> \ingroup single_lin 153* 154* ===================================================================== 155 SUBROUTINE SLAVSY_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, 156 $ B, LDB, INFO ) 157* 158* -- LAPACK test routine -- 159* -- LAPACK is a software package provided by Univ. of Tennessee, -- 160* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 161* 162* .. Scalar Arguments .. 163 CHARACTER DIAG, TRANS, UPLO 164 INTEGER INFO, LDA, LDB, N, NRHS 165* .. 166* .. Array Arguments .. 167 INTEGER IPIV( * ) 168 REAL A( LDA, * ), B( LDB, * ) 169* .. 170* 171* ===================================================================== 172* 173* .. Parameters .. 174 REAL ONE 175 PARAMETER ( ONE = 1.0E+0 ) 176* .. 177* .. Local Scalars .. 178 LOGICAL NOUNIT 179 INTEGER J, K, KP 180 REAL D11, D12, D21, D22, T1, T2 181* .. 182* .. External Functions .. 183 LOGICAL LSAME 184 EXTERNAL LSAME 185* .. 186* .. External Subroutines .. 187 EXTERNAL SGEMV, SGER, SSCAL, SSWAP, XERBLA 188* .. 189* .. Intrinsic Functions .. 190 INTRINSIC ABS, MAX 191* .. 192* .. Executable Statements .. 193* 194* Test the input parameters. 195* 196 INFO = 0 197 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 198 INFO = -1 199 ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT. 200 $ LSAME( TRANS, 'T' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 201 INFO = -2 202 ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) ) 203 $ THEN 204 INFO = -3 205 ELSE IF( N.LT.0 ) THEN 206 INFO = -4 207 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 208 INFO = -6 209 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 210 INFO = -9 211 END IF 212 IF( INFO.NE.0 ) THEN 213 CALL XERBLA( 'SLAVSY_ROOK ', -INFO ) 214 RETURN 215 END IF 216* 217* Quick return if possible. 218* 219 IF( N.EQ.0 ) 220 $ RETURN 221* 222 NOUNIT = LSAME( DIAG, 'N' ) 223*------------------------------------------ 224* 225* Compute B := A * B (No transpose) 226* 227*------------------------------------------ 228 IF( LSAME( TRANS, 'N' ) ) THEN 229* 230* Compute B := U*B 231* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) 232* 233 IF( LSAME( UPLO, 'U' ) ) THEN 234* 235* Loop forward applying the transformations. 236* 237 K = 1 238 10 CONTINUE 239 IF( K.GT.N ) 240 $ GO TO 30 241 IF( IPIV( K ).GT.0 ) THEN 242* 243* 1 x 1 pivot block 244* 245* Multiply by the diagonal element if forming U * D. 246* 247 IF( NOUNIT ) 248 $ CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 249* 250* Multiply by P(K) * inv(U(K)) if K > 1. 251* 252 IF( K.GT.1 ) THEN 253* 254* Apply the transformation. 255* 256 CALL SGER( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ), 257 $ LDB, B( 1, 1 ), LDB ) 258* 259* Interchange if P(K) .ne. I. 260* 261 KP = IPIV( K ) 262 IF( KP.NE.K ) 263 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 264 END IF 265 K = K + 1 266 ELSE 267* 268* 2 x 2 pivot block 269* 270* Multiply by the diagonal block if forming U * D. 271* 272 IF( NOUNIT ) THEN 273 D11 = A( K, K ) 274 D22 = A( K+1, K+1 ) 275 D12 = A( K, K+1 ) 276 D21 = D12 277 DO 20 J = 1, NRHS 278 T1 = B( K, J ) 279 T2 = B( K+1, J ) 280 B( K, J ) = D11*T1 + D12*T2 281 B( K+1, J ) = D21*T1 + D22*T2 282 20 CONTINUE 283 END IF 284* 285* Multiply by P(K) * inv(U(K)) if K > 1. 286* 287 IF( K.GT.1 ) THEN 288* 289* Apply the transformations. 290* 291 CALL SGER( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ), 292 $ LDB, B( 1, 1 ), LDB ) 293 CALL SGER( K-1, NRHS, ONE, A( 1, K+1 ), 1, 294 $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB ) 295* 296* Interchange if a permutation was applied at the 297* K-th step of the factorization. 298* 299* Swap the first of pair with IMAXth 300* 301 KP = ABS( IPIV( K ) ) 302 IF( KP.NE.K ) 303 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 304* 305* NOW swap the first of pair with Pth 306* 307 KP = ABS( IPIV( K+1 ) ) 308 IF( KP.NE.K+1 ) 309 $ CALL SSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), 310 $ LDB ) 311 END IF 312 K = K + 2 313 END IF 314 GO TO 10 315 30 CONTINUE 316* 317* Compute B := L*B 318* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) . 319* 320 ELSE 321* 322* Loop backward applying the transformations to B. 323* 324 K = N 325 40 CONTINUE 326 IF( K.LT.1 ) 327 $ GO TO 60 328* 329* Test the pivot index. If greater than zero, a 1 x 1 330* pivot was used, otherwise a 2 x 2 pivot was used. 331* 332 IF( IPIV( K ).GT.0 ) THEN 333* 334* 1 x 1 pivot block: 335* 336* Multiply by the diagonal element if forming L * D. 337* 338 IF( NOUNIT ) 339 $ CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 340* 341* Multiply by P(K) * inv(L(K)) if K < N. 342* 343 IF( K.NE.N ) THEN 344 KP = IPIV( K ) 345* 346* Apply the transformation. 347* 348 CALL SGER( N-K, NRHS, ONE, A( K+1, K ), 1, B( K, 1 ), 349 $ LDB, B( K+1, 1 ), LDB ) 350* 351* Interchange if a permutation was applied at the 352* K-th step of the factorization. 353* 354 IF( KP.NE.K ) 355 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 356 END IF 357 K = K - 1 358* 359 ELSE 360* 361* 2 x 2 pivot block: 362* 363* Multiply by the diagonal block if forming L * D. 364* 365 IF( NOUNIT ) THEN 366 D11 = A( K-1, K-1 ) 367 D22 = A( K, K ) 368 D21 = A( K, K-1 ) 369 D12 = D21 370 DO 50 J = 1, NRHS 371 T1 = B( K-1, J ) 372 T2 = B( K, J ) 373 B( K-1, J ) = D11*T1 + D12*T2 374 B( K, J ) = D21*T1 + D22*T2 375 50 CONTINUE 376 END IF 377* 378* Multiply by P(K) * inv(L(K)) if K < N. 379* 380 IF( K.NE.N ) THEN 381* 382* Apply the transformation. 383* 384 CALL SGER( N-K, NRHS, ONE, A( K+1, K ), 1, B( K, 1 ), 385 $ LDB, B( K+1, 1 ), LDB ) 386 CALL SGER( N-K, NRHS, ONE, A( K+1, K-1 ), 1, 387 $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB ) 388* 389* Interchange if a permutation was applied at the 390* K-th step of the factorization. 391* 392* Swap the second of pair with IMAXth 393* 394 KP = ABS( IPIV( K ) ) 395 IF( KP.NE.K ) 396 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 397* 398* NOW swap the first of pair with Pth 399* 400 KP = ABS( IPIV( K-1 ) ) 401 IF( KP.NE.K-1 ) 402 $ CALL SSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), 403 $ LDB ) 404 END IF 405 K = K - 2 406 END IF 407 GO TO 40 408 60 CONTINUE 409 END IF 410*---------------------------------------- 411* 412* Compute B := A' * B (transpose) 413* 414*---------------------------------------- 415 ELSE 416* 417* Form B := U'*B 418* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) 419* and U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m) 420* 421 IF( LSAME( UPLO, 'U' ) ) THEN 422* 423* Loop backward applying the transformations. 424* 425 K = N 426 70 CONTINUE 427 IF( K.LT.1 ) 428 $ GO TO 90 429* 430* 1 x 1 pivot block. 431* 432 IF( IPIV( K ).GT.0 ) THEN 433 IF( K.GT.1 ) THEN 434* 435* Interchange if P(K) .ne. I. 436* 437 KP = IPIV( K ) 438 IF( KP.NE.K ) 439 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 440* 441* Apply the transformation 442* 443 CALL SGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB, 444 $ A( 1, K ), 1, ONE, B( K, 1 ), LDB ) 445 END IF 446 IF( NOUNIT ) 447 $ CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 448 K = K - 1 449* 450* 2 x 2 pivot block. 451* 452 ELSE 453 IF( K.GT.2 ) THEN 454* 455* Swap the second of pair with Pth 456* 457 KP = ABS( IPIV( K ) ) 458 IF( KP.NE.K ) 459 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 460* 461* Now swap the first of pair with IMAX(r)th 462* 463 KP = ABS( IPIV( K-1 ) ) 464 IF( KP.NE.K-1 ) 465 $ CALL SSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), 466 $ LDB ) 467* 468* Apply the transformations 469* 470 CALL SGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB, 471 $ A( 1, K ), 1, ONE, B( K, 1 ), LDB ) 472 CALL SGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB, 473 $ A( 1, K-1 ), 1, ONE, B( K-1, 1 ), LDB ) 474 END IF 475* 476* Multiply by the diagonal block if non-unit. 477* 478 IF( NOUNIT ) THEN 479 D11 = A( K-1, K-1 ) 480 D22 = A( K, K ) 481 D12 = A( K-1, K ) 482 D21 = D12 483 DO 80 J = 1, NRHS 484 T1 = B( K-1, J ) 485 T2 = B( K, J ) 486 B( K-1, J ) = D11*T1 + D12*T2 487 B( K, J ) = D21*T1 + D22*T2 488 80 CONTINUE 489 END IF 490 K = K - 2 491 END IF 492 GO TO 70 493 90 CONTINUE 494* 495* Form B := L'*B 496* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) 497* and L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1) 498* 499 ELSE 500* 501* Loop forward applying the L-transformations. 502* 503 K = 1 504 100 CONTINUE 505 IF( K.GT.N ) 506 $ GO TO 120 507* 508* 1 x 1 pivot block 509* 510 IF( IPIV( K ).GT.0 ) THEN 511 IF( K.LT.N ) THEN 512* 513* Interchange if P(K) .ne. I. 514* 515 KP = IPIV( K ) 516 IF( KP.NE.K ) 517 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 518* 519* Apply the transformation 520* 521 CALL SGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ), 522 $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB ) 523 END IF 524 IF( NOUNIT ) 525 $ CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 526 K = K + 1 527* 528* 2 x 2 pivot block. 529* 530 ELSE 531 IF( K.LT.N-1 ) THEN 532* 533* Swap the first of pair with Pth 534* 535 KP = ABS( IPIV( K ) ) 536 IF( KP.NE.K ) 537 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 538* 539* Now swap the second of pair with IMAX(r)th 540* 541 KP = ABS( IPIV( K+1 ) ) 542 IF( KP.NE.K+1 ) 543 $ CALL SSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), 544 $ LDB ) 545* 546* Apply the transformation 547* 548 CALL SGEMV( 'Transpose', N-K-1, NRHS, ONE, 549 $ B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, ONE, 550 $ B( K+1, 1 ), LDB ) 551 CALL SGEMV( 'Transpose', N-K-1, NRHS, ONE, 552 $ B( K+2, 1 ), LDB, A( K+2, K ), 1, ONE, 553 $ B( K, 1 ), LDB ) 554 END IF 555* 556* Multiply by the diagonal block if non-unit. 557* 558 IF( NOUNIT ) THEN 559 D11 = A( K, K ) 560 D22 = A( K+1, K+1 ) 561 D21 = A( K+1, K ) 562 D12 = D21 563 DO 110 J = 1, NRHS 564 T1 = B( K, J ) 565 T2 = B( K+1, J ) 566 B( K, J ) = D11*T1 + D12*T2 567 B( K+1, J ) = D21*T1 + D22*T2 568 110 CONTINUE 569 END IF 570 K = K + 2 571 END IF 572 GO TO 100 573 120 CONTINUE 574 END IF 575* 576 END IF 577 RETURN 578* 579* End of SLAVSY_ROOK 580* 581 END 582