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