LAPACK: TESTING/LIN/zlavhe.f Source File

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