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