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