![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSYR2K 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) 00012 * 00013 * .. Scalar Arguments .. 00014 * DOUBLE PRECISION ALPHA,BETA 00015 * INTEGER K,LDA,LDB,LDC,N 00016 * CHARACTER TRANS,UPLO 00017 * .. 00018 * .. Array Arguments .. 00019 * DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) 00020 * .. 00021 * 00022 * 00023 *> \par Purpose: 00024 * ============= 00025 *> 00026 *> \verbatim 00027 *> 00028 *> DSYR2K performs one of the symmetric rank 2k operations 00029 *> 00030 *> C := alpha*A*B**T + alpha*B*A**T + beta*C, 00031 *> 00032 *> or 00033 *> 00034 *> C := alpha*A**T*B + alpha*B**T*A + beta*C, 00035 *> 00036 *> where alpha and beta are scalars, C is an n by n symmetric matrix 00037 *> and A and B are n by k matrices in the first case and k by n 00038 *> matrices in the second case. 00039 *> \endverbatim 00040 * 00041 * Arguments: 00042 * ========== 00043 * 00044 *> \param[in] UPLO 00045 *> \verbatim 00046 *> UPLO is 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 *> \endverbatim 00057 *> 00058 *> \param[in] TRANS 00059 *> \verbatim 00060 *> TRANS is CHARACTER*1 00061 *> On entry, TRANS specifies the operation to be performed as 00062 *> follows: 00063 *> 00064 *> TRANS = 'N' or 'n' C := alpha*A*B**T + alpha*B*A**T + 00065 *> beta*C. 00066 *> 00067 *> TRANS = 'T' or 't' C := alpha*A**T*B + alpha*B**T*A + 00068 *> beta*C. 00069 *> 00070 *> TRANS = 'C' or 'c' C := alpha*A**T*B + alpha*B**T*A + 00071 *> beta*C. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] N 00075 *> \verbatim 00076 *> N is INTEGER 00077 *> On entry, N specifies the order of the matrix C. N must be 00078 *> at least zero. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] K 00082 *> \verbatim 00083 *> K is INTEGER 00084 *> On entry with TRANS = 'N' or 'n', K specifies the number 00085 *> of columns of the matrices A and B, and on entry with 00086 *> TRANS = 'T' or 't' or 'C' or 'c', K specifies the number 00087 *> of rows of the matrices A and B. K must be at least zero. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] ALPHA 00091 *> \verbatim 00092 *> ALPHA is DOUBLE PRECISION. 00093 *> On entry, ALPHA specifies the scalar alpha. 00094 *> \endverbatim 00095 *> 00096 *> \param[in] A 00097 *> \verbatim 00098 *> A is DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is 00099 *> k when TRANS = 'N' or 'n', and is n otherwise. 00100 *> Before entry with TRANS = 'N' or 'n', the leading n by k 00101 *> part of the array A must contain the matrix A, otherwise 00102 *> the leading k by n part of the array A must contain the 00103 *> matrix A. 00104 *> \endverbatim 00105 *> 00106 *> \param[in] LDA 00107 *> \verbatim 00108 *> LDA is INTEGER 00109 *> On entry, LDA specifies the first dimension of A as declared 00110 *> in the calling (sub) program. When TRANS = 'N' or 'n' 00111 *> then LDA must be at least max( 1, n ), otherwise LDA must 00112 *> be at least max( 1, k ). 00113 *> \endverbatim 00114 *> 00115 *> \param[in] B 00116 *> \verbatim 00117 *> B is DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is 00118 *> k when TRANS = 'N' or 'n', and is n otherwise. 00119 *> Before entry with TRANS = 'N' or 'n', the leading n by k 00120 *> part of the array B must contain the matrix B, otherwise 00121 *> the leading k by n part of the array B must contain the 00122 *> matrix B. 00123 *> \endverbatim 00124 *> 00125 *> \param[in] LDB 00126 *> \verbatim 00127 *> LDB is INTEGER 00128 *> On entry, LDB specifies the first dimension of B as declared 00129 *> in the calling (sub) program. When TRANS = 'N' or 'n' 00130 *> then LDB must be at least max( 1, n ), otherwise LDB must 00131 *> be at least max( 1, k ). 00132 *> \endverbatim 00133 *> 00134 *> \param[in] BETA 00135 *> \verbatim 00136 *> BETA is DOUBLE PRECISION. 00137 *> On entry, BETA specifies the scalar beta. 00138 *> \endverbatim 00139 *> 00140 *> \param[in,out] C 00141 *> \verbatim 00142 *> C is DOUBLE PRECISION array of DIMENSION ( LDC, n ). 00143 *> Before entry with UPLO = 'U' or 'u', the leading n by n 00144 *> upper triangular part of the array C must contain the upper 00145 *> triangular part of the symmetric matrix and the strictly 00146 *> lower triangular part of C is not referenced. On exit, the 00147 *> upper triangular part of the array C is overwritten by the 00148 *> upper triangular part of the updated matrix. 00149 *> Before entry with UPLO = 'L' or 'l', the leading n by n 00150 *> lower triangular part of the array C must contain the lower 00151 *> triangular part of the symmetric matrix and the strictly 00152 *> upper triangular part of C is not referenced. On exit, the 00153 *> lower triangular part of the array C is overwritten by the 00154 *> lower triangular part of the updated matrix. 00155 *> \endverbatim 00156 *> 00157 *> \param[in] LDC 00158 *> \verbatim 00159 *> LDC is INTEGER 00160 *> On entry, LDC specifies the first dimension of C as declared 00161 *> in the calling (sub) program. LDC must be at least 00162 *> max( 1, n ). 00163 *> \endverbatim 00164 * 00165 * Authors: 00166 * ======== 00167 * 00168 *> \author Univ. of Tennessee 00169 *> \author Univ. of California Berkeley 00170 *> \author Univ. of Colorado Denver 00171 *> \author NAG Ltd. 00172 * 00173 *> \date November 2011 00174 * 00175 *> \ingroup double_blas_level3 00176 * 00177 *> \par Further Details: 00178 * ===================== 00179 *> 00180 *> \verbatim 00181 *> 00182 *> Level 3 Blas routine. 00183 *> 00184 *> 00185 *> -- Written on 8-February-1989. 00186 *> Jack Dongarra, Argonne National Laboratory. 00187 *> Iain Duff, AERE Harwell. 00188 *> Jeremy Du Croz, Numerical Algorithms Group Ltd. 00189 *> Sven Hammarling, Numerical Algorithms Group Ltd. 00190 *> \endverbatim 00191 *> 00192 * ===================================================================== 00193 SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) 00194 * 00195 * -- Reference BLAS level3 routine (version 3.4.0) -- 00196 * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- 00197 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00198 * November 2011 00199 * 00200 * .. Scalar Arguments .. 00201 DOUBLE PRECISION ALPHA,BETA 00202 INTEGER K,LDA,LDB,LDC,N 00203 CHARACTER TRANS,UPLO 00204 * .. 00205 * .. Array Arguments .. 00206 DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) 00207 * .. 00208 * 00209 * ===================================================================== 00210 * 00211 * .. External Functions .. 00212 LOGICAL LSAME 00213 EXTERNAL LSAME 00214 * .. 00215 * .. External Subroutines .. 00216 EXTERNAL XERBLA 00217 * .. 00218 * .. Intrinsic Functions .. 00219 INTRINSIC MAX 00220 * .. 00221 * .. Local Scalars .. 00222 DOUBLE PRECISION TEMP1,TEMP2 00223 INTEGER I,INFO,J,L,NROWA 00224 LOGICAL UPPER 00225 * .. 00226 * .. Parameters .. 00227 DOUBLE PRECISION ONE,ZERO 00228 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) 00229 * .. 00230 * 00231 * Test the input parameters. 00232 * 00233 IF (LSAME(TRANS,'N')) THEN 00234 NROWA = N 00235 ELSE 00236 NROWA = K 00237 END IF 00238 UPPER = LSAME(UPLO,'U') 00239 * 00240 INFO = 0 00241 IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN 00242 INFO = 1 00243 ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND. 00244 + (.NOT.LSAME(TRANS,'T')) .AND. 00245 + (.NOT.LSAME(TRANS,'C'))) THEN 00246 INFO = 2 00247 ELSE IF (N.LT.0) THEN 00248 INFO = 3 00249 ELSE IF (K.LT.0) THEN 00250 INFO = 4 00251 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN 00252 INFO = 7 00253 ELSE IF (LDB.LT.MAX(1,NROWA)) THEN 00254 INFO = 9 00255 ELSE IF (LDC.LT.MAX(1,N)) THEN 00256 INFO = 12 00257 END IF 00258 IF (INFO.NE.0) THEN 00259 CALL XERBLA('DSYR2K',INFO) 00260 RETURN 00261 END IF 00262 * 00263 * Quick return if possible. 00264 * 00265 IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR. 00266 + (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN 00267 * 00268 * And when alpha.eq.zero. 00269 * 00270 IF (ALPHA.EQ.ZERO) THEN 00271 IF (UPPER) THEN 00272 IF (BETA.EQ.ZERO) THEN 00273 DO 20 J = 1,N 00274 DO 10 I = 1,J 00275 C(I,J) = ZERO 00276 10 CONTINUE 00277 20 CONTINUE 00278 ELSE 00279 DO 40 J = 1,N 00280 DO 30 I = 1,J 00281 C(I,J) = BETA*C(I,J) 00282 30 CONTINUE 00283 40 CONTINUE 00284 END IF 00285 ELSE 00286 IF (BETA.EQ.ZERO) THEN 00287 DO 60 J = 1,N 00288 DO 50 I = J,N 00289 C(I,J) = ZERO 00290 50 CONTINUE 00291 60 CONTINUE 00292 ELSE 00293 DO 80 J = 1,N 00294 DO 70 I = J,N 00295 C(I,J) = BETA*C(I,J) 00296 70 CONTINUE 00297 80 CONTINUE 00298 END IF 00299 END IF 00300 RETURN 00301 END IF 00302 * 00303 * Start the operations. 00304 * 00305 IF (LSAME(TRANS,'N')) THEN 00306 * 00307 * Form C := alpha*A*B**T + alpha*B*A**T + C. 00308 * 00309 IF (UPPER) THEN 00310 DO 130 J = 1,N 00311 IF (BETA.EQ.ZERO) THEN 00312 DO 90 I = 1,J 00313 C(I,J) = ZERO 00314 90 CONTINUE 00315 ELSE IF (BETA.NE.ONE) THEN 00316 DO 100 I = 1,J 00317 C(I,J) = BETA*C(I,J) 00318 100 CONTINUE 00319 END IF 00320 DO 120 L = 1,K 00321 IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN 00322 TEMP1 = ALPHA*B(J,L) 00323 TEMP2 = ALPHA*A(J,L) 00324 DO 110 I = 1,J 00325 C(I,J) = C(I,J) + A(I,L)*TEMP1 + 00326 + B(I,L)*TEMP2 00327 110 CONTINUE 00328 END IF 00329 120 CONTINUE 00330 130 CONTINUE 00331 ELSE 00332 DO 180 J = 1,N 00333 IF (BETA.EQ.ZERO) THEN 00334 DO 140 I = J,N 00335 C(I,J) = ZERO 00336 140 CONTINUE 00337 ELSE IF (BETA.NE.ONE) THEN 00338 DO 150 I = J,N 00339 C(I,J) = BETA*C(I,J) 00340 150 CONTINUE 00341 END IF 00342 DO 170 L = 1,K 00343 IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN 00344 TEMP1 = ALPHA*B(J,L) 00345 TEMP2 = ALPHA*A(J,L) 00346 DO 160 I = J,N 00347 C(I,J) = C(I,J) + A(I,L)*TEMP1 + 00348 + B(I,L)*TEMP2 00349 160 CONTINUE 00350 END IF 00351 170 CONTINUE 00352 180 CONTINUE 00353 END IF 00354 ELSE 00355 * 00356 * Form C := alpha*A**T*B + alpha*B**T*A + C. 00357 * 00358 IF (UPPER) THEN 00359 DO 210 J = 1,N 00360 DO 200 I = 1,J 00361 TEMP1 = ZERO 00362 TEMP2 = ZERO 00363 DO 190 L = 1,K 00364 TEMP1 = TEMP1 + A(L,I)*B(L,J) 00365 TEMP2 = TEMP2 + B(L,I)*A(L,J) 00366 190 CONTINUE 00367 IF (BETA.EQ.ZERO) THEN 00368 C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2 00369 ELSE 00370 C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 + 00371 + ALPHA*TEMP2 00372 END IF 00373 200 CONTINUE 00374 210 CONTINUE 00375 ELSE 00376 DO 240 J = 1,N 00377 DO 230 I = J,N 00378 TEMP1 = ZERO 00379 TEMP2 = ZERO 00380 DO 220 L = 1,K 00381 TEMP1 = TEMP1 + A(L,I)*B(L,J) 00382 TEMP2 = TEMP2 + B(L,I)*A(L,J) 00383 220 CONTINUE 00384 IF (BETA.EQ.ZERO) THEN 00385 C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2 00386 ELSE 00387 C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 + 00388 + ALPHA*TEMP2 00389 END IF 00390 230 CONTINUE 00391 240 CONTINUE 00392 END IF 00393 END IF 00394 * 00395 RETURN 00396 * 00397 * End of DSYR2K. 00398 * 00399 END