![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLAVSY 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 SLAVSY( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B, 00012 * LDB, INFO ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER DIAG, TRANS, UPLO 00016 * INTEGER INFO, LDA, LDB, N, NRHS 00017 * .. 00018 * .. Array Arguments .. 00019 * INTEGER IPIV( * ) 00020 * REAL A( LDA, * ), B( LDB, * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> SLAVSY 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 SSYTRF. 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 REAL array, dimension (LDA,N) 00086 *> The block diagonal matrix D and the multipliers used to 00087 *> obtain the factor U or L as computed by SSYTRF. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] LDA 00091 *> \verbatim 00092 *> LDA is INTEGER 00093 *> The leading dimension of the array A. LDA >= max(1,N). 00094 *> \endverbatim 00095 *> 00096 *> \param[in] IPIV 00097 *> \verbatim 00098 *> IPIV is INTEGER array, dimension (N) 00099 *> Details of the interchanges and the block structure of D, 00100 *> as determined by SSYTRF. 00101 *> 00102 *> If UPLO = 'U': 00103 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) 00104 *> were interchanged and D(k,k) is a 1-by-1 diagonal block. 00105 *> (If IPIV( k ) = k, no interchange was done). 00106 *> 00107 *> If IPIV(k) = IPIV(k-1) < 0, then rows and 00108 *> columns k-1 and -IPIV(k) were interchanged, 00109 *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. 00110 *> 00111 *> If UPLO = 'L': 00112 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) 00113 *> were interchanged and D(k,k) is a 1-by-1 diagonal block. 00114 *> (If IPIV( k ) = k, no interchange was done). 00115 *> 00116 *> If IPIV(k) = IPIV(k+1) < 0, then rows and 00117 *> columns k+1 and -IPIV(k) were interchanged, 00118 *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 00119 *> \endverbatim 00120 *> 00121 *> \param[in,out] B 00122 *> \verbatim 00123 *> B is REAL array, dimension (LDB,NRHS) 00124 *> On entry, B contains NRHS vectors of length N. 00125 *> On exit, B is overwritten with the product A * B. 00126 *> \endverbatim 00127 *> 00128 *> \param[in] LDB 00129 *> \verbatim 00130 *> LDB is INTEGER 00131 *> The leading dimension of the array B. LDB >= max(1,N). 00132 *> \endverbatim 00133 *> 00134 *> \param[out] INFO 00135 *> \verbatim 00136 *> INFO is INTEGER 00137 *> = 0: successful exit 00138 *> < 0: if INFO = -k, the k-th argument had an illegal value 00139 *> \endverbatim 00140 * 00141 * Authors: 00142 * ======== 00143 * 00144 *> \author Univ. of Tennessee 00145 *> \author Univ. of California Berkeley 00146 *> \author Univ. of Colorado Denver 00147 *> \author NAG Ltd. 00148 * 00149 *> \date April 2012 00150 * 00151 *> \ingroup single_lin 00152 * 00153 * ===================================================================== 00154 SUBROUTINE SLAVSY( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B, 00155 $ LDB, INFO ) 00156 * 00157 * -- LAPACK test routine (version 3.4.1) -- 00158 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00159 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00160 * April 2012 00161 * 00162 * .. Scalar Arguments .. 00163 CHARACTER DIAG, TRANS, UPLO 00164 INTEGER INFO, LDA, LDB, N, NRHS 00165 * .. 00166 * .. Array Arguments .. 00167 INTEGER IPIV( * ) 00168 REAL A( LDA, * ), B( LDB, * ) 00169 * .. 00170 * 00171 * ===================================================================== 00172 * 00173 * .. Parameters .. 00174 REAL ONE 00175 PARAMETER ( ONE = 1.0E+0 ) 00176 * .. 00177 * .. Local Scalars .. 00178 LOGICAL NOUNIT 00179 INTEGER J, K, KP 00180 REAL D11, D12, D21, D22, T1, T2 00181 * .. 00182 * .. External Functions .. 00183 LOGICAL LSAME 00184 EXTERNAL LSAME 00185 * .. 00186 * .. External Subroutines .. 00187 EXTERNAL SGEMV, SGER, SSCAL, SSWAP, XERBLA 00188 * .. 00189 * .. Intrinsic Functions .. 00190 INTRINSIC ABS, MAX 00191 * .. 00192 * .. Executable Statements .. 00193 * 00194 * Test the input parameters. 00195 * 00196 INFO = 0 00197 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00198 INFO = -1 00199 ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT. 00200 $ LSAME( TRANS, 'T' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 00201 INFO = -2 00202 ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) ) 00203 $ THEN 00204 INFO = -3 00205 ELSE IF( N.LT.0 ) THEN 00206 INFO = -4 00207 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00208 INFO = -6 00209 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00210 INFO = -9 00211 END IF 00212 IF( INFO.NE.0 ) THEN 00213 CALL XERBLA( 'SLAVSY ', -INFO ) 00214 RETURN 00215 END IF 00216 * 00217 * Quick return if possible. 00218 * 00219 IF( N.EQ.0 ) 00220 $ RETURN 00221 * 00222 NOUNIT = LSAME( DIAG, 'N' ) 00223 *------------------------------------------ 00224 * 00225 * Compute B := A * B (No transpose) 00226 * 00227 *------------------------------------------ 00228 IF( LSAME( TRANS, 'N' ) ) THEN 00229 * 00230 * Compute B := U*B 00231 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) 00232 * 00233 IF( LSAME( UPLO, 'U' ) ) THEN 00234 * 00235 * Loop forward applying the transformations. 00236 * 00237 K = 1 00238 10 CONTINUE 00239 IF( K.GT.N ) 00240 $ GO TO 30 00241 IF( IPIV( K ).GT.0 ) THEN 00242 * 00243 * 1 x 1 pivot block 00244 * 00245 * Multiply by the diagonal element if forming U * D. 00246 * 00247 IF( NOUNIT ) 00248 $ CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 00249 * 00250 * Multiply by P(K) * inv(U(K)) if K > 1. 00251 * 00252 IF( K.GT.1 ) THEN 00253 * 00254 * Apply the transformation. 00255 * 00256 CALL SGER( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ), 00257 $ LDB, B( 1, 1 ), LDB ) 00258 * 00259 * Interchange if P(K) .ne. I. 00260 * 00261 KP = IPIV( K ) 00262 IF( KP.NE.K ) 00263 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00264 END IF 00265 K = K + 1 00266 ELSE 00267 * 00268 * 2 x 2 pivot block 00269 * 00270 * Multiply by the diagonal block if forming U * D. 00271 * 00272 IF( NOUNIT ) THEN 00273 D11 = A( K, K ) 00274 D22 = A( K+1, K+1 ) 00275 D12 = A( K, K+1 ) 00276 D21 = D12 00277 DO 20 J = 1, NRHS 00278 T1 = B( K, J ) 00279 T2 = B( K+1, J ) 00280 B( K, J ) = D11*T1 + D12*T2 00281 B( K+1, J ) = D21*T1 + D22*T2 00282 20 CONTINUE 00283 END IF 00284 * 00285 * Multiply by P(K) * inv(U(K)) if K > 1. 00286 * 00287 IF( K.GT.1 ) THEN 00288 * 00289 * Apply the transformations. 00290 * 00291 CALL SGER( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ), 00292 $ LDB, B( 1, 1 ), LDB ) 00293 CALL SGER( K-1, NRHS, ONE, A( 1, K+1 ), 1, 00294 $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB ) 00295 * 00296 * Interchange if P(K) .ne. I. 00297 * 00298 KP = ABS( IPIV( K ) ) 00299 IF( KP.NE.K ) 00300 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00301 END IF 00302 K = K + 2 00303 END IF 00304 GO TO 10 00305 30 CONTINUE 00306 * 00307 * Compute B := L*B 00308 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) . 00309 * 00310 ELSE 00311 * 00312 * Loop backward applying the transformations to B. 00313 * 00314 K = N 00315 40 CONTINUE 00316 IF( K.LT.1 ) 00317 $ GO TO 60 00318 * 00319 * Test the pivot index. If greater than zero, a 1 x 1 00320 * pivot was used, otherwise a 2 x 2 pivot was used. 00321 * 00322 IF( IPIV( K ).GT.0 ) THEN 00323 * 00324 * 1 x 1 pivot block: 00325 * 00326 * Multiply by the diagonal element if forming L * D. 00327 * 00328 IF( NOUNIT ) 00329 $ CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 00330 * 00331 * Multiply by P(K) * inv(L(K)) if K < N. 00332 * 00333 IF( K.NE.N ) THEN 00334 KP = IPIV( K ) 00335 * 00336 * Apply the transformation. 00337 * 00338 CALL SGER( N-K, NRHS, ONE, A( K+1, K ), 1, B( K, 1 ), 00339 $ LDB, B( K+1, 1 ), LDB ) 00340 * 00341 * Interchange if a permutation was applied at the 00342 * K-th step of the factorization. 00343 * 00344 IF( KP.NE.K ) 00345 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00346 END IF 00347 K = K - 1 00348 * 00349 ELSE 00350 * 00351 * 2 x 2 pivot block: 00352 * 00353 * Multiply by the diagonal block if forming L * D. 00354 * 00355 IF( NOUNIT ) THEN 00356 D11 = A( K-1, K-1 ) 00357 D22 = A( K, K ) 00358 D21 = A( K, K-1 ) 00359 D12 = D21 00360 DO 50 J = 1, NRHS 00361 T1 = B( K-1, J ) 00362 T2 = B( K, J ) 00363 B( K-1, J ) = D11*T1 + D12*T2 00364 B( K, J ) = D21*T1 + D22*T2 00365 50 CONTINUE 00366 END IF 00367 * 00368 * Multiply by P(K) * inv(L(K)) if K < N. 00369 * 00370 IF( K.NE.N ) THEN 00371 * 00372 * Apply the transformation. 00373 * 00374 CALL SGER( N-K, NRHS, ONE, A( K+1, K ), 1, B( K, 1 ), 00375 $ LDB, B( K+1, 1 ), LDB ) 00376 CALL SGER( N-K, NRHS, ONE, A( K+1, K-1 ), 1, 00377 $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB ) 00378 * 00379 * Interchange if a permutation was applied at the 00380 * K-th step of the factorization. 00381 * 00382 KP = ABS( IPIV( K ) ) 00383 IF( KP.NE.K ) 00384 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00385 END IF 00386 K = K - 2 00387 END IF 00388 GO TO 40 00389 60 CONTINUE 00390 END IF 00391 *---------------------------------------- 00392 * 00393 * Compute B := A' * B (transpose) 00394 * 00395 *---------------------------------------- 00396 ELSE 00397 * 00398 * Form B := U'*B 00399 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) 00400 * and U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m) 00401 * 00402 IF( LSAME( UPLO, 'U' ) ) THEN 00403 * 00404 * Loop backward applying the transformations. 00405 * 00406 K = N 00407 70 CONTINUE 00408 IF( K.LT.1 ) 00409 $ GO TO 90 00410 * 00411 * 1 x 1 pivot block. 00412 * 00413 IF( IPIV( K ).GT.0 ) THEN 00414 IF( K.GT.1 ) THEN 00415 * 00416 * Interchange if P(K) .ne. I. 00417 * 00418 KP = IPIV( K ) 00419 IF( KP.NE.K ) 00420 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00421 * 00422 * Apply the transformation 00423 * 00424 CALL SGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB, 00425 $ A( 1, K ), 1, ONE, B( K, 1 ), LDB ) 00426 END IF 00427 IF( NOUNIT ) 00428 $ CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 00429 K = K - 1 00430 * 00431 * 2 x 2 pivot block. 00432 * 00433 ELSE 00434 IF( K.GT.2 ) THEN 00435 * 00436 * Interchange if P(K) .ne. I. 00437 * 00438 KP = ABS( IPIV( K ) ) 00439 IF( KP.NE.K-1 ) 00440 $ CALL SSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), 00441 $ LDB ) 00442 * 00443 * Apply the transformations 00444 * 00445 CALL SGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB, 00446 $ A( 1, K ), 1, ONE, B( K, 1 ), LDB ) 00447 CALL SGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB, 00448 $ A( 1, K-1 ), 1, ONE, B( K-1, 1 ), LDB ) 00449 END IF 00450 * 00451 * Multiply by the diagonal block if non-unit. 00452 * 00453 IF( NOUNIT ) THEN 00454 D11 = A( K-1, K-1 ) 00455 D22 = A( K, K ) 00456 D12 = A( K-1, K ) 00457 D21 = D12 00458 DO 80 J = 1, NRHS 00459 T1 = B( K-1, J ) 00460 T2 = B( K, J ) 00461 B( K-1, J ) = D11*T1 + D12*T2 00462 B( K, J ) = D21*T1 + D22*T2 00463 80 CONTINUE 00464 END IF 00465 K = K - 2 00466 END IF 00467 GO TO 70 00468 90 CONTINUE 00469 * 00470 * Form B := L'*B 00471 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) 00472 * and L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1) 00473 * 00474 ELSE 00475 * 00476 * Loop forward applying the L-transformations. 00477 * 00478 K = 1 00479 100 CONTINUE 00480 IF( K.GT.N ) 00481 $ GO TO 120 00482 * 00483 * 1 x 1 pivot block 00484 * 00485 IF( IPIV( K ).GT.0 ) THEN 00486 IF( K.LT.N ) THEN 00487 * 00488 * Interchange if P(K) .ne. I. 00489 * 00490 KP = IPIV( K ) 00491 IF( KP.NE.K ) 00492 $ CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) 00493 * 00494 * Apply the transformation 00495 * 00496 CALL SGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ), 00497 $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB ) 00498 END IF 00499 IF( NOUNIT ) 00500 $ CALL SSCAL( NRHS, A( K, K ), B( K, 1 ), LDB ) 00501 K = K + 1 00502 * 00503 * 2 x 2 pivot block. 00504 * 00505 ELSE 00506 IF( K.LT.N-1 ) THEN 00507 * 00508 * Interchange if P(K) .ne. I. 00509 * 00510 KP = ABS( IPIV( K ) ) 00511 IF( KP.NE.K+1 ) 00512 $ CALL SSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), 00513 $ LDB ) 00514 * 00515 * Apply the transformation 00516 * 00517 CALL SGEMV( 'Transpose', N-K-1, NRHS, ONE, 00518 $ B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, ONE, 00519 $ B( K+1, 1 ), LDB ) 00520 CALL SGEMV( 'Transpose', N-K-1, NRHS, ONE, 00521 $ B( K+2, 1 ), LDB, A( K+2, K ), 1, ONE, 00522 $ B( K, 1 ), LDB ) 00523 END IF 00524 * 00525 * Multiply by the diagonal block if non-unit. 00526 * 00527 IF( NOUNIT ) THEN 00528 D11 = A( K, K ) 00529 D22 = A( K+1, K+1 ) 00530 D21 = A( K+1, K ) 00531 D12 = D21 00532 DO 110 J = 1, NRHS 00533 T1 = B( K, J ) 00534 T2 = B( K+1, J ) 00535 B( K, J ) = D11*T1 + D12*T2 00536 B( K+1, J ) = D21*T1 + D22*T2 00537 110 CONTINUE 00538 END IF 00539 K = K + 2 00540 END IF 00541 GO TO 100 00542 120 CONTINUE 00543 END IF 00544 * 00545 END IF 00546 RETURN 00547 * 00548 * End of SLAVSY 00549 * 00550 END