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