LAPACK: TESTING/LIN/zlavsy.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLAVSY( 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 *     ZLAVSY  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 ZSYTRF.
00024 *     ZSYTRF 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', ZLAVSY multiplies either by U or U * D
00035 *     (or L or L * D).
00036 *     If TRANS = 'T' or 't', ZLAVSY 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*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, ZSCAL, ZSWAP
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( 'ZLAVSY ', -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 = 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 = 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' * 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       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 *
00366                   CALL ZGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB,
00367      $                        A( 1, K ), 1, ONE, B( K, 1 ), LDB )
00368                END IF
00369                IF( NOUNIT )
00370      $            CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
00371                K = K - 1
00372 *
00373 *           2 x 2 pivot block.
00374 *
00375             ELSE
00376                IF( K.GT.2 ) THEN
00377 *
00378 *                 Interchange if P(K) != I.
00379 *
00380                   KP = ABS( IPIV( K ) )
00381                   IF( KP.NE.K-1 )
00382      $               CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
00383      $                           LDB )
00384 *
00385 *                 Apply the transformations
00386 *
00387                   CALL ZGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
00388      $                        A( 1, K ), 1, ONE, B( K, 1 ), LDB )
00389                   CALL ZGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
00390      $                        A( 1, K-1 ), 1, ONE, B( K-1, 1 ), LDB )
00391                END IF
00392 *
00393 *              Multiply by the diagonal block if non-unit.
00394 *
00395                IF( NOUNIT ) THEN
00396                   D11 = A( K-1, K-1 )
00397                   D22 = A( K, K )
00398                   D12 = A( K-1, K )
00399                   D21 = D12
00400                   DO 80 J = 1, NRHS
00401                      T1 = B( K-1, J )
00402                      T2 = B( K, J )
00403                      B( K-1, J ) = D11*T1 + D12*T2
00404                      B( K, J ) = D21*T1 + D22*T2
00405    80             CONTINUE
00406                END IF
00407                K = K - 2
00408             END IF
00409             GO TO 70
00410    90       CONTINUE
00411 *
00412 *        Form  B := L'*B
00413 *        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
00414 *        and   L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1)
00415 *
00416          ELSE
00417 *
00418 *           Loop forward applying the L-transformations.
00419 *
00420             K = 1
00421   100       CONTINUE
00422             IF( K.GT.N )
00423      $         GO TO 120
00424 *
00425 *           1 x 1 pivot block
00426 *
00427             IF( IPIV( K ).GT.0 ) THEN
00428                IF( K.LT.N ) THEN
00429 *
00430 *                 Interchange if P(K) != I.
00431 *
00432                   KP = IPIV( K )
00433                   IF( KP.NE.K )
00434      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00435 *
00436 *                 Apply the transformation
00437 *
00438                   CALL ZGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ),
00439      $                        LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
00440                END IF
00441                IF( NOUNIT )
00442      $            CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
00443                K = K + 1
00444 *
00445 *           2 x 2 pivot block.
00446 *
00447             ELSE
00448                IF( K.LT.N-1 ) THEN
00449 *
00450 *              Interchange if P(K) != I.
00451 *
00452                   KP = ABS( IPIV( K ) )
00453                   IF( KP.NE.K+1 )
00454      $               CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
00455      $                           LDB )
00456 *
00457 *                 Apply the transformation
00458 *
00459                   CALL ZGEMV( 'Transpose', N-K-1, NRHS, ONE,
00460      $                        B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, ONE,
00461      $                        B( K+1, 1 ), LDB )
00462                   CALL ZGEMV( 'Transpose', N-K-1, NRHS, ONE,
00463      $                        B( K+2, 1 ), LDB, A( K+2, K ), 1, ONE,
00464      $                        B( K, 1 ), LDB )
00465                END IF
00466 *
00467 *              Multiply by the diagonal block if non-unit.
00468 *
00469                IF( NOUNIT ) THEN
00470                   D11 = A( K, K )
00471                   D22 = A( K+1, K+1 )
00472                   D21 = A( K+1, K )
00473                   D12 = D21
00474                   DO 110 J = 1, NRHS
00475                      T1 = B( K, J )
00476                      T2 = B( K+1, J )
00477                      B( K, J ) = D11*T1 + D12*T2
00478                      B( K+1, J ) = D21*T1 + D22*T2
00479   110             CONTINUE
00480                END IF
00481                K = K + 2
00482             END IF
00483             GO TO 100
00484   120       CONTINUE
00485          END IF
00486       END IF
00487       RETURN
00488 *
00489 *     End of ZLAVSY
00490 *
00491       END