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