LAPACK: TESTING/LIN/zlavhp.f Source File

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