LAPACK: TESTING/LIN/dlavsp.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAVSP( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
00002      $                   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, LDB, N, NRHS
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IPIV( * )
00014       DOUBLE PRECISION   A( * ), B( LDB, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DLAVSP  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 block U*D*U' or L*D*L' factorization computed by DSPTRF.
00024 *
00025 *  If TRANS = 'N', multiplies by U  or U * D  (or L  or L * D)
00026 *  If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L' )
00027 *  If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L' )
00028 *
00029 *  Arguments
00030 *  ==========
00031 *
00032 *  UPLO    (input) CHARACTER*1
00033 *          Specifies whether the factor stored in A is upper or lower
00034 *          triangular.
00035 *          = 'U':  Upper triangular
00036 *          = 'L':  Lower triangular
00037 *
00038 *  TRANS   (input) CHARACTER*1
00039 *          Specifies the operation to be performed:
00040 *          = 'N':  x := A*x
00041 *          = 'T':  x := A'*x
00042 *          = 'C':  x := A'*x
00043 *
00044 *  DIAG    (input) CHARACTER*1
00045 *          Specifies whether or not the diagonal blocks are unit
00046 *          matrices.  If the diagonal blocks are assumed to be unit,
00047 *          then A = U or A = L, otherwise A = U*D or A = L*D.
00048 *          = 'U':  Diagonal blocks are assumed to be unit matrices.
00049 *          = 'N':  Diagonal blocks are assumed to be non-unit matrices.
00050 *
00051 *  N       (input) INTEGER
00052 *          The number of rows and columns of the matrix A.  N >= 0.
00053 *
00054 *  NRHS    (input) INTEGER
00055 *          The number of right hand sides, i.e., the number of vectors
00056 *          x to be multiplied by A.  NRHS >= 0.
00057 *
00058 *  A       (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00059 *          The block diagonal matrix D and the multipliers used to
00060 *          obtain the factor U or L, stored as a packed triangular
00061 *          matrix as computed by DSPTRF.
00062 *
00063 *  IPIV    (input) INTEGER array, dimension (N)
00064 *          The pivot indices from DSPTRF.
00065 *
00066 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
00067 *          On entry, B contains NRHS vectors of length N.
00068 *          On exit, B is overwritten with the product A * B.
00069 *
00070 *  LDB     (input) INTEGER
00071 *          The leading dimension of the array B.  LDB >= max(1,N).
00072 *
00073 *  INFO    (output) INTEGER
00074 *          = 0: successful exit
00075 *          < 0: if INFO = -k, the k-th argument had an illegal value
00076 *
00077 *  =====================================================================
00078 *
00079 *     .. Parameters ..
00080       DOUBLE PRECISION   ONE
00081       PARAMETER          ( ONE = 1.0D+0 )
00082 *     ..
00083 *     .. Local Scalars ..
00084       LOGICAL            NOUNIT
00085       INTEGER            J, K, KC, KCNEXT, KP
00086       DOUBLE PRECISION   D11, D12, D21, D22, T1, T2
00087 *     ..
00088 *     .. External Functions ..
00089       LOGICAL            LSAME
00090       EXTERNAL           LSAME
00091 *     ..
00092 *     .. External Subroutines ..
00093       EXTERNAL           DGEMV, DGER, DSCAL, DSWAP, XERBLA
00094 *     ..
00095 *     .. Intrinsic Functions ..
00096       INTRINSIC          ABS, MAX
00097 *     ..
00098 *     .. Executable Statements ..
00099 *
00100 *     Test the input parameters.
00101 *
00102       INFO = 0
00103       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00104          INFO = -1
00105       ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.
00106      $         LSAME( TRANS, 'T' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
00107          INFO = -2
00108       ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
00109      $          THEN
00110          INFO = -3
00111       ELSE IF( N.LT.0 ) THEN
00112          INFO = -4
00113       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00114          INFO = -8
00115       END IF
00116       IF( INFO.NE.0 ) THEN
00117          CALL XERBLA( 'DLAVSP ', -INFO )
00118          RETURN
00119       END IF
00120 *
00121 *     Quick return if possible.
00122 *
00123       IF( N.EQ.0 )
00124      $   RETURN
00125 *
00126       NOUNIT = LSAME( DIAG, 'N' )
00127 *------------------------------------------
00128 *
00129 *     Compute  B := A * B  (No transpose)
00130 *
00131 *------------------------------------------
00132       IF( LSAME( TRANS, 'N' ) ) THEN
00133 *
00134 *        Compute  B := U*B
00135 *        where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
00136 *
00137          IF( LSAME( UPLO, 'U' ) ) THEN
00138 *
00139 *        Loop forward applying the transformations.
00140 *
00141             K = 1
00142             KC = 1
00143    10       CONTINUE
00144             IF( K.GT.N )
00145      $         GO TO 30
00146 *
00147 *           1 x 1 pivot block
00148 *
00149             IF( IPIV( K ).GT.0 ) THEN
00150 *
00151 *              Multiply by the diagonal element if forming U * D.
00152 *
00153                IF( NOUNIT )
00154      $            CALL DSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
00155 *
00156 *              Multiply by P(K) * inv(U(K))  if K > 1.
00157 *
00158                IF( K.GT.1 ) THEN
00159 *
00160 *                 Apply the transformation.
00161 *
00162                   CALL DGER( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ), LDB,
00163      $                       B( 1, 1 ), LDB )
00164 *
00165 *                 Interchange if P(K) != I.
00166 *
00167                   KP = IPIV( K )
00168                   IF( KP.NE.K )
00169      $               CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00170                END IF
00171                KC = KC + K
00172                K = K + 1
00173             ELSE
00174 *
00175 *              2 x 2 pivot block
00176 *
00177                KCNEXT = KC + K
00178 *
00179 *              Multiply by the diagonal block if forming U * D.
00180 *
00181                IF( NOUNIT ) THEN
00182                   D11 = A( KCNEXT-1 )
00183                   D22 = A( KCNEXT+K )
00184                   D12 = A( KCNEXT+K-1 )
00185                   D21 = D12
00186                   DO 20 J = 1, NRHS
00187                      T1 = B( K, J )
00188                      T2 = B( K+1, J )
00189                      B( K, J ) = D11*T1 + D12*T2
00190                      B( K+1, J ) = D21*T1 + D22*T2
00191    20             CONTINUE
00192                END IF
00193 *
00194 *              Multiply by  P(K) * inv(U(K))  if K > 1.
00195 *
00196                IF( K.GT.1 ) THEN
00197 *
00198 *                 Apply the transformations.
00199 *
00200                   CALL DGER( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ), LDB,
00201      $                       B( 1, 1 ), LDB )
00202                   CALL DGER( K-1, NRHS, ONE, A( KCNEXT ), 1,
00203      $                       B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
00204 *
00205 *                 Interchange if P(K) != I.
00206 *
00207                   KP = ABS( IPIV( K ) )
00208                   IF( KP.NE.K )
00209      $               CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00210                END IF
00211                KC = KCNEXT + K + 1
00212                K = K + 2
00213             END IF
00214             GO TO 10
00215    30       CONTINUE
00216 *
00217 *        Compute  B := L*B
00218 *        where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
00219 *
00220          ELSE
00221 *
00222 *           Loop backward applying the transformations to B.
00223 *
00224             K = N
00225             KC = N*( N+1 ) / 2 + 1
00226    40       CONTINUE
00227             IF( K.LT.1 )
00228      $         GO TO 60
00229             KC = KC - ( N-K+1 )
00230 *
00231 *           Test the pivot index.  If greater than zero, a 1 x 1
00232 *           pivot was used, otherwise a 2 x 2 pivot was used.
00233 *
00234             IF( IPIV( K ).GT.0 ) THEN
00235 *
00236 *              1 x 1 pivot block:
00237 *
00238 *              Multiply by the diagonal element if forming L * D.
00239 *
00240                IF( NOUNIT )
00241      $            CALL DSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
00242 *
00243 *              Multiply by  P(K) * inv(L(K))  if K < N.
00244 *
00245                IF( K.NE.N ) THEN
00246                   KP = IPIV( K )
00247 *
00248 *                 Apply the transformation.
00249 *
00250                   CALL DGER( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
00251      $                       LDB, B( K+1, 1 ), LDB )
00252 *
00253 *                 Interchange if a permutation was applied at the
00254 *                 K-th step of the factorization.
00255 *
00256                   IF( KP.NE.K )
00257      $               CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00258                END IF
00259                K = K - 1
00260 *
00261             ELSE
00262 *
00263 *              2 x 2 pivot block:
00264 *
00265                KCNEXT = KC - ( N-K+2 )
00266 *
00267 *              Multiply by the diagonal block if forming L * D.
00268 *
00269                IF( NOUNIT ) THEN
00270                   D11 = A( KCNEXT )
00271                   D22 = A( KC )
00272                   D21 = A( KCNEXT+1 )
00273                   D12 = D21
00274                   DO 50 J = 1, NRHS
00275                      T1 = B( K-1, J )
00276                      T2 = B( K, J )
00277                      B( K-1, J ) = D11*T1 + D12*T2
00278                      B( K, J ) = D21*T1 + D22*T2
00279    50             CONTINUE
00280                END IF
00281 *
00282 *              Multiply by  P(K) * inv(L(K))  if K < N.
00283 *
00284                IF( K.NE.N ) THEN
00285 *
00286 *                 Apply the transformation.
00287 *
00288                   CALL DGER( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
00289      $                       LDB, B( K+1, 1 ), LDB )
00290                   CALL DGER( N-K, NRHS, ONE, A( KCNEXT+2 ), 1,
00291      $                       B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
00292 *
00293 *                 Interchange if a permutation was applied at the
00294 *                 K-th step of the factorization.
00295 *
00296                   KP = ABS( IPIV( K ) )
00297                   IF( KP.NE.K )
00298      $               CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00299                END IF
00300                KC = KCNEXT
00301                K = K - 2
00302             END IF
00303             GO TO 40
00304    60       CONTINUE
00305          END IF
00306 *----------------------------------------
00307 *
00308 *     Compute  B := A' * B  (transpose)
00309 *
00310 *----------------------------------------
00311       ELSE
00312 *
00313 *        Form  B := U'*B
00314 *        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
00315 *        and   U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m)
00316 *
00317          IF( LSAME( UPLO, 'U' ) ) THEN
00318 *
00319 *           Loop backward applying the transformations.
00320 *
00321             K = N
00322             KC = N*( N+1 ) / 2 + 1
00323    70       CONTINUE
00324             IF( K.LT.1 )
00325      $         GO TO 90
00326             KC = KC - K
00327 *
00328 *           1 x 1 pivot block.
00329 *
00330             IF( IPIV( K ).GT.0 ) THEN
00331                IF( K.GT.1 ) THEN
00332 *
00333 *                 Interchange if P(K) != I.
00334 *
00335                   KP = IPIV( K )
00336                   IF( KP.NE.K )
00337      $               CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00338 *
00339 *                 Apply the transformation
00340 *
00341                   CALL DGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB,
00342      $                        A( KC ), 1, ONE, B( K, 1 ), LDB )
00343                END IF
00344                IF( NOUNIT )
00345      $            CALL DSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
00346                K = K - 1
00347 *
00348 *           2 x 2 pivot block.
00349 *
00350             ELSE
00351                KCNEXT = KC - ( K-1 )
00352                IF( K.GT.2 ) THEN
00353 *
00354 *                 Interchange if P(K) != I.
00355 *
00356                   KP = ABS( IPIV( K ) )
00357                   IF( KP.NE.K-1 )
00358      $               CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
00359      $                           LDB )
00360 *
00361 *                 Apply the transformations
00362 *
00363                   CALL DGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
00364      $                        A( KC ), 1, ONE, B( K, 1 ), LDB )
00365                   CALL DGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
00366      $                        A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
00367                END IF
00368 *
00369 *              Multiply by the diagonal block if non-unit.
00370 *
00371                IF( NOUNIT ) THEN
00372                   D11 = A( KC-1 )
00373                   D22 = A( KC+K-1 )
00374                   D12 = A( KC+K-2 )
00375                   D21 = D12
00376                   DO 80 J = 1, NRHS
00377                      T1 = B( K-1, J )
00378                      T2 = B( K, J )
00379                      B( K-1, J ) = D11*T1 + D12*T2
00380                      B( K, J ) = D21*T1 + D22*T2
00381    80             CONTINUE
00382                END IF
00383                KC = KCNEXT
00384                K = K - 2
00385             END IF
00386             GO TO 70
00387    90       CONTINUE
00388 *
00389 *        Form  B := L'*B
00390 *        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
00391 *        and   L' = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
00392 *
00393          ELSE
00394 *
00395 *           Loop forward applying the L-transformations.
00396 *
00397             K = 1
00398             KC = 1
00399   100       CONTINUE
00400             IF( K.GT.N )
00401      $         GO TO 120
00402 *
00403 *           1 x 1 pivot block
00404 *
00405             IF( IPIV( K ).GT.0 ) THEN
00406                IF( K.LT.N ) THEN
00407 *
00408 *                 Interchange if P(K) != I.
00409 *
00410                   KP = IPIV( K )
00411                   IF( KP.NE.K )
00412      $               CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00413 *
00414 *                 Apply the transformation
00415 *
00416                   CALL DGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ),
00417      $                        LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
00418                END IF
00419                IF( NOUNIT )
00420      $            CALL DSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
00421                KC = KC + N - K + 1
00422                K = K + 1
00423 *
00424 *           2 x 2 pivot block.
00425 *
00426             ELSE
00427                KCNEXT = KC + N - K + 1
00428                IF( K.LT.N-1 ) THEN
00429 *
00430 *              Interchange if P(K) != I.
00431 *
00432                   KP = ABS( IPIV( K ) )
00433                   IF( KP.NE.K+1 )
00434      $               CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
00435      $                           LDB )
00436 *
00437 *                 Apply the transformation
00438 *
00439                   CALL DGEMV( 'Transpose', N-K-1, NRHS, ONE,
00440      $                        B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
00441      $                        B( K+1, 1 ), LDB )
00442                   CALL DGEMV( 'Transpose', N-K-1, NRHS, ONE,
00443      $                        B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
00444      $                        B( K, 1 ), LDB )
00445                END IF
00446 *
00447 *              Multiply by the diagonal block if non-unit.
00448 *
00449                IF( NOUNIT ) THEN
00450                   D11 = A( KC )
00451                   D22 = A( KCNEXT )
00452                   D21 = A( KC+1 )
00453                   D12 = D21
00454                   DO 110 J = 1, NRHS
00455                      T1 = B( K, J )
00456                      T2 = B( K+1, J )
00457                      B( K, J ) = D11*T1 + D12*T2
00458                      B( K+1, J ) = D21*T1 + D22*T2
00459   110             CONTINUE
00460                END IF
00461                KC = KCNEXT + ( N-K )
00462                K = K + 2
00463             END IF
00464             GO TO 100
00465   120       CONTINUE
00466          END IF
00467 *
00468       END IF
00469       RETURN
00470 *
00471 *     End of DLAVSP
00472 *
00473       END