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