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