![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLAVSP 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 DLAVSP( 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 * DOUBLE PRECISION A( * ), B( LDB, * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> DLAVSP 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 DSPTRF. 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 DOUBLE PRECISION 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 DSPTRF. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] IPIV 00092 *> \verbatim 00093 *> IPIV is INTEGER array, dimension (N) 00094 *> The pivot indices from DSPTRF. 00095 *> \endverbatim 00096 *> 00097 *> \param[in,out] B 00098 *> \verbatim 00099 *> B is DOUBLE PRECISION 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 double_lin 00128 * 00129 * ===================================================================== 00130 SUBROUTINE DLAVSP( 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 DOUBLE PRECISION A( * ), B( LDB, * ) 00145 * .. 00146 * 00147 * ===================================================================== 00148 * 00149 * .. Parameters .. 00150 DOUBLE PRECISION ONE 00151 PARAMETER ( ONE = 1.0D+0 ) 00152 * .. 00153 * .. Local Scalars .. 00154 LOGICAL NOUNIT 00155 INTEGER J, K, KC, KCNEXT, KP 00156 DOUBLE PRECISION D11, D12, D21, D22, T1, T2 00157 * .. 00158 * .. External Functions .. 00159 LOGICAL LSAME 00160 EXTERNAL LSAME 00161 * .. 00162 * .. External Subroutines .. 00163 EXTERNAL DGEMV, DGER, DSCAL, DSWAP, 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( 'DLAVSP ', -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 DSCAL( 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 DGER( 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 DSWAP( 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 DGER( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ), LDB, 00271 $ B( 1, 1 ), LDB ) 00272 CALL DGER( 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 DSWAP( 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 DSCAL( 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 DGER( 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 DSWAP( 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 DGER( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ), 00359 $ LDB, B( K+1, 1 ), LDB ) 00360 CALL DGER( 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 DSWAP( 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 DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00408 * 00409 * Apply the transformation 00410 * 00411 CALL DGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB, 00412 $ A( KC ), 1, ONE, B( K, 1 ), LDB ) 00413 END IF 00414 IF( NOUNIT ) 00415 $ CALL DSCAL( 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 DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), 00429 $ LDB ) 00430 * 00431 * Apply the transformations 00432 * 00433 CALL DGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB, 00434 $ A( KC ), 1, ONE, B( K, 1 ), LDB ) 00435 CALL DGEMV( '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 DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00483 * 00484 * Apply the transformation 00485 * 00486 CALL DGEMV( '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 DSCAL( 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 DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), 00505 $ LDB ) 00506 * 00507 * Apply the transformation 00508 * 00509 CALL DGEMV( '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 DGEMV( '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 DLAVSP 00542 * 00543 END