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