|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
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
1.7.3 
