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