LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zlavhp.f
Go to the documentation of this file.
00001 *> \brief \b ZLAVHP
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE ZLAVHP( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
00012 *                          INFO )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          DIAG, TRANS, UPLO
00016 *       INTEGER            INFO, LDB, N, NRHS
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       INTEGER            IPIV( * )
00020 *       COMPLEX*16         A( * ), B( LDB, * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *>    ZLAVHP  performs one of the matrix-vector operations
00030 *>       x := A*x  or  x := A^H*x,
00031 *>    where x is an N element vector and  A is one of the factors
00032 *>    from the symmetric factorization computed by ZHPTRF.
00033 *>    ZHPTRF produces a factorization of the form
00034 *>         U * D * U^H     or     L * D * L^H,
00035 *>    where U (or L) is a product of permutation and unit upper (lower)
00036 *>    triangular matrices, U^H (or L^H) is the conjugate transpose of
00037 *>    U (or L), and D is Hermitian and block diagonal with 1 x 1 and
00038 *>    2 x 2 diagonal blocks.  The multipliers for the transformations
00039 *>    and the upper or lower triangular parts of the diagonal blocks
00040 *>    are stored columnwise in packed format in the linear array A.
00041 *>
00042 *>    If TRANS = 'N' or 'n', ZLAVHP multiplies either by U or U * D
00043 *>    (or L or L * D).
00044 *>    If TRANS = 'C' or 'c', ZLAVHP multiplies either by U^H or D * U^H
00045 *>    (or L^H or D * L^H ).
00046 *> \endverbatim
00047 *
00048 *  Arguments:
00049 *  ==========
00050 *
00051 *> \verbatim
00052 *>  UPLO   - CHARACTER*1
00053 *>           On entry, UPLO specifies whether the triangular matrix
00054 *>           stored in A is upper or lower triangular.
00055 *>              UPLO = 'U' or 'u'   The matrix is upper triangular.
00056 *>              UPLO = 'L' or 'l'   The matrix is lower triangular.
00057 *>           Unchanged on exit.
00058 *>
00059 *>  TRANS  - CHARACTER*1
00060 *>           On entry, TRANS specifies the operation to be performed as
00061 *>           follows:
00062 *>              TRANS = 'N' or 'n'   x := A*x.
00063 *>              TRANS = 'C' or 'c'   x := A^H*x.
00064 *>           Unchanged on exit.
00065 *>
00066 *>  DIAG   - CHARACTER*1
00067 *>           On entry, DIAG specifies whether the diagonal blocks are
00068 *>           assumed to be unit matrices, as follows:
00069 *>              DIAG = 'U' or 'u'   Diagonal blocks are unit matrices.
00070 *>              DIAG = 'N' or 'n'   Diagonal blocks are non-unit.
00071 *>           Unchanged on exit.
00072 *>
00073 *>  N      - INTEGER
00074 *>           On entry, N specifies the order of the matrix A.
00075 *>           N must be at least zero.
00076 *>           Unchanged on exit.
00077 *>
00078 *>  NRHS   - INTEGER
00079 *>           On entry, NRHS specifies the number of right hand sides,
00080 *>           i.e., the number of vectors x to be multiplied by A.
00081 *>           NRHS must be at least zero.
00082 *>           Unchanged on exit.
00083 *>
00084 *>  A      - COMPLEX*16 array, dimension( N*(N+1)/2 )
00085 *>           On entry, A contains a block diagonal matrix and the
00086 *>           multipliers of the transformations used to obtain it,
00087 *>           stored as a packed triangular matrix.
00088 *>           Unchanged on exit.
00089 *>
00090 *>  IPIV   - INTEGER array, dimension( N )
00091 *>           On entry, IPIV contains the vector of pivot indices as
00092 *>           determined by ZSPTRF or ZHPTRF.
00093 *>           If IPIV( K ) = K, no interchange was done.
00094 *>           If IPIV( K ) <> K but IPIV( K ) > 0, then row K was inter-
00095 *>           changed with row IPIV( K ) and a 1 x 1 pivot block was used.
00096 *>           If IPIV( K ) < 0 and UPLO = 'U', then row K-1 was exchanged
00097 *>           with row | IPIV( K ) | and a 2 x 2 pivot block was used.
00098 *>           If IPIV( K ) < 0 and UPLO = 'L', then row K+1 was exchanged
00099 *>           with row | IPIV( K ) | and a 2 x 2 pivot block was used.
00100 *>
00101 *>  B      - COMPLEX*16 array, dimension( LDB, NRHS )
00102 *>           On entry, B contains NRHS vectors of length N.
00103 *>           On exit, B is overwritten with the product A * B.
00104 *>
00105 *>  LDB    - INTEGER
00106 *>           On entry, LDB contains the leading dimension of B as
00107 *>           declared in the calling program.  LDB must be at least
00108 *>           max( 1, N ).
00109 *>           Unchanged on exit.
00110 *>
00111 *>  INFO   - INTEGER
00112 *>           INFO is the error flag.
00113 *>           On exit, a value of 0 indicates a successful exit.
00114 *>           A negative value, say -K, indicates that the K-th argument
00115 *>           has an illegal value.
00116 *> \endverbatim
00117 *
00118 *  Authors:
00119 *  ========
00120 *
00121 *> \author Univ. of Tennessee 
00122 *> \author Univ. of California Berkeley 
00123 *> \author Univ. of Colorado Denver 
00124 *> \author NAG Ltd. 
00125 *
00126 *> \date November 2011
00127 *
00128 *> \ingroup complex16_lin
00129 *
00130 *  =====================================================================
00131       SUBROUTINE ZLAVHP( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
00132      $                   INFO )
00133 *
00134 *  -- LAPACK test routine (version 3.4.0) --
00135 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00136 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00137 *     November 2011
00138 *
00139 *     .. Scalar Arguments ..
00140       CHARACTER          DIAG, TRANS, UPLO
00141       INTEGER            INFO, LDB, N, NRHS
00142 *     ..
00143 *     .. Array Arguments ..
00144       INTEGER            IPIV( * )
00145       COMPLEX*16         A( * ), B( LDB, * )
00146 *     ..
00147 *
00148 *  =====================================================================
00149 *
00150 *     .. Parameters ..
00151       COMPLEX*16         ONE
00152       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
00153 *     ..
00154 *     .. Local Scalars ..
00155       LOGICAL            NOUNIT
00156       INTEGER            J, K, KC, KCNEXT, KP
00157       COMPLEX*16         D11, D12, D21, D22, T1, T2
00158 *     ..
00159 *     .. External Functions ..
00160       LOGICAL            LSAME
00161       EXTERNAL           LSAME
00162 *     ..
00163 *     .. External Subroutines ..
00164       EXTERNAL           XERBLA, ZGEMV, ZGERU, ZLACGV, ZSCAL, ZSWAP
00165 *     ..
00166 *     .. Intrinsic Functions ..
00167       INTRINSIC          ABS, DCONJG, MAX
00168 *     ..
00169 *     .. Executable Statements ..
00170 *
00171 *     Test the input parameters.
00172 *
00173       INFO = 0
00174       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00175          INFO = -1
00176       ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'C' ) )
00177      $          THEN
00178          INFO = -2
00179       ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
00180      $          THEN
00181          INFO = -3
00182       ELSE IF( N.LT.0 ) THEN
00183          INFO = -4
00184       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00185          INFO = -8
00186       END IF
00187       IF( INFO.NE.0 ) THEN
00188          CALL XERBLA( 'ZLAVHP ', -INFO )
00189          RETURN
00190       END IF
00191 *
00192 *     Quick return if possible.
00193 *
00194       IF( N.EQ.0 )
00195      $   RETURN
00196 *
00197       NOUNIT = LSAME( DIAG, 'N' )
00198 *------------------------------------------
00199 *
00200 *     Compute  B := A * B  (No transpose)
00201 *
00202 *------------------------------------------
00203       IF( LSAME( TRANS, 'N' ) ) THEN
00204 *
00205 *        Compute  B := U*B
00206 *        where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
00207 *
00208          IF( LSAME( UPLO, 'U' ) ) THEN
00209 *
00210 *        Loop forward applying the transformations.
00211 *
00212             K = 1
00213             KC = 1
00214    10       CONTINUE
00215             IF( K.GT.N )
00216      $         GO TO 30
00217 *
00218 *           1 x 1 pivot block
00219 *
00220             IF( IPIV( K ).GT.0 ) THEN
00221 *
00222 *              Multiply by the diagonal element if forming U * D.
00223 *
00224                IF( NOUNIT )
00225      $            CALL ZSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
00226 *
00227 *              Multiply by P(K) * inv(U(K))  if K > 1.
00228 *
00229                IF( K.GT.1 ) THEN
00230 *
00231 *                 Apply the transformation.
00232 *
00233                   CALL ZGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
00234      $                        LDB, B( 1, 1 ), LDB )
00235 *
00236 *                 Interchange if P(K) != I.
00237 *
00238                   KP = IPIV( K )
00239                   IF( KP.NE.K )
00240      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00241                END IF
00242                KC = KC + K
00243                K = K + 1
00244             ELSE
00245 *
00246 *              2 x 2 pivot block
00247 *
00248                KCNEXT = KC + K
00249 *
00250 *              Multiply by the diagonal block if forming U * D.
00251 *
00252                IF( NOUNIT ) THEN
00253                   D11 = A( KCNEXT-1 )
00254                   D22 = A( KCNEXT+K )
00255                   D12 = A( KCNEXT+K-1 )
00256                   D21 = DCONJG( D12 )
00257                   DO 20 J = 1, NRHS
00258                      T1 = B( K, J )
00259                      T2 = B( K+1, J )
00260                      B( K, J ) = D11*T1 + D12*T2
00261                      B( K+1, J ) = D21*T1 + D22*T2
00262    20             CONTINUE
00263                END IF
00264 *
00265 *              Multiply by  P(K) * inv(U(K))  if K > 1.
00266 *
00267                IF( K.GT.1 ) THEN
00268 *
00269 *                 Apply the transformations.
00270 *
00271                   CALL ZGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
00272      $                        LDB, B( 1, 1 ), LDB )
00273                   CALL ZGERU( K-1, NRHS, ONE, A( KCNEXT ), 1,
00274      $                        B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
00275 *
00276 *                 Interchange if P(K) != I.
00277 *
00278                   KP = ABS( IPIV( K ) )
00279                   IF( KP.NE.K )
00280      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00281                END IF
00282                KC = KCNEXT + K + 1
00283                K = K + 2
00284             END IF
00285             GO TO 10
00286    30       CONTINUE
00287 *
00288 *        Compute  B := L*B
00289 *        where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
00290 *
00291          ELSE
00292 *
00293 *           Loop backward applying the transformations to B.
00294 *
00295             K = N
00296             KC = N*( N+1 ) / 2 + 1
00297    40       CONTINUE
00298             IF( K.LT.1 )
00299      $         GO TO 60
00300             KC = KC - ( N-K+1 )
00301 *
00302 *           Test the pivot index.  If greater than zero, a 1 x 1
00303 *           pivot was used, otherwise a 2 x 2 pivot was used.
00304 *
00305             IF( IPIV( K ).GT.0 ) THEN
00306 *
00307 *              1 x 1 pivot block:
00308 *
00309 *              Multiply by the diagonal element if forming L * D.
00310 *
00311                IF( NOUNIT )
00312      $            CALL ZSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
00313 *
00314 *              Multiply by  P(K) * inv(L(K))  if K < N.
00315 *
00316                IF( K.NE.N ) THEN
00317                   KP = IPIV( K )
00318 *
00319 *                 Apply the transformation.
00320 *
00321                   CALL ZGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
00322      $                        LDB, B( K+1, 1 ), LDB )
00323 *
00324 *                 Interchange if a permutation was applied at the
00325 *                 K-th step of the factorization.
00326 *
00327                   IF( KP.NE.K )
00328      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00329                END IF
00330                K = K - 1
00331 *
00332             ELSE
00333 *
00334 *              2 x 2 pivot block:
00335 *
00336                KCNEXT = KC - ( N-K+2 )
00337 *
00338 *              Multiply by the diagonal block if forming L * D.
00339 *
00340                IF( NOUNIT ) THEN
00341                   D11 = A( KCNEXT )
00342                   D22 = A( KC )
00343                   D21 = A( KCNEXT+1 )
00344                   D12 = DCONJG( D21 )
00345                   DO 50 J = 1, NRHS
00346                      T1 = B( K-1, J )
00347                      T2 = B( K, J )
00348                      B( K-1, J ) = D11*T1 + D12*T2
00349                      B( K, J ) = D21*T1 + D22*T2
00350    50             CONTINUE
00351                END IF
00352 *
00353 *              Multiply by  P(K) * inv(L(K))  if K < N.
00354 *
00355                IF( K.NE.N ) THEN
00356 *
00357 *                 Apply the transformation.
00358 *
00359                   CALL ZGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
00360      $                        LDB, B( K+1, 1 ), LDB )
00361                   CALL ZGERU( N-K, NRHS, ONE, A( KCNEXT+2 ), 1,
00362      $                        B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
00363 *
00364 *                 Interchange if a permutation was applied at the
00365 *                 K-th step of the factorization.
00366 *
00367                   KP = ABS( IPIV( K ) )
00368                   IF( KP.NE.K )
00369      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00370                END IF
00371                KC = KCNEXT
00372                K = K - 2
00373             END IF
00374             GO TO 40
00375    60       CONTINUE
00376          END IF
00377 *-------------------------------------------------
00378 *
00379 *     Compute  B := A^H * B  (conjugate transpose)
00380 *
00381 *-------------------------------------------------
00382       ELSE
00383 *
00384 *        Form  B := U^H*B
00385 *        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
00386 *        and   U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
00387 *
00388          IF( LSAME( UPLO, 'U' ) ) THEN
00389 *
00390 *           Loop backward applying the transformations.
00391 *
00392             K = N
00393             KC = N*( N+1 ) / 2 + 1
00394    70       CONTINUE
00395             IF( K.LT.1 )
00396      $         GO TO 90
00397             KC = KC - K
00398 *
00399 *           1 x 1 pivot block.
00400 *
00401             IF( IPIV( K ).GT.0 ) THEN
00402                IF( K.GT.1 ) THEN
00403 *
00404 *                 Interchange if P(K) != I.
00405 *
00406                   KP = IPIV( K )
00407                   IF( KP.NE.K )
00408      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00409 *
00410 *                 Apply the transformation:
00411 *                    y := y - B' * conjg(x)
00412 *                 where x is a column of A and y is a row of B.
00413 *
00414                   CALL ZLACGV( NRHS, B( K, 1 ), LDB )
00415                   CALL ZGEMV( 'Conjugate', K-1, NRHS, ONE, B, LDB,
00416      $                        A( KC ), 1, ONE, B( K, 1 ), LDB )
00417                   CALL ZLACGV( NRHS, B( K, 1 ), LDB )
00418                END IF
00419                IF( NOUNIT )
00420      $            CALL ZSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
00421                K = K - 1
00422 *
00423 *           2 x 2 pivot block.
00424 *
00425             ELSE
00426                KCNEXT = KC - ( K-1 )
00427                IF( K.GT.2 ) THEN
00428 *
00429 *                 Interchange if P(K) != I.
00430 *
00431                   KP = ABS( IPIV( K ) )
00432                   IF( KP.NE.K-1 )
00433      $               CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
00434      $                           LDB )
00435 *
00436 *                 Apply the transformations.
00437 *
00438                   CALL ZLACGV( NRHS, B( K, 1 ), LDB )
00439                   CALL ZGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
00440      $                        A( KC ), 1, ONE, B( K, 1 ), LDB )
00441                   CALL ZLACGV( NRHS, B( K, 1 ), LDB )
00442 *
00443                   CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
00444                   CALL ZGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
00445      $                        A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
00446                   CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
00447                END IF
00448 *
00449 *              Multiply by the diagonal block if non-unit.
00450 *
00451                IF( NOUNIT ) THEN
00452                   D11 = A( KC-1 )
00453                   D22 = A( KC+K-1 )
00454                   D12 = A( KC+K-2 )
00455                   D21 = DCONJG( D12 )
00456                   DO 80 J = 1, NRHS
00457                      T1 = B( K-1, J )
00458                      T2 = B( K, J )
00459                      B( K-1, J ) = D11*T1 + D12*T2
00460                      B( K, J ) = D21*T1 + D22*T2
00461    80             CONTINUE
00462                END IF
00463                KC = KCNEXT
00464                K = K - 2
00465             END IF
00466             GO TO 70
00467    90       CONTINUE
00468 *
00469 *        Form  B := L^H*B
00470 *        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
00471 *        and   L^H = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
00472 *
00473          ELSE
00474 *
00475 *           Loop forward applying the L-transformations.
00476 *
00477             K = 1
00478             KC = 1
00479   100       CONTINUE
00480             IF( K.GT.N )
00481      $         GO TO 120
00482 *
00483 *           1 x 1 pivot block
00484 *
00485             IF( IPIV( K ).GT.0 ) THEN
00486                IF( K.LT.N ) THEN
00487 *
00488 *                 Interchange if P(K) != I.
00489 *
00490                   KP = IPIV( K )
00491                   IF( KP.NE.K )
00492      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00493 *
00494 *                 Apply the transformation
00495 *
00496                   CALL ZLACGV( NRHS, B( K, 1 ), LDB )
00497                   CALL ZGEMV( 'Conjugate', N-K, NRHS, ONE, B( K+1, 1 ),
00498      $                        LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
00499                   CALL ZLACGV( NRHS, B( K, 1 ), LDB )
00500                END IF
00501                IF( NOUNIT )
00502      $            CALL ZSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
00503                KC = KC + N - K + 1
00504                K = K + 1
00505 *
00506 *           2 x 2 pivot block.
00507 *
00508             ELSE
00509                KCNEXT = KC + N - K + 1
00510                IF( K.LT.N-1 ) THEN
00511 *
00512 *              Interchange if P(K) != I.
00513 *
00514                   KP = ABS( IPIV( K ) )
00515                   IF( KP.NE.K+1 )
00516      $               CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
00517      $                           LDB )
00518 *
00519 *                 Apply the transformation
00520 *
00521                   CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
00522                   CALL ZGEMV( 'Conjugate', N-K-1, NRHS, ONE,
00523      $                        B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
00524      $                        B( K+1, 1 ), LDB )
00525                   CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
00526 *
00527                   CALL ZLACGV( NRHS, B( K, 1 ), LDB )
00528                   CALL ZGEMV( 'Conjugate', N-K-1, NRHS, ONE,
00529      $                        B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
00530      $                        B( K, 1 ), LDB )
00531                   CALL ZLACGV( NRHS, B( K, 1 ), LDB )
00532                END IF
00533 *
00534 *              Multiply by the diagonal block if non-unit.
00535 *
00536                IF( NOUNIT ) THEN
00537                   D11 = A( KC )
00538                   D22 = A( KCNEXT )
00539                   D21 = A( KC+1 )
00540                   D12 = DCONJG( D21 )
00541                   DO 110 J = 1, NRHS
00542                      T1 = B( K, J )
00543                      T2 = B( K+1, J )
00544                      B( K, J ) = D11*T1 + D12*T2
00545                      B( K+1, J ) = D21*T1 + D22*T2
00546   110             CONTINUE
00547                END IF
00548                KC = KCNEXT + ( N-K )
00549                K = K + 2
00550             END IF
00551             GO TO 100
00552   120       CONTINUE
00553          END IF
00554 *
00555       END IF
00556       RETURN
00557 *
00558 *     End of ZLAVHP
00559 *
00560       END
 All Files Functions