TESTING/LIN/clavhp.f

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