 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZHFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, 00002 $ C ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * 00006 * -- Contributed by Julien Langou of the Univ. of Colorado Denver -- 00007 * -- April 2011 -- 00008 * 00009 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00010 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00011 * 00012 * .. 00013 * .. Scalar Arguments .. 00014 DOUBLE PRECISION ALPHA, BETA 00015 INTEGER K, LDA, N 00016 CHARACTER TRANS, TRANSR, UPLO 00017 * .. 00018 * .. Array Arguments .. 00019 COMPLEX*16 A( LDA, * ), C( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * Level 3 BLAS like routine for C in RFP Format. 00026 * 00027 * ZHFRK performs one of the Hermitian rank--k operations 00028 * 00029 * C := alpha*A*A**H + beta*C, 00030 * 00031 * or 00032 * 00033 * C := alpha*A**H*A + beta*C, 00034 * 00035 * where alpha and beta are real scalars, C is an n--by--n Hermitian 00036 * matrix and A is an n--by--k matrix in the first case and a k--by--n 00037 * matrix in the second case. 00038 * 00039 * Arguments 00040 * ========== 00041 * 00042 * TRANSR (input) CHARACTER*1 00043 * = 'N': The Normal Form of RFP A is stored; 00044 * = 'C': The Conjugate-transpose Form of RFP A is stored. 00045 * 00046 * UPLO (input) CHARACTER*1 00047 * On entry, UPLO specifies whether the upper or lower 00048 * triangular part of the array C is to be referenced as 00049 * follows: 00050 * 00051 * UPLO = 'U' or 'u' Only the upper triangular part of C 00052 * is to be referenced. 00053 * 00054 * UPLO = 'L' or 'l' Only the lower triangular part of C 00055 * is to be referenced. 00056 * 00057 * Unchanged on exit. 00058 * 00059 * TRANS (input) CHARACTER*1 00060 * On entry, TRANS specifies the operation to be performed as 00061 * follows: 00062 * 00063 * TRANS = 'N' or 'n' C := alpha*A*A**H + beta*C. 00064 * 00065 * TRANS = 'C' or 'c' C := alpha*A**H*A + beta*C. 00066 * 00067 * Unchanged on exit. 00068 * 00069 * N (input) INTEGER 00070 * On entry, N specifies the order of the matrix C. N must be 00071 * at least zero. 00072 * Unchanged on exit. 00073 * 00074 * K (input) INTEGER 00075 * On entry with TRANS = 'N' or 'n', K specifies the number 00076 * of columns of the matrix A, and on entry with 00077 * TRANS = 'C' or 'c', K specifies the number of rows of the 00078 * matrix A. K must be at least zero. 00079 * Unchanged on exit. 00080 * 00081 * ALPHA (input) DOUBLE PRECISION 00082 * On entry, ALPHA specifies the scalar alpha. 00083 * Unchanged on exit. 00084 * 00085 * A (input) COMPLEX*16 array of DIMENSION (LDA,ka) 00086 * where KA 00087 * is K when TRANS = 'N' or 'n', and is N otherwise. Before 00088 * entry with TRANS = 'N' or 'n', the leading N--by--K part of 00089 * the array A must contain the matrix A, otherwise the leading 00090 * K--by--N part of the array A must contain the matrix A. 00091 * Unchanged on exit. 00092 * 00093 * LDA (input) INTEGER 00094 * On entry, LDA specifies the first dimension of A as declared 00095 * in the calling (sub) program. When TRANS = 'N' or 'n' 00096 * then LDA must be at least max( 1, n ), otherwise LDA must 00097 * be at least max( 1, k ). 00098 * Unchanged on exit. 00099 * 00100 * BETA (input) DOUBLE PRECISION 00101 * On entry, BETA specifies the scalar beta. 00102 * Unchanged on exit. 00103 * 00104 * C (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) 00105 * On entry, the matrix A in RFP Format. RFP Format is 00106 * described by TRANSR, UPLO and N. Note that the imaginary 00107 * parts of the diagonal elements need not be set, they are 00108 * assumed to be zero, and on exit they are set to zero. 00109 * 00110 * ===================================================================== 00111 * 00112 * .. Parameters .. 00113 DOUBLE PRECISION ONE, ZERO 00114 COMPLEX*16 CZERO 00115 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00116 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) 00117 * .. 00118 * .. Local Scalars .. 00119 LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS 00120 INTEGER INFO, NROWA, J, NK, N1, N2 00121 COMPLEX*16 CALPHA, CBETA 00122 * .. 00123 * .. External Functions .. 00124 LOGICAL LSAME 00125 EXTERNAL LSAME 00126 * .. 00127 * .. External Subroutines .. 00128 EXTERNAL XERBLA, ZGEMM, ZHERK 00129 * .. 00130 * .. Intrinsic Functions .. 00131 INTRINSIC MAX, DCMPLX 00132 * .. 00133 * .. Executable Statements .. 00134 * 00135 * 00136 * Test the input parameters. 00137 * 00138 INFO = 0 00139 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00140 LOWER = LSAME( UPLO, 'L' ) 00141 NOTRANS = LSAME( TRANS, 'N' ) 00142 * 00143 IF( NOTRANS ) THEN 00144 NROWA = N 00145 ELSE 00146 NROWA = K 00147 END IF 00148 * 00149 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 00150 INFO = -1 00151 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00152 INFO = -2 00153 ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 00154 INFO = -3 00155 ELSE IF( N.LT.0 ) THEN 00156 INFO = -4 00157 ELSE IF( K.LT.0 ) THEN 00158 INFO = -5 00159 ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN 00160 INFO = -8 00161 END IF 00162 IF( INFO.NE.0 ) THEN 00163 CALL XERBLA( 'ZHFRK ', -INFO ) 00164 RETURN 00165 END IF 00166 * 00167 * Quick return if possible. 00168 * 00169 * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not 00170 * done (it is in ZHERK for example) and left in the general case. 00171 * 00172 IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND. 00173 $ ( BETA.EQ.ONE ) ) )RETURN 00174 * 00175 IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN 00176 DO J = 1, ( ( N*( N+1 ) ) / 2 ) 00177 C( J ) = CZERO 00178 END DO 00179 RETURN 00180 END IF 00181 * 00182 CALPHA = DCMPLX( ALPHA, ZERO ) 00183 CBETA = DCMPLX( BETA, ZERO ) 00184 * 00185 * C is N-by-N. 00186 * If N is odd, set NISODD = .TRUE., and N1 and N2. 00187 * If N is even, NISODD = .FALSE., and NK. 00188 * 00189 IF( MOD( N, 2 ).EQ.0 ) THEN 00190 NISODD = .FALSE. 00191 NK = N / 2 00192 ELSE 00193 NISODD = .TRUE. 00194 IF( LOWER ) THEN 00195 N2 = N / 2 00196 N1 = N - N2 00197 ELSE 00198 N1 = N / 2 00199 N2 = N - N1 00200 END IF 00201 END IF 00202 * 00203 IF( NISODD ) THEN 00204 * 00205 * N is odd 00206 * 00207 IF( NORMALTRANSR ) THEN 00208 * 00209 * N is odd and TRANSR = 'N' 00210 * 00211 IF( LOWER ) THEN 00212 * 00213 * N is odd, TRANSR = 'N', and UPLO = 'L' 00214 * 00215 IF( NOTRANS ) THEN 00216 * 00217 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' 00218 * 00219 CALL ZHERK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00220 $ BETA, C( 1 ), N ) 00221 CALL ZHERK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 00222 $ BETA, C( N+1 ), N ) 00223 CALL ZGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ), 00224 $ LDA, A( 1, 1 ), LDA, CBETA, C( N1+1 ), N ) 00225 * 00226 ELSE 00227 * 00228 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'C' 00229 * 00230 CALL ZHERK( 'L', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 00231 $ BETA, C( 1 ), N ) 00232 CALL ZHERK( 'U', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA, 00233 $ BETA, C( N+1 ), N ) 00234 CALL ZGEMM( 'C', 'N', N2, N1, K, CALPHA, A( 1, N1+1 ), 00235 $ LDA, A( 1, 1 ), LDA, CBETA, C( N1+1 ), N ) 00236 * 00237 END IF 00238 * 00239 ELSE 00240 * 00241 * N is odd, TRANSR = 'N', and UPLO = 'U' 00242 * 00243 IF( NOTRANS ) THEN 00244 * 00245 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' 00246 * 00247 CALL ZHERK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00248 $ BETA, C( N2+1 ), N ) 00249 CALL ZHERK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA, 00250 $ BETA, C( N1+1 ), N ) 00251 CALL ZGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ), 00252 $ LDA, A( N2, 1 ), LDA, CBETA, C( 1 ), N ) 00253 * 00254 ELSE 00255 * 00256 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'C' 00257 * 00258 CALL ZHERK( 'L', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 00259 $ BETA, C( N2+1 ), N ) 00260 CALL ZHERK( 'U', 'C', N2, K, ALPHA, A( 1, N2 ), LDA, 00261 $ BETA, C( N1+1 ), N ) 00262 CALL ZGEMM( 'C', 'N', N1, N2, K, CALPHA, A( 1, 1 ), 00263 $ LDA, A( 1, N2 ), LDA, CBETA, C( 1 ), N ) 00264 * 00265 END IF 00266 * 00267 END IF 00268 * 00269 ELSE 00270 * 00271 * N is odd, and TRANSR = 'C' 00272 * 00273 IF( LOWER ) THEN 00274 * 00275 * N is odd, TRANSR = 'C', and UPLO = 'L' 00276 * 00277 IF( NOTRANS ) THEN 00278 * 00279 * N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'N' 00280 * 00281 CALL ZHERK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00282 $ BETA, C( 1 ), N1 ) 00283 CALL ZHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 00284 $ BETA, C( 2 ), N1 ) 00285 CALL ZGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ), 00286 $ LDA, A( N1+1, 1 ), LDA, CBETA, 00287 $ C( N1*N1+1 ), N1 ) 00288 * 00289 ELSE 00290 * 00291 * N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'C' 00292 * 00293 CALL ZHERK( 'U', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 00294 $ BETA, C( 1 ), N1 ) 00295 CALL ZHERK( 'L', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA, 00296 $ BETA, C( 2 ), N1 ) 00297 CALL ZGEMM( 'C', 'N', N1, N2, K, CALPHA, A( 1, 1 ), 00298 $ LDA, A( 1, N1+1 ), LDA, CBETA, 00299 $ C( N1*N1+1 ), N1 ) 00300 * 00301 END IF 00302 * 00303 ELSE 00304 * 00305 * N is odd, TRANSR = 'C', and UPLO = 'U' 00306 * 00307 IF( NOTRANS ) THEN 00308 * 00309 * N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'N' 00310 * 00311 CALL ZHERK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00312 $ BETA, C( N2*N2+1 ), N2 ) 00313 CALL ZHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 00314 $ BETA, C( N1*N2+1 ), N2 ) 00315 CALL ZGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ), 00316 $ LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), N2 ) 00317 * 00318 ELSE 00319 * 00320 * N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'C' 00321 * 00322 CALL ZHERK( 'U', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 00323 $ BETA, C( N2*N2+1 ), N2 ) 00324 CALL ZHERK( 'L', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA, 00325 $ BETA, C( N1*N2+1 ), N2 ) 00326 CALL ZGEMM( 'C', 'N', N2, N1, K, CALPHA, A( 1, N1+1 ), 00327 $ LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), N2 ) 00328 * 00329 END IF 00330 * 00331 END IF 00332 * 00333 END IF 00334 * 00335 ELSE 00336 * 00337 * N is even 00338 * 00339 IF( NORMALTRANSR ) THEN 00340 * 00341 * N is even and TRANSR = 'N' 00342 * 00343 IF( LOWER ) THEN 00344 * 00345 * N is even, TRANSR = 'N', and UPLO = 'L' 00346 * 00347 IF( NOTRANS ) THEN 00348 * 00349 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' 00350 * 00351 CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00352 $ BETA, C( 2 ), N+1 ) 00353 CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00354 $ BETA, C( 1 ), N+1 ) 00355 CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ), 00356 $ LDA, A( 1, 1 ), LDA, CBETA, C( NK+2 ), 00357 $ N+1 ) 00358 * 00359 ELSE 00360 * 00361 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'C' 00362 * 00363 CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 00364 $ BETA, C( 2 ), N+1 ) 00365 CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00366 $ BETA, C( 1 ), N+1 ) 00367 CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, NK+1 ), 00368 $ LDA, A( 1, 1 ), LDA, CBETA, C( NK+2 ), 00369 $ N+1 ) 00370 * 00371 END IF 00372 * 00373 ELSE 00374 * 00375 * N is even, TRANSR = 'N', and UPLO = 'U' 00376 * 00377 IF( NOTRANS ) THEN 00378 * 00379 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' 00380 * 00381 CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00382 $ BETA, C( NK+2 ), N+1 ) 00383 CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00384 $ BETA, C( NK+1 ), N+1 ) 00385 CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ), 00386 $ LDA, A( NK+1, 1 ), LDA, CBETA, C( 1 ), 00387 $ N+1 ) 00388 * 00389 ELSE 00390 * 00391 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'C' 00392 * 00393 CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 00394 $ BETA, C( NK+2 ), N+1 ) 00395 CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00396 $ BETA, C( NK+1 ), N+1 ) 00397 CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, 1 ), 00398 $ LDA, A( 1, NK+1 ), LDA, CBETA, C( 1 ), 00399 $ N+1 ) 00400 * 00401 END IF 00402 * 00403 END IF 00404 * 00405 ELSE 00406 * 00407 * N is even, and TRANSR = 'C' 00408 * 00409 IF( LOWER ) THEN 00410 * 00411 * N is even, TRANSR = 'C', and UPLO = 'L' 00412 * 00413 IF( NOTRANS ) THEN 00414 * 00415 * N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'N' 00416 * 00417 CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00418 $ BETA, C( NK+1 ), NK ) 00419 CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00420 $ BETA, C( 1 ), NK ) 00421 CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ), 00422 $ LDA, A( NK+1, 1 ), LDA, CBETA, 00423 $ C( ( ( NK+1 )*NK )+1 ), NK ) 00424 * 00425 ELSE 00426 * 00427 * N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'C' 00428 * 00429 CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 00430 $ BETA, C( NK+1 ), NK ) 00431 CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00432 $ BETA, C( 1 ), NK ) 00433 CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, 1 ), 00434 $ LDA, A( 1, NK+1 ), LDA, CBETA, 00435 $ C( ( ( NK+1 )*NK )+1 ), NK ) 00436 * 00437 END IF 00438 * 00439 ELSE 00440 * 00441 * N is even, TRANSR = 'C', and UPLO = 'U' 00442 * 00443 IF( NOTRANS ) THEN 00444 * 00445 * N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'N' 00446 * 00447 CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00448 $ BETA, C( NK*( NK+1 )+1 ), NK ) 00449 CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00450 $ BETA, C( NK*NK+1 ), NK ) 00451 CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ), 00452 $ LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), NK ) 00453 * 00454 ELSE 00455 * 00456 * N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'C' 00457 * 00458 CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 00459 $ BETA, C( NK*( NK+1 )+1 ), NK ) 00460 CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00461 $ BETA, C( NK*NK+1 ), NK ) 00462 CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, NK+1 ), 00463 $ LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), NK ) 00464 * 00465 END IF 00466 * 00467 END IF 00468 * 00469 END IF 00470 * 00471 END IF 00472 * 00473 RETURN 00474 * 00475 * End of ZHFRK 00476 * 00477 END
1.7.3 
