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