LAPACK: TESTING/LIN/clavsy.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLAVSY( 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 *     CLAVSY  performs one of the matrix-vector operations
00021 *        x := A*x  or  x := A'*x,
00022 *     where x is an N element vector and  A is one of the factors
00023 *     from the symmetric factorization computed by CSYTRF.
00024 *     CSYTRF produces a factorization of the form
00025 *          U * D * U'      or     L * D * L' ,
00026 *     where U (or L) is a product of permutation and unit upper (lower)
00027 *     triangular matrices, U' (or L') is the transpose of
00028 *     U (or L), and D is symmetric 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', CLAVSY multiplies either by U or U * D
00035 *     (or L or L * D).
00036 *     If TRANS = 'T' or 't', CLAVSY multiplies either by U' or D * U'
00037 *     (or L' or D * L' ).
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 = 'T' or 't'   x := A'*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, CSCAL, CSWAP, XERBLA
00130 *     ..
00131 *     .. Intrinsic Functions ..
00132       INTRINSIC          ABS, 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, 'T' ) )
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( 'CLAVSY ', -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 = 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 = 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' * B  (transpose)
00336 *
00337 *----------------------------------------
00338       ELSE IF( LSAME( TRANS, 'T' ) ) THEN
00339 *
00340 *        Form  B := U'*B
00341 *        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
00342 *        and   U' = inv(U'(1))*P(1)* ... *inv(U'(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 *
00365                   CALL CGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB,
00366      $                        A( 1, K ), 1, ONE, B( K, 1 ), LDB )
00367                END IF
00368                IF( NOUNIT )
00369      $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
00370                K = K - 1
00371 *
00372 *           2 x 2 pivot block.
00373 *
00374             ELSE
00375                IF( K.GT.2 ) THEN
00376 *
00377 *                 Interchange if P(K) != I.
00378 *
00379                   KP = ABS( IPIV( K ) )
00380                   IF( KP.NE.K-1 )
00381      $               CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
00382      $                           LDB )
00383 *
00384 *                 Apply the transformations
00385 *
00386                   CALL CGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
00387      $                        A( 1, K ), 1, ONE, B( K, 1 ), LDB )
00388                   CALL CGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
00389      $                        A( 1, K-1 ), 1, ONE, B( K-1, 1 ), LDB )
00390                END IF
00391 *
00392 *              Multiply by the diagonal block if non-unit.
00393 *
00394                IF( NOUNIT ) THEN
00395                   D11 = A( K-1, K-1 )
00396                   D22 = A( K, K )
00397                   D12 = A( K-1, K )
00398                   D21 = D12
00399                   DO 80 J = 1, NRHS
00400                      T1 = B( K-1, J )
00401                      T2 = B( K, J )
00402                      B( K-1, J ) = D11*T1 + D12*T2
00403                      B( K, J ) = D21*T1 + D22*T2
00404    80             CONTINUE
00405                END IF
00406                K = K - 2
00407             END IF
00408             GO TO 70
00409    90       CONTINUE
00410 *
00411 *        Form  B := L'*B
00412 *        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
00413 *        and   L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1)
00414 *
00415          ELSE
00416 *
00417 *           Loop forward applying the L-transformations.
00418 *
00419             K = 1
00420   100       CONTINUE
00421             IF( K.GT.N )
00422      $         GO TO 120
00423 *
00424 *           1 x 1 pivot block
00425 *
00426             IF( IPIV( K ).GT.0 ) THEN
00427                IF( K.LT.N ) THEN
00428 *
00429 *                 Interchange if P(K) != I.
00430 *
00431                   KP = IPIV( K )
00432                   IF( KP.NE.K )
00433      $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00434 *
00435 *                 Apply the transformation
00436 *
00437                   CALL CGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ),
00438      $                        LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
00439                END IF
00440                IF( NOUNIT )
00441      $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
00442                K = K + 1
00443 *
00444 *           2 x 2 pivot block.
00445 *
00446             ELSE
00447                IF( K.LT.N-1 ) THEN
00448 *
00449 *              Interchange if P(K) != I.
00450 *
00451                   KP = ABS( IPIV( K ) )
00452                   IF( KP.NE.K+1 )
00453      $               CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
00454      $                           LDB )
00455 *
00456 *                 Apply the transformation
00457 *
00458                   CALL CGEMV( 'Transpose', N-K-1, NRHS, ONE,
00459      $                        B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, ONE,
00460      $                        B( K+1, 1 ), LDB )
00461                   CALL CGEMV( 'Transpose', N-K-1, NRHS, ONE,
00462      $                        B( K+2, 1 ), LDB, A( K+2, K ), 1, ONE,
00463      $                        B( K, 1 ), LDB )
00464                END IF
00465 *
00466 *              Multiply by the diagonal block if non-unit.
00467 *
00468                IF( NOUNIT ) THEN
00469                   D11 = A( K, K )
00470                   D22 = A( K+1, K+1 )
00471                   D21 = A( K+1, K )
00472                   D12 = D21
00473                   DO 110 J = 1, NRHS
00474                      T1 = B( K, J )
00475                      T2 = B( K+1, J )
00476                      B( K, J ) = D11*T1 + D12*T2
00477                      B( K+1, J ) = D21*T1 + D22*T2
00478   110             CONTINUE
00479                END IF
00480                K = K + 2
00481             END IF
00482             GO TO 100
00483   120       CONTINUE
00484          END IF
00485       END IF
00486       RETURN
00487 *
00488 *     End of CLAVSY
00489 *
00490       END