 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLAVSP( 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 REAL A( * ), B( LDB, * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * SLAVSP 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 SSPTRF. 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) REAL 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 SSPTRF. 00062 * 00063 * IPIV (input) INTEGER array, dimension (N) 00064 * The pivot indices from SSPTRF. 00065 * 00066 * B (input/output) REAL 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 REAL ONE 00081 PARAMETER ( ONE = 1.0E+0 ) 00082 * .. 00083 * .. Local Scalars .. 00084 LOGICAL NOUNIT 00085 INTEGER J, K, KC, KCNEXT, KP 00086 REAL D11, D12, D21, D22, T1, T2 00087 * .. 00088 * .. External Functions .. 00089 LOGICAL LSAME 00090 EXTERNAL LSAME 00091 * .. 00092 * .. External Subroutines .. 00093 EXTERNAL SGEMV, SGER, SSCAL, SSWAP, 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( 'SLAVSP ', -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 SSCAL( 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 SGER( 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 SSWAP( 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 SGER( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ), LDB, 00201 $ B( 1, 1 ), LDB ) 00202 CALL SGER( 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 SSWAP( 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 SSCAL( 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 SGER( 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 SSWAP( 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 SGER( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ), 00289 $ LDB, B( K+1, 1 ), LDB ) 00290 CALL SGER( 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 SSWAP( 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 SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00338 * 00339 * Apply the transformation 00340 * 00341 CALL SGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB, 00342 $ A( KC ), 1, ONE, B( K, 1 ), LDB ) 00343 END IF 00344 IF( NOUNIT ) 00345 $ CALL SSCAL( 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 SSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), 00359 $ LDB ) 00360 * 00361 * Apply the transformations 00362 * 00363 CALL SGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB, 00364 $ A( KC ), 1, ONE, B( K, 1 ), LDB ) 00365 CALL SGEMV( '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 SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00413 * 00414 * Apply the transformation 00415 * 00416 CALL SGEMV( '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 SSCAL( 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 SSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), 00435 $ LDB ) 00436 * 00437 * Apply the transformation 00438 * 00439 CALL SGEMV( '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 SGEMV( '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 SLAVSP 00472 * 00473 END
1.7.3 
