![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSFRK 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SSFRK + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssfrk.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssfrk.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssfrk.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SSFRK( TRANSR, UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, 00022 * C ) 00023 * 00024 * .. Scalar Arguments .. 00025 * REAL ALPHA, BETA 00026 * INTEGER K, LDA, N 00027 * CHARACTER TRANS, TRANSR, UPLO 00028 * .. 00029 * .. Array Arguments .. 00030 * REAL 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 *> SSFRK 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 REAL 00114 *> On entry, ALPHA specifies the scalar alpha. 00115 *> Unchanged on exit. 00116 *> \endverbatim 00117 *> 00118 *> \param[in] A 00119 *> \verbatim 00120 *> A is REAL array of 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 REAL 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 REAL 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 realOTHERcomputational 00164 * 00165 * ===================================================================== 00166 SUBROUTINE SSFRK( 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 REAL ALPHA, BETA 00176 INTEGER K, LDA, N 00177 CHARACTER TRANS, TRANSR, UPLO 00178 * .. 00179 * .. Array Arguments .. 00180 REAL A( LDA, * ), C( * ) 00181 * .. 00182 * 00183 * ===================================================================== 00184 * 00185 * .. Parameters .. 00186 REAL ONE, ZERO 00187 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00188 * .. 00189 * .. Local Scalars .. 00190 LOGICAL LOWER, NORMALTRANSR, NISODD, NOTRANS 00191 INTEGER INFO, NROWA, J, NK, N1, N2 00192 * .. 00193 * .. External Functions .. 00194 LOGICAL LSAME 00195 EXTERNAL LSAME 00196 * .. 00197 * .. External Subroutines .. 00198 EXTERNAL SGEMM, SSYRK, XERBLA 00199 * .. 00200 * .. Intrinsic Functions .. 00201 INTRINSIC MAX 00202 * .. 00203 * .. Executable Statements .. 00204 * 00205 * Test the input parameters. 00206 * 00207 INFO = 0 00208 NORMALTRANSR = LSAME( TRANSR, 'N' ) 00209 LOWER = LSAME( UPLO, 'L' ) 00210 NOTRANS = LSAME( TRANS, 'N' ) 00211 * 00212 IF( NOTRANS ) THEN 00213 NROWA = N 00214 ELSE 00215 NROWA = K 00216 END IF 00217 * 00218 IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN 00219 INFO = -1 00220 ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN 00221 INFO = -2 00222 ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN 00223 INFO = -3 00224 ELSE IF( N.LT.0 ) THEN 00225 INFO = -4 00226 ELSE IF( K.LT.0 ) THEN 00227 INFO = -5 00228 ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THEN 00229 INFO = -8 00230 END IF 00231 IF( INFO.NE.0 ) THEN 00232 CALL XERBLA( 'SSFRK ', -INFO ) 00233 RETURN 00234 END IF 00235 * 00236 * Quick return if possible. 00237 * 00238 * The quick return case: ((ALPHA.EQ.0).AND.(BETA.NE.ZERO)) is not 00239 * done (it is in SSYRK for example) and left in the general case. 00240 * 00241 IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND. 00242 $ ( BETA.EQ.ONE ) ) )RETURN 00243 * 00244 IF( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ZERO ) ) THEN 00245 DO J = 1, ( ( N*( N+1 ) ) / 2 ) 00246 C( J ) = ZERO 00247 END DO 00248 RETURN 00249 END IF 00250 * 00251 * C is N-by-N. 00252 * If N is odd, set NISODD = .TRUE., and N1 and N2. 00253 * If N is even, NISODD = .FALSE., and NK. 00254 * 00255 IF( MOD( N, 2 ).EQ.0 ) THEN 00256 NISODD = .FALSE. 00257 NK = N / 2 00258 ELSE 00259 NISODD = .TRUE. 00260 IF( LOWER ) THEN 00261 N2 = N / 2 00262 N1 = N - N2 00263 ELSE 00264 N1 = N / 2 00265 N2 = N - N1 00266 END IF 00267 END IF 00268 * 00269 IF( NISODD ) THEN 00270 * 00271 * N is odd 00272 * 00273 IF( NORMALTRANSR ) THEN 00274 * 00275 * N is odd and TRANSR = 'N' 00276 * 00277 IF( LOWER ) THEN 00278 * 00279 * N is odd, TRANSR = 'N', and UPLO = 'L' 00280 * 00281 IF( NOTRANS ) THEN 00282 * 00283 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' 00284 * 00285 CALL SSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00286 $ BETA, C( 1 ), N ) 00287 CALL SSYRK( 'U', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 00288 $ BETA, C( N+1 ), N ) 00289 CALL SGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ), 00290 $ LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N ) 00291 * 00292 ELSE 00293 * 00294 * N is odd, TRANSR = 'N', UPLO = 'L', and TRANS = 'T' 00295 * 00296 CALL SSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA, 00297 $ BETA, C( 1 ), N ) 00298 CALL SSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA, 00299 $ BETA, C( N+1 ), N ) 00300 CALL SGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ), 00301 $ LDA, A( 1, 1 ), LDA, BETA, C( N1+1 ), N ) 00302 * 00303 END IF 00304 * 00305 ELSE 00306 * 00307 * N is odd, TRANSR = 'N', and UPLO = 'U' 00308 * 00309 IF( NOTRANS ) THEN 00310 * 00311 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' 00312 * 00313 CALL SSYRK( 'L', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00314 $ BETA, C( N2+1 ), N ) 00315 CALL SSYRK( 'U', 'N', N2, K, ALPHA, A( N2, 1 ), LDA, 00316 $ BETA, C( N1+1 ), N ) 00317 CALL SGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ), 00318 $ LDA, A( N2, 1 ), LDA, BETA, C( 1 ), N ) 00319 * 00320 ELSE 00321 * 00322 * N is odd, TRANSR = 'N', UPLO = 'U', and TRANS = 'T' 00323 * 00324 CALL SSYRK( 'L', 'T', N1, K, ALPHA, A( 1, 1 ), LDA, 00325 $ BETA, C( N2+1 ), N ) 00326 CALL SSYRK( 'U', 'T', N2, K, ALPHA, A( 1, N2 ), LDA, 00327 $ BETA, C( N1+1 ), N ) 00328 CALL SGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ), 00329 $ LDA, A( 1, N2 ), LDA, BETA, C( 1 ), N ) 00330 * 00331 END IF 00332 * 00333 END IF 00334 * 00335 ELSE 00336 * 00337 * N is odd, and TRANSR = 'T' 00338 * 00339 IF( LOWER ) THEN 00340 * 00341 * N is odd, TRANSR = 'T', and UPLO = 'L' 00342 * 00343 IF( NOTRANS ) THEN 00344 * 00345 * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'N' 00346 * 00347 CALL SSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00348 $ BETA, C( 1 ), N1 ) 00349 CALL SSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 00350 $ BETA, C( 2 ), N1 ) 00351 CALL SGEMM( 'N', 'T', N1, N2, K, ALPHA, A( 1, 1 ), 00352 $ LDA, A( N1+1, 1 ), LDA, BETA, 00353 $ C( N1*N1+1 ), N1 ) 00354 * 00355 ELSE 00356 * 00357 * N is odd, TRANSR = 'T', UPLO = 'L', and TRANS = 'T' 00358 * 00359 CALL SSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA, 00360 $ BETA, C( 1 ), N1 ) 00361 CALL SSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA, 00362 $ BETA, C( 2 ), N1 ) 00363 CALL SGEMM( 'T', 'N', N1, N2, K, ALPHA, A( 1, 1 ), 00364 $ LDA, A( 1, N1+1 ), LDA, BETA, 00365 $ C( N1*N1+1 ), N1 ) 00366 * 00367 END IF 00368 * 00369 ELSE 00370 * 00371 * N is odd, TRANSR = 'T', and UPLO = 'U' 00372 * 00373 IF( NOTRANS ) THEN 00374 * 00375 * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'N' 00376 * 00377 CALL SSYRK( 'U', 'N', N1, K, ALPHA, A( 1, 1 ), LDA, 00378 $ BETA, C( N2*N2+1 ), N2 ) 00379 CALL SSYRK( 'L', 'N', N2, K, ALPHA, A( N1+1, 1 ), LDA, 00380 $ BETA, C( N1*N2+1 ), N2 ) 00381 CALL SGEMM( 'N', 'T', N2, N1, K, ALPHA, A( N1+1, 1 ), 00382 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 ) 00383 * 00384 ELSE 00385 * 00386 * N is odd, TRANSR = 'T', UPLO = 'U', and TRANS = 'T' 00387 * 00388 CALL SSYRK( 'U', 'T', N1, K, ALPHA, A( 1, 1 ), LDA, 00389 $ BETA, C( N2*N2+1 ), N2 ) 00390 CALL SSYRK( 'L', 'T', N2, K, ALPHA, A( 1, N1+1 ), LDA, 00391 $ BETA, C( N1*N2+1 ), N2 ) 00392 CALL SGEMM( 'T', 'N', N2, N1, K, ALPHA, A( 1, N1+1 ), 00393 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), N2 ) 00394 * 00395 END IF 00396 * 00397 END IF 00398 * 00399 END IF 00400 * 00401 ELSE 00402 * 00403 * N is even 00404 * 00405 IF( NORMALTRANSR ) THEN 00406 * 00407 * N is even and TRANSR = 'N' 00408 * 00409 IF( LOWER ) THEN 00410 * 00411 * N is even, TRANSR = 'N', and UPLO = 'L' 00412 * 00413 IF( NOTRANS ) THEN 00414 * 00415 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'N' 00416 * 00417 CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00418 $ BETA, C( 2 ), N+1 ) 00419 CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00420 $ BETA, C( 1 ), N+1 ) 00421 CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ), 00422 $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ), 00423 $ N+1 ) 00424 * 00425 ELSE 00426 * 00427 * N is even, TRANSR = 'N', UPLO = 'L', and TRANS = 'T' 00428 * 00429 CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA, 00430 $ BETA, C( 2 ), N+1 ) 00431 CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00432 $ BETA, C( 1 ), N+1 ) 00433 CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ), 00434 $ LDA, A( 1, 1 ), LDA, BETA, C( NK+2 ), 00435 $ N+1 ) 00436 * 00437 END IF 00438 * 00439 ELSE 00440 * 00441 * N is even, TRANSR = 'N', and UPLO = 'U' 00442 * 00443 IF( NOTRANS ) THEN 00444 * 00445 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'N' 00446 * 00447 CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00448 $ BETA, C( NK+2 ), N+1 ) 00449 CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00450 $ BETA, C( NK+1 ), N+1 ) 00451 CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ), 00452 $ LDA, A( NK+1, 1 ), LDA, BETA, C( 1 ), 00453 $ N+1 ) 00454 * 00455 ELSE 00456 * 00457 * N is even, TRANSR = 'N', UPLO = 'U', and TRANS = 'T' 00458 * 00459 CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, 1 ), LDA, 00460 $ BETA, C( NK+2 ), N+1 ) 00461 CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00462 $ BETA, C( NK+1 ), N+1 ) 00463 CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ), 00464 $ LDA, A( 1, NK+1 ), LDA, BETA, C( 1 ), 00465 $ N+1 ) 00466 * 00467 END IF 00468 * 00469 END IF 00470 * 00471 ELSE 00472 * 00473 * N is even, and TRANSR = 'T' 00474 * 00475 IF( LOWER ) THEN 00476 * 00477 * N is even, TRANSR = 'T', and UPLO = 'L' 00478 * 00479 IF( NOTRANS ) THEN 00480 * 00481 * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'N' 00482 * 00483 CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00484 $ BETA, C( NK+1 ), NK ) 00485 CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00486 $ BETA, C( 1 ), NK ) 00487 CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( 1, 1 ), 00488 $ LDA, A( NK+1, 1 ), LDA, BETA, 00489 $ C( ( ( NK+1 )*NK )+1 ), NK ) 00490 * 00491 ELSE 00492 * 00493 * N is even, TRANSR = 'T', UPLO = 'L', and TRANS = 'T' 00494 * 00495 CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA, 00496 $ BETA, C( NK+1 ), NK ) 00497 CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00498 $ BETA, C( 1 ), NK ) 00499 CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, 1 ), 00500 $ LDA, A( 1, NK+1 ), LDA, BETA, 00501 $ C( ( ( NK+1 )*NK )+1 ), NK ) 00502 * 00503 END IF 00504 * 00505 ELSE 00506 * 00507 * N is even, TRANSR = 'T', and UPLO = 'U' 00508 * 00509 IF( NOTRANS ) THEN 00510 * 00511 * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'N' 00512 * 00513 CALL SSYRK( 'U', 'N', NK, K, ALPHA, A( 1, 1 ), LDA, 00514 $ BETA, C( NK*( NK+1 )+1 ), NK ) 00515 CALL SSYRK( 'L', 'N', NK, K, ALPHA, A( NK+1, 1 ), LDA, 00516 $ BETA, C( NK*NK+1 ), NK ) 00517 CALL SGEMM( 'N', 'T', NK, NK, K, ALPHA, A( NK+1, 1 ), 00518 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK ) 00519 * 00520 ELSE 00521 * 00522 * N is even, TRANSR = 'T', UPLO = 'U', and TRANS = 'T' 00523 * 00524 CALL SSYRK( 'U', 'T', NK, K, ALPHA, A( 1, 1 ), LDA, 00525 $ BETA, C( NK*( NK+1 )+1 ), NK ) 00526 CALL SSYRK( 'L', 'T', NK, K, ALPHA, A( 1, NK+1 ), LDA, 00527 $ BETA, C( NK*NK+1 ), NK ) 00528 CALL SGEMM( 'T', 'N', NK, NK, K, ALPHA, A( 1, NK+1 ), 00529 $ LDA, A( 1, 1 ), LDA, BETA, C( 1 ), NK ) 00530 * 00531 END IF 00532 * 00533 END IF 00534 * 00535 END IF 00536 * 00537 END IF 00538 * 00539 RETURN 00540 * 00541 * End of SSFRK 00542 * 00543 END