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