1*> \brief \b CLAVHE 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 CLAVHE( 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 A( LDA, * ), B( LDB, * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> CLAVHE performs one of the matrix-vector operations 30*> x := A*x or x := A^H*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 CHETRF. 33*> 34*> If TRANS = 'N', multiplies by U or U * D (or L or L * D) 35*> If TRANS = 'C', 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*> = 'C': x := A^H*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 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 CHETRF_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 CHETRF. 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) = IPIV(k-1) < 0, then rows and 107*> columns k-1 and -IPIV(k) were interchanged, 108*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. 109*> 110*> If UPLO = 'L': 111*> If IPIV(k) > 0, then rows and columns k and IPIV(k) 112*> were interchanged and D(k,k) is a 1-by-1 diagonal block. 113*> (If IPIV( k ) = k, no interchange was done). 114*> 115*> If IPIV(k) = IPIV(k+1) < 0, then rows and 116*> columns k+1 and -IPIV(k) were interchanged, 117*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 118*> \endverbatim 119*> 120*> \param[in,out] B 121*> \verbatim 122*> B is COMPLEX array, dimension (LDB,NRHS) 123*> On entry, B contains NRHS vectors of length N. 124*> On exit, B is overwritten with the product A * B. 125*> \endverbatim 126*> 127*> \param[in] LDB 128*> \verbatim 129*> LDB is INTEGER 130*> The leading dimension of the array B. LDB >= max(1,N). 131*> \endverbatim 132*> 133*> \param[out] INFO 134*> \verbatim 135*> INFO is INTEGER 136*> = 0: successful exit 137*> < 0: if INFO = -k, the k-th argument had an illegal value 138*> \endverbatim 139* 140* Authors: 141* ======== 142* 143*> \author Univ. of Tennessee 144*> \author Univ. of California Berkeley 145*> \author Univ. of Colorado Denver 146*> \author NAG Ltd. 147* 148*> \ingroup complex_lin 149* 150* ===================================================================== 151 SUBROUTINE CLAVHE( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B, 152 $ LDB, INFO ) 153* 154* -- LAPACK test routine -- 155* -- LAPACK is a software package provided by Univ. of Tennessee, -- 156* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 157* 158* .. Scalar Arguments .. 159 CHARACTER DIAG, TRANS, UPLO 160 INTEGER INFO, LDA, LDB, N, NRHS 161* .. 162* .. Array Arguments .. 163 INTEGER IPIV( * ) 164 COMPLEX A( LDA, * ), B( LDB, * ) 165* .. 166* 167* ===================================================================== 168* 169* .. Parameters .. 170 COMPLEX ONE 171 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) 172* .. 173* .. Local Scalars .. 174 LOGICAL NOUNIT 175 INTEGER J, K, KP 176 COMPLEX D11, D12, D21, D22, T1, T2 177* .. 178* .. External Functions .. 179 LOGICAL LSAME 180 EXTERNAL LSAME 181* .. 182* .. External Subroutines .. 183 EXTERNAL CGEMV, CGERU, CLACGV, CSCAL, CSWAP, XERBLA 184* .. 185* .. Intrinsic Functions .. 186 INTRINSIC ABS, CONJG, MAX 187* .. 188* .. Executable Statements .. 189* 190* Test the input parameters. 191* 192 INFO = 0 193 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 194 INFO = -1 195 ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) 196 $ THEN 197 INFO = -2 198 ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) ) 199 $ THEN 200 INFO = -3 201 ELSE IF( N.LT.0 ) THEN 202 INFO = -4 203 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 204 INFO = -6 205 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 206 INFO = -9 207 END IF 208 IF( INFO.NE.0 ) THEN 209 CALL XERBLA( 'CLAVHE ', -INFO ) 210 RETURN 211 END IF 212* 213* Quick return if possible. 214* 215 IF( N.EQ.0 ) 216 $ RETURN 217* 218 NOUNIT = LSAME( DIAG, 'N' ) 219*------------------------------------------ 220* 221* Compute B := A * B (No transpose) 222* 223*------------------------------------------ 224 IF( LSAME( TRANS, 'N' ) ) THEN 225* 226* Compute B := U*B 227* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) 228* 229 IF( LSAME( UPLO, 'U' ) ) THEN 230* 231* Loop forward applying the transformations. 232* 233 K = 1 234 10 CONTINUE 235 IF( K.GT.N ) 236 $ GO TO 30 237 IF( IPIV( K ).GT.0 ) THEN 238* 239* 1 x 1 pivot block 240* 241* Multiply by the diagonal element if forming U * D. 242* 243 IF( NOUNIT ) 244 $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 245* 246* Multiply by P(K) * inv(U(K)) if K > 1. 247* 248 IF( K.GT.1 ) THEN 249* 250* Apply the transformation. 251* 252 CALL CGERU( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ), 253 $ LDB, B( 1, 1 ), LDB ) 254* 255* Interchange if P(K) != I. 256* 257 KP = IPIV( K ) 258 IF( KP.NE.K ) 259 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 260 END IF 261 K = K + 1 262 ELSE 263* 264* 2 x 2 pivot block 265* 266* Multiply by the diagonal block if forming U * D. 267* 268 IF( NOUNIT ) THEN 269 D11 = A( K, K ) 270 D22 = A( K+1, K+1 ) 271 D12 = A( K, K+1 ) 272 D21 = CONJG( D12 ) 273 DO 20 J = 1, NRHS 274 T1 = B( K, J ) 275 T2 = B( K+1, J ) 276 B( K, J ) = D11*T1 + D12*T2 277 B( K+1, J ) = D21*T1 + D22*T2 278 20 CONTINUE 279 END IF 280* 281* Multiply by P(K) * inv(U(K)) if K > 1. 282* 283 IF( K.GT.1 ) THEN 284* 285* Apply the transformations. 286* 287 CALL CGERU( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ), 288 $ LDB, B( 1, 1 ), LDB ) 289 CALL CGERU( K-1, NRHS, ONE, A( 1, K+1 ), 1, 290 $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB ) 291* 292* Interchange if P(K) != I. 293* 294 KP = ABS( IPIV( K ) ) 295 IF( KP.NE.K ) 296 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 297 END IF 298 K = K + 2 299 END IF 300 GO TO 10 301 30 CONTINUE 302* 303* Compute B := L*B 304* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) . 305* 306 ELSE 307* 308* Loop backward applying the transformations to B. 309* 310 K = N 311 40 CONTINUE 312 IF( K.LT.1 ) 313 $ GO TO 60 314* 315* Test the pivot index. If greater than zero, a 1 x 1 316* pivot was used, otherwise a 2 x 2 pivot was used. 317* 318 IF( IPIV( K ).GT.0 ) THEN 319* 320* 1 x 1 pivot block: 321* 322* Multiply by the diagonal element if forming L * D. 323* 324 IF( NOUNIT ) 325 $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 326* 327* Multiply by P(K) * inv(L(K)) if K < N. 328* 329 IF( K.NE.N ) THEN 330 KP = IPIV( K ) 331* 332* Apply the transformation. 333* 334 CALL CGERU( N-K, NRHS, ONE, A( K+1, K ), 1, 335 $ B( K, 1 ), LDB, B( K+1, 1 ), LDB ) 336* 337* Interchange if a permutation was applied at the 338* K-th step of the factorization. 339* 340 IF( KP.NE.K ) 341 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 342 END IF 343 K = K - 1 344* 345 ELSE 346* 347* 2 x 2 pivot block: 348* 349* Multiply by the diagonal block if forming L * D. 350* 351 IF( NOUNIT ) THEN 352 D11 = A( K-1, K-1 ) 353 D22 = A( K, K ) 354 D21 = A( K, K-1 ) 355 D12 = CONJG( D21 ) 356 DO 50 J = 1, NRHS 357 T1 = B( K-1, J ) 358 T2 = B( K, J ) 359 B( K-1, J ) = D11*T1 + D12*T2 360 B( K, J ) = D21*T1 + D22*T2 361 50 CONTINUE 362 END IF 363* 364* Multiply by P(K) * inv(L(K)) if K < N. 365* 366 IF( K.NE.N ) THEN 367* 368* Apply the transformation. 369* 370 CALL CGERU( N-K, NRHS, ONE, A( K+1, K ), 1, 371 $ B( K, 1 ), LDB, B( K+1, 1 ), LDB ) 372 CALL CGERU( N-K, NRHS, ONE, A( K+1, K-1 ), 1, 373 $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB ) 374* 375* Interchange if a permutation was applied at the 376* K-th step of the factorization. 377* 378 KP = ABS( IPIV( K ) ) 379 IF( KP.NE.K ) 380 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 381 END IF 382 K = K - 2 383 END IF 384 GO TO 40 385 60 CONTINUE 386 END IF 387*-------------------------------------------------- 388* 389* Compute B := A^H * B (conjugate transpose) 390* 391*-------------------------------------------------- 392 ELSE 393* 394* Form B := U^H*B 395* where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) 396* and U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m) 397* 398 IF( LSAME( UPLO, 'U' ) ) THEN 399* 400* Loop backward applying the transformations. 401* 402 K = N 403 70 IF( K.LT.1 ) 404 $ GO TO 90 405* 406* 1 x 1 pivot block. 407* 408 IF( IPIV( K ).GT.0 ) THEN 409 IF( K.GT.1 ) THEN 410* 411* Interchange if P(K) != I. 412* 413 KP = IPIV( K ) 414 IF( KP.NE.K ) 415 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 416* 417* Apply the transformation 418* y = y - B' conjg(x), 419* where x is a column of A and y is a row of B. 420* 421 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 422 CALL CGEMV( 'Conjugate', K-1, NRHS, ONE, B, LDB, 423 $ A( 1, K ), 1, ONE, B( K, 1 ), LDB ) 424 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 425 END IF 426 IF( NOUNIT ) 427 $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 428 K = K - 1 429* 430* 2 x 2 pivot block. 431* 432 ELSE 433 IF( K.GT.2 ) THEN 434* 435* Interchange if P(K) != I. 436* 437 KP = ABS( IPIV( K ) ) 438 IF( KP.NE.K-1 ) 439 $ CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), 440 $ LDB ) 441* 442* Apply the transformations 443* y = y - B' conjg(x), 444* where x is a block column of A and y is a block 445* row of B. 446* 447 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 448 CALL CGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB, 449 $ A( 1, K ), 1, ONE, B( K, 1 ), LDB ) 450 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 451* 452 CALL CLACGV( NRHS, B( K-1, 1 ), LDB ) 453 CALL CGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB, 454 $ A( 1, K-1 ), 1, ONE, B( K-1, 1 ), LDB ) 455 CALL CLACGV( NRHS, B( K-1, 1 ), LDB ) 456 END IF 457* 458* Multiply by the diagonal block if non-unit. 459* 460 IF( NOUNIT ) THEN 461 D11 = A( K-1, K-1 ) 462 D22 = A( K, K ) 463 D12 = A( K-1, K ) 464 D21 = CONJG( D12 ) 465 DO 80 J = 1, NRHS 466 T1 = B( K-1, J ) 467 T2 = B( K, J ) 468 B( K-1, J ) = D11*T1 + D12*T2 469 B( K, J ) = D21*T1 + D22*T2 470 80 CONTINUE 471 END IF 472 K = K - 2 473 END IF 474 GO TO 70 475 90 CONTINUE 476* 477* Form B := L^H*B 478* where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) 479* and L^H = inv(L^H(m))*P(m)* ... *inv(L^H(1))*P(1) 480* 481 ELSE 482* 483* Loop forward applying the L-transformations. 484* 485 K = 1 486 100 CONTINUE 487 IF( K.GT.N ) 488 $ GO TO 120 489* 490* 1 x 1 pivot block 491* 492 IF( IPIV( K ).GT.0 ) THEN 493 IF( K.LT.N ) THEN 494* 495* Interchange if P(K) != I. 496* 497 KP = IPIV( K ) 498 IF( KP.NE.K ) 499 $ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 500* 501* Apply the transformation 502* 503 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 504 CALL CGEMV( 'Conjugate', N-K, NRHS, ONE, B( K+1, 1 ), 505 $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB ) 506 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 507 END IF 508 IF( NOUNIT ) 509 $ CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 510 K = K + 1 511* 512* 2 x 2 pivot block. 513* 514 ELSE 515 IF( K.LT.N-1 ) THEN 516* 517* Interchange if P(K) != I. 518* 519 KP = ABS( IPIV( K ) ) 520 IF( KP.NE.K+1 ) 521 $ CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), 522 $ LDB ) 523* 524* Apply the transformation 525* 526 CALL CLACGV( NRHS, B( K+1, 1 ), LDB ) 527 CALL CGEMV( 'Conjugate', N-K-1, NRHS, ONE, 528 $ B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, ONE, 529 $ B( K+1, 1 ), LDB ) 530 CALL CLACGV( NRHS, B( K+1, 1 ), LDB ) 531* 532 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 533 CALL CGEMV( 'Conjugate', N-K-1, NRHS, ONE, 534 $ B( K+2, 1 ), LDB, A( K+2, K ), 1, ONE, 535 $ B( K, 1 ), LDB ) 536 CALL CLACGV( NRHS, B( K, 1 ), LDB ) 537 END IF 538* 539* Multiply by the diagonal block if non-unit. 540* 541 IF( NOUNIT ) THEN 542 D11 = A( K, K ) 543 D22 = A( K+1, K+1 ) 544 D21 = A( K+1, K ) 545 D12 = CONJG( D21 ) 546 DO 110 J = 1, NRHS 547 T1 = B( K, J ) 548 T2 = B( K+1, J ) 549 B( K, J ) = D11*T1 + D12*T2 550 B( K+1, J ) = D21*T1 + D22*T2 551 110 CONTINUE 552 END IF 553 K = K + 2 554 END IF 555 GO TO 100 556 120 CONTINUE 557 END IF 558* 559 END IF 560 RETURN 561* 562* End of CLAVHE 563* 564 END 565