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