LAPACK: TESTING/LIN/clavhe.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLAVHE( 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            A( LDA, * ), B( LDB, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *     CLAVHE  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 CHETRF.
00024 *     CHETRF 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', CLAVHE multiplies either by U or U * D
00035 *     (or L or L * D).
00036 *     If TRANS = 'C' or 'c', CLAVHE 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 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 CSYTRF or CHETRF.
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 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            ONE
00117       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
00118 *     ..
00119 *     .. Local Scalars ..
00120       LOGICAL            NOUNIT
00121       INTEGER            J, K, KP
00122       COMPLEX            D11, D12, D21, D22, T1, T2
00123 *     ..
00124 *     .. External Functions ..
00125       LOGICAL            LSAME
00126       EXTERNAL           LSAME
00127 *     ..
00128 *     .. External Subroutines ..
00129       EXTERNAL           CGEMV, CGERU, CLACGV, CSCAL, CSWAP, XERBLA
00130 *     ..
00131 *     .. Intrinsic Functions ..
00132       INTRINSIC          ABS, CONJG, 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( 'CLAVHE ', -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 CSCAL( 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 CGERU( 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 CSWAP( 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 = CONJG( 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 CGERU( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ),
00234      $                        LDB, B( 1, 1 ), LDB )
00235                   CALL CGERU( 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 CSWAP( 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 CSCAL( 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 CGERU( N-K, NRHS, ONE, A( K+1, K ), 1,
00281      $                        B( K, 1 ), 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 CSWAP( 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 = CONJG( 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 CGERU( N-K, NRHS, ONE, A( K+1, K ), 1,
00317      $                        B( K, 1 ), LDB, B( K+1, 1 ), LDB )
00318                   CALL CGERU( 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 CSWAP( 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       IF( K.LT.1 )
00350      $         GO TO 90
00351 *
00352 *           1 x 1 pivot block.
00353 *
00354             IF( IPIV( K ).GT.0 ) THEN
00355                IF( K.GT.1 ) THEN
00356 *
00357 *                 Interchange if P(K) != I.
00358 *
00359                   KP = IPIV( K )
00360                   IF( KP.NE.K )
00361      $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00362 *
00363 *                 Apply the transformation
00364 *                    y = y - B' conjg(x),
00365 *                 where x is a column of A and y is a row of B.
00366 *
00367                   CALL CLACGV( NRHS, B( K, 1 ), LDB )
00368                   CALL CGEMV( 'Conjugate', K-1, NRHS, ONE, B, LDB,
00369      $                        A( 1, K ), 1, ONE, B( K, 1 ), LDB )
00370                   CALL CLACGV( NRHS, B( K, 1 ), LDB )
00371                END IF
00372                IF( NOUNIT )
00373      $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
00374                K = K - 1
00375 *
00376 *           2 x 2 pivot block.
00377 *
00378             ELSE
00379                IF( K.GT.2 ) THEN
00380 *
00381 *                 Interchange if P(K) != I.
00382 *
00383                   KP = ABS( IPIV( K ) )
00384                   IF( KP.NE.K-1 )
00385      $               CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
00386      $                           LDB )
00387 *
00388 *                 Apply the transformations
00389 *                    y = y - B' conjg(x),
00390 *                 where x is a block column of A and y is a block
00391 *                 row of B.
00392 *
00393                   CALL CLACGV( NRHS, B( K, 1 ), LDB )
00394                   CALL CGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
00395      $                        A( 1, K ), 1, ONE, B( K, 1 ), LDB )
00396                   CALL CLACGV( NRHS, B( K, 1 ), LDB )
00397 *
00398                   CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
00399                   CALL CGEMV( 'Conjugate', K-2, NRHS, ONE, B, LDB,
00400      $                        A( 1, K-1 ), 1, ONE, B( K-1, 1 ), LDB )
00401                   CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
00402                END IF
00403 *
00404 *              Multiply by the diagonal block if non-unit.
00405 *
00406                IF( NOUNIT ) THEN
00407                   D11 = A( K-1, K-1 )
00408                   D22 = A( K, K )
00409                   D12 = A( K-1, K )
00410                   D21 = CONJG( D12 )
00411                   DO 80 J = 1, NRHS
00412                      T1 = B( K-1, J )
00413                      T2 = B( K, J )
00414                      B( K-1, J ) = D11*T1 + D12*T2
00415                      B( K, J ) = D21*T1 + D22*T2
00416    80             CONTINUE
00417                END IF
00418                K = K - 2
00419             END IF
00420             GO TO 70
00421    90       CONTINUE
00422 *
00423 *        Form  B := L^H*B
00424 *        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
00425 *        and   L^H = inv(L^H(m))*P(m)* ... *inv(L^H(1))*P(1)
00426 *
00427          ELSE
00428 *
00429 *           Loop forward applying the L-transformations.
00430 *
00431             K = 1
00432   100       CONTINUE
00433             IF( K.GT.N )
00434      $         GO TO 120
00435 *
00436 *           1 x 1 pivot block
00437 *
00438             IF( IPIV( K ).GT.0 ) THEN
00439                IF( K.LT.N ) THEN
00440 *
00441 *                 Interchange if P(K) != I.
00442 *
00443                   KP = IPIV( K )
00444                   IF( KP.NE.K )
00445      $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00446 *
00447 *                 Apply the transformation
00448 *
00449                   CALL CLACGV( NRHS, B( K, 1 ), LDB )
00450                   CALL CGEMV( 'Conjugate', N-K, NRHS, ONE, B( K+1, 1 ),
00451      $                        LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
00452                   CALL CLACGV( NRHS, B( K, 1 ), LDB )
00453                END IF
00454                IF( NOUNIT )
00455      $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
00456                K = K + 1
00457 *
00458 *           2 x 2 pivot block.
00459 *
00460             ELSE
00461                IF( K.LT.N-1 ) THEN
00462 *
00463 *              Interchange if P(K) != I.
00464 *
00465                   KP = ABS( IPIV( K ) )
00466                   IF( KP.NE.K+1 )
00467      $               CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
00468      $                           LDB )
00469 *
00470 *                 Apply the transformation
00471 *
00472                   CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
00473                   CALL CGEMV( 'Conjugate', N-K-1, NRHS, ONE,
00474      $                        B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, ONE,
00475      $                        B( K+1, 1 ), LDB )
00476                   CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
00477 *
00478                   CALL CLACGV( NRHS, B( K, 1 ), LDB )
00479                   CALL CGEMV( 'Conjugate', N-K-1, NRHS, ONE,
00480      $                        B( K+2, 1 ), LDB, A( K+2, K ), 1, ONE,
00481      $                        B( K, 1 ), LDB )
00482                   CALL CLACGV( NRHS, B( K, 1 ), LDB )
00483                END IF
00484 *
00485 *              Multiply by the diagonal block if non-unit.
00486 *
00487                IF( NOUNIT ) THEN
00488                   D11 = A( K, K )
00489                   D22 = A( K+1, K+1 )
00490                   D21 = A( K+1, K )
00491                   D12 = CONJG( D21 )
00492                   DO 110 J = 1, NRHS
00493                      T1 = B( K, J )
00494                      T2 = B( K+1, J )
00495                      B( K, J ) = D11*T1 + D12*T2
00496                      B( K+1, J ) = D21*T1 + D22*T2
00497   110             CONTINUE
00498                END IF
00499                K = K + 2
00500             END IF
00501             GO TO 100
00502   120       CONTINUE
00503          END IF
00504 *
00505       END IF
00506       RETURN
00507 *
00508 *     End of CLAVHE
00509 *
00510       END