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