|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZHFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, 00002 + C ) 00003 * 00004 * -- LAPACK routine (version 3.3.0) -- 00005 * 00006 * -- Contributed by Julien Langou of the Univ. of Colorado Denver -- 00007 * November 2010 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*conjg( A' ) + beta*C, 00030 * 00031 * or 00032 * 00033 * C := alpha*conjg( A' )*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*conjg( A' ) + beta*C. 00064 * 00065 * TRANS = 'C' or 'c' C := alpha*conjg( A' )*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 * Arguments 00111 * ========== 00112 * 00113 * .. 00114 * .. Parameters .. 00115 DOUBLE PRECISION ONE, ZERO 00116 COMPLEX*16 CZERO 00117 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00118 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) 00119 * .. 00120 * .. Local Scalars .. 00121 LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS 00122 INTEGER INFO, NROWA, J, NK, N1, N2 00123 COMPLEX*16 CALPHA, CBETA 00124 * .. 00125 * .. External Functions .. 00126 LOGICAL LSAME 00127 EXTERNAL LSAME 00128 * .. 00129 * .. External Subroutines .. 00130 EXTERNAL XERBLA, ZGEMM, ZHERK 00131 * .. 00132 * .. Intrinsic Functions .. 00133 INTRINSIC MAX, DCMPLX 00134 * .. 00135 * .. Executable Statements .. 00136 * 00137 * 00138 * Test the input parameters. 00139 * 00140 INFO = 0 00141 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00142 LOWER = LSAME( UPLO, 'L' ) 00143 NOTRANS = LSAME( TRANS, 'N' ) 00144 * 00145 IF( NOTRANS ) THEN 00146 NROWA = N 00147 ELSE 00148 NROWA = K 00149 END IF 00150 * 00151 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'C' ) ) THEN 00152 INFO = -1 00153 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00154 INFO = -2 00155 ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 00156 INFO = -3 00157 ELSE IF( N.LT.0 ) THEN 00158 INFO = -4 00159 ELSE IF( K.LT.0 ) THEN 00160 INFO = -5 00161 ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN 00162 INFO = -8 00163 END IF 00164 IF( INFO.NE.0 ) THEN 00165 CALL XERBLA( 'ZHFRK ', -INFO ) 00166 RETURN 00167 END IF 00168 * 00169 * Quick return if possible. 00170 * 00171 * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not 00172 * done (it is in ZHERK for example) and left in the general case. 00173 * 00174 IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND. 00175 + ( BETA.EQ.ONE ) ) )RETURN 00176 * 00177 IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN 00178 DO J = 1, ( ( N*( N+1 ) ) / 2 ) 00179 C( J ) = CZERO 00180 END DO 00181 RETURN 00182 END IF 00183 * 00184 CALPHA = DCMPLX( ALPHA, ZERO ) 00185 CBETA = DCMPLX( BETA, ZERO ) 00186 * 00187 * C is N-by-N. 00188 * If N is odd, set NISODD = .TRUE., and N1 and N2. 00189 * If N is even, NISODD = .FALSE., and NK. 00190 * 00191 IF( MOD( N, 2 ).EQ.0 ) THEN 00192 NISODD = .FALSE. 00193 NK = N / 2 00194 ELSE 00195 NISODD = .TRUE. 00196 IF( LOWER ) THEN 00197 N2 = N / 2 00198 N1 = N - N2 00199 ELSE 00200 N1 = N / 2 00201 N2 = N - N1 00202 END IF 00203 END IF 00204 * 00205 IF( NISODD ) THEN 00206 * 00207 * N is odd 00208 * 00209 IF( NORMALTRANSR ) THEN 00210 * 00211 * N is odd and TRANSR = 'N' 00212 * 00213 IF( LOWER ) THEN 00214 * 00215 * N is odd, TRANSR = 'N', and UPLO = 'L' 00216 * 00217 IF( NOTRANS ) THEN 00218 * 00219 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' 00220 * 00221 CALL ZHERK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00222 + BETA, C( 1 ), N ) 00223 CALL ZHERK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 00224 + BETA, C( N+1 ), N ) 00225 CALL ZGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ), 00226 + LDA, A( 1, 1 ), LDA, CBETA, C( N1+1 ), N ) 00227 * 00228 ELSE 00229 * 00230 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'C' 00231 * 00232 CALL ZHERK( 'L', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 00233 + BETA, C( 1 ), N ) 00234 CALL ZHERK( 'U', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA, 00235 + BETA, C( N+1 ), N ) 00236 CALL ZGEMM( 'C', 'N', N2, N1, K, CALPHA, A( 1, N1+1 ), 00237 + LDA, A( 1, 1 ), LDA, CBETA, C( N1+1 ), N ) 00238 * 00239 END IF 00240 * 00241 ELSE 00242 * 00243 * N is odd, TRANSR = 'N', and UPLO = 'U' 00244 * 00245 IF( NOTRANS ) THEN 00246 * 00247 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' 00248 * 00249 CALL ZHERK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00250 + BETA, C( N2+1 ), N ) 00251 CALL ZHERK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA, 00252 + BETA, C( N1+1 ), N ) 00253 CALL ZGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ), 00254 + LDA, A( N2, 1 ), LDA, CBETA, C( 1 ), N ) 00255 * 00256 ELSE 00257 * 00258 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'C' 00259 * 00260 CALL ZHERK( 'L', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 00261 + BETA, C( N2+1 ), N ) 00262 CALL ZHERK( 'U', 'C', N2, K, ALPHA, A( 1, N2 ), LDA, 00263 + BETA, C( N1+1 ), N ) 00264 CALL ZGEMM( 'C', 'N', N1, N2, K, CALPHA, A( 1, 1 ), 00265 + LDA, A( 1, N2 ), LDA, CBETA, C( 1 ), N ) 00266 * 00267 END IF 00268 * 00269 END IF 00270 * 00271 ELSE 00272 * 00273 * N is odd, and TRANSR = 'C' 00274 * 00275 IF( LOWER ) THEN 00276 * 00277 * N is odd, TRANSR = 'C', and UPLO = 'L' 00278 * 00279 IF( NOTRANS ) THEN 00280 * 00281 * N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'N' 00282 * 00283 CALL ZHERK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00284 + BETA, C( 1 ), N1 ) 00285 CALL ZHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 00286 + BETA, C( 2 ), N1 ) 00287 CALL ZGEMM( 'N', 'C', N1, N2, K, CALPHA, A( 1, 1 ), 00288 + LDA, A( N1+1, 1 ), LDA, CBETA, 00289 + C( N1*N1+1 ), N1 ) 00290 * 00291 ELSE 00292 * 00293 * N is odd, TRANSR = 'C', UPLO = 'L', and TRANS = 'C' 00294 * 00295 CALL ZHERK( 'U', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 00296 + BETA, C( 1 ), N1 ) 00297 CALL ZHERK( 'L', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA, 00298 + BETA, C( 2 ), N1 ) 00299 CALL ZGEMM( 'C', 'N', N1, N2, K, CALPHA, A( 1, 1 ), 00300 + LDA, A( 1, N1+1 ), LDA, CBETA, 00301 + C( N1*N1+1 ), N1 ) 00302 * 00303 END IF 00304 * 00305 ELSE 00306 * 00307 * N is odd, TRANSR = 'C', and UPLO = 'U' 00308 * 00309 IF( NOTRANS ) THEN 00310 * 00311 * N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'N' 00312 * 00313 CALL ZHERK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00314 + BETA, C( N2*N2+1 ), N2 ) 00315 CALL ZHERK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 00316 + BETA, C( N1*N2+1 ), N2 ) 00317 CALL ZGEMM( 'N', 'C', N2, N1, K, CALPHA, A( N1+1, 1 ), 00318 + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), N2 ) 00319 * 00320 ELSE 00321 * 00322 * N is odd, TRANSR = 'C', UPLO = 'U', and TRANS = 'C' 00323 * 00324 CALL ZHERK( 'U', 'C', N1, K, ALPHA, A( 1, 1 ), LDA, 00325 + BETA, C( N2*N2+1 ), N2 ) 00326 CALL ZHERK( 'L', 'C', N2, K, ALPHA, A( 1, N1+1 ), LDA, 00327 + BETA, C( N1*N2+1 ), N2 ) 00328 CALL ZGEMM( 'C', 'N', N2, N1, K, CALPHA, A( 1, N1+1 ), 00329 + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), N2 ) 00330 * 00331 END IF 00332 * 00333 END IF 00334 * 00335 END IF 00336 * 00337 ELSE 00338 * 00339 * N is even 00340 * 00341 IF( NORMALTRANSR ) THEN 00342 * 00343 * N is even and TRANSR = 'N' 00344 * 00345 IF( LOWER ) THEN 00346 * 00347 * N is even, TRANSR = 'N', and UPLO = 'L' 00348 * 00349 IF( NOTRANS ) THEN 00350 * 00351 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' 00352 * 00353 CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00354 + BETA, C( 2 ), N+1 ) 00355 CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00356 + BETA, C( 1 ), N+1 ) 00357 CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ), 00358 + LDA, A( 1, 1 ), LDA, CBETA, C( NK+2 ), 00359 + N+1 ) 00360 * 00361 ELSE 00362 * 00363 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'C' 00364 * 00365 CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 00366 + BETA, C( 2 ), N+1 ) 00367 CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00368 + BETA, C( 1 ), N+1 ) 00369 CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, NK+1 ), 00370 + LDA, A( 1, 1 ), LDA, CBETA, C( NK+2 ), 00371 + N+1 ) 00372 * 00373 END IF 00374 * 00375 ELSE 00376 * 00377 * N is even, TRANSR = 'N', and UPLO = 'U' 00378 * 00379 IF( NOTRANS ) THEN 00380 * 00381 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' 00382 * 00383 CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00384 + BETA, C( NK+2 ), N+1 ) 00385 CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00386 + BETA, C( NK+1 ), N+1 ) 00387 CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ), 00388 + LDA, A( NK+1, 1 ), LDA, CBETA, C( 1 ), 00389 + N+1 ) 00390 * 00391 ELSE 00392 * 00393 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'C' 00394 * 00395 CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 00396 + BETA, C( NK+2 ), N+1 ) 00397 CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00398 + BETA, C( NK+1 ), N+1 ) 00399 CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, 1 ), 00400 + LDA, A( 1, NK+1 ), LDA, CBETA, C( 1 ), 00401 + N+1 ) 00402 * 00403 END IF 00404 * 00405 END IF 00406 * 00407 ELSE 00408 * 00409 * N is even, and TRANSR = 'C' 00410 * 00411 IF( LOWER ) THEN 00412 * 00413 * N is even, TRANSR = 'C', and UPLO = 'L' 00414 * 00415 IF( NOTRANS ) THEN 00416 * 00417 * N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'N' 00418 * 00419 CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00420 + BETA, C( NK+1 ), NK ) 00421 CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00422 + BETA, C( 1 ), NK ) 00423 CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( 1, 1 ), 00424 + LDA, A( NK+1, 1 ), LDA, CBETA, 00425 + C( ( ( NK+1 )*NK )+1 ), NK ) 00426 * 00427 ELSE 00428 * 00429 * N is even, TRANSR = 'C', UPLO = 'L', and TRANS = 'C' 00430 * 00431 CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 00432 + BETA, C( NK+1 ), NK ) 00433 CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00434 + BETA, C( 1 ), NK ) 00435 CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, 1 ), 00436 + LDA, A( 1, NK+1 ), LDA, CBETA, 00437 + C( ( ( NK+1 )*NK )+1 ), NK ) 00438 * 00439 END IF 00440 * 00441 ELSE 00442 * 00443 * N is even, TRANSR = 'C', and UPLO = 'U' 00444 * 00445 IF( NOTRANS ) THEN 00446 * 00447 * N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'N' 00448 * 00449 CALL ZHERK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00450 + BETA, C( NK*( NK+1 )+1 ), NK ) 00451 CALL ZHERK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00452 + BETA, C( NK*NK+1 ), NK ) 00453 CALL ZGEMM( 'N', 'C', NK, NK, K, CALPHA, A( NK+1, 1 ), 00454 + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), NK ) 00455 * 00456 ELSE 00457 * 00458 * N is even, TRANSR = 'C', UPLO = 'U', and TRANS = 'C' 00459 * 00460 CALL ZHERK( 'U', 'C', NK, K, ALPHA, A( 1, 1 ), LDA, 00461 + BETA, C( NK*( NK+1 )+1 ), NK ) 00462 CALL ZHERK( 'L', 'C', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00463 + BETA, C( NK*NK+1 ), NK ) 00464 CALL ZGEMM( 'C', 'N', NK, NK, K, CALPHA, A( 1, NK+1 ), 00465 + LDA, A( 1, 1 ), LDA, CBETA, C( 1 ), NK ) 00466 * 00467 END IF 00468 * 00469 END IF 00470 * 00471 END IF 00472 * 00473 END IF 00474 * 00475 RETURN 00476 * 00477 * End of ZHFRK 00478 * 00479 END
1.7.3 
