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