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