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