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