LAPACK: TESTING/LIN/dlavsy.f Source File

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