|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) 00002 * .. Scalar Arguments .. 00003 DOUBLE PRECISION ALPHA,BETA 00004 INTEGER K,LDA,LDB,LDC,N 00005 CHARACTER TRANS,UPLO 00006 * .. 00007 * .. Array Arguments .. 00008 DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) 00009 * .. 00010 * 00011 * Purpose 00012 * ======= 00013 * 00014 * DSYR2K performs one of the symmetric rank 2k operations 00015 * 00016 * C := alpha*A*B' + alpha*B*A' + beta*C, 00017 * 00018 * or 00019 * 00020 * C := alpha*A'*B + alpha*B'*A + beta*C, 00021 * 00022 * where alpha and beta are scalars, C is an n by n symmetric matrix 00023 * and A and B are n by k matrices in the first case and k by n 00024 * matrices in the second case. 00025 * 00026 * Arguments 00027 * ========== 00028 * 00029 * UPLO - CHARACTER*1. 00030 * On entry, UPLO specifies whether the upper or lower 00031 * triangular part of the array C is to be referenced as 00032 * follows: 00033 * 00034 * UPLO = 'U' or 'u' Only the upper triangular part of C 00035 * is to be referenced. 00036 * 00037 * UPLO = 'L' or 'l' Only the lower triangular part of C 00038 * is to be referenced. 00039 * 00040 * Unchanged on exit. 00041 * 00042 * TRANS - CHARACTER*1. 00043 * On entry, TRANS specifies the operation to be performed as 00044 * follows: 00045 * 00046 * TRANS = 'N' or 'n' C := alpha*A*B' + alpha*B*A' + 00047 * beta*C. 00048 * 00049 * TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A + 00050 * beta*C. 00051 * 00052 * TRANS = 'C' or 'c' C := alpha*A'*B + alpha*B'*A + 00053 * beta*C. 00054 * 00055 * Unchanged on exit. 00056 * 00057 * N - INTEGER. 00058 * On entry, N specifies the order of the matrix C. N must be 00059 * at least zero. 00060 * Unchanged on exit. 00061 * 00062 * K - INTEGER. 00063 * On entry with TRANS = 'N' or 'n', K specifies the number 00064 * of columns of the matrices A and B, and on entry with 00065 * TRANS = 'T' or 't' or 'C' or 'c', K specifies the number 00066 * of rows of the matrices A and B. K must be at least zero. 00067 * Unchanged on exit. 00068 * 00069 * ALPHA - DOUBLE PRECISION. 00070 * On entry, ALPHA specifies the scalar alpha. 00071 * Unchanged on exit. 00072 * 00073 * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is 00074 * k when TRANS = 'N' or 'n', and is n otherwise. 00075 * Before entry with TRANS = 'N' or 'n', the leading n by k 00076 * part of the array A must contain the matrix A, otherwise 00077 * the leading k by n part of the array A must contain the 00078 * matrix A. 00079 * Unchanged on exit. 00080 * 00081 * LDA - INTEGER. 00082 * On entry, LDA specifies the first dimension of A as declared 00083 * in the calling (sub) program. When TRANS = 'N' or 'n' 00084 * then LDA must be at least max( 1, n ), otherwise LDA must 00085 * be at least max( 1, k ). 00086 * Unchanged on exit. 00087 * 00088 * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is 00089 * k when TRANS = 'N' or 'n', and is n otherwise. 00090 * Before entry with TRANS = 'N' or 'n', the leading n by k 00091 * part of the array B must contain the matrix B, otherwise 00092 * the leading k by n part of the array B must contain the 00093 * matrix B. 00094 * Unchanged on exit. 00095 * 00096 * LDB - INTEGER. 00097 * On entry, LDB specifies the first dimension of B as declared 00098 * in the calling (sub) program. When TRANS = 'N' or 'n' 00099 * then LDB must be at least max( 1, n ), otherwise LDB must 00100 * be at least max( 1, k ). 00101 * Unchanged on exit. 00102 * 00103 * BETA - DOUBLE PRECISION. 00104 * On entry, BETA specifies the scalar beta. 00105 * Unchanged on exit. 00106 * 00107 * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). 00108 * Before entry with UPLO = 'U' or 'u', the leading n by n 00109 * upper triangular part of the array C must contain the upper 00110 * triangular part of the symmetric matrix and the strictly 00111 * lower triangular part of C is not referenced. On exit, the 00112 * upper triangular part of the array C is overwritten by the 00113 * upper triangular part of the updated matrix. 00114 * Before entry with UPLO = 'L' or 'l', the leading n by n 00115 * lower triangular part of the array C must contain the lower 00116 * triangular part of the symmetric matrix and the strictly 00117 * upper triangular part of C is not referenced. On exit, the 00118 * lower triangular part of the array C is overwritten by the 00119 * lower triangular part of the updated matrix. 00120 * 00121 * LDC - INTEGER. 00122 * On entry, LDC specifies the first dimension of C as declared 00123 * in the calling (sub) program. LDC must be at least 00124 * max( 1, n ). 00125 * Unchanged on exit. 00126 * 00127 * Further Details 00128 * =============== 00129 * 00130 * Level 3 Blas routine. 00131 * 00132 * 00133 * -- Written on 8-February-1989. 00134 * Jack Dongarra, Argonne National Laboratory. 00135 * Iain Duff, AERE Harwell. 00136 * Jeremy Du Croz, Numerical Algorithms Group Ltd. 00137 * Sven Hammarling, Numerical Algorithms Group Ltd. 00138 * 00139 * ===================================================================== 00140 * 00141 * .. External Functions .. 00142 LOGICAL LSAME 00143 EXTERNAL LSAME 00144 * .. 00145 * .. External Subroutines .. 00146 EXTERNAL XERBLA 00147 * .. 00148 * .. Intrinsic Functions .. 00149 INTRINSIC MAX 00150 * .. 00151 * .. Local Scalars .. 00152 DOUBLE PRECISION TEMP1,TEMP2 00153 INTEGER I,INFO,J,L,NROWA 00154 LOGICAL UPPER 00155 * .. 00156 * .. Parameters .. 00157 DOUBLE PRECISION ONE,ZERO 00158 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) 00159 * .. 00160 * 00161 * Test the input parameters. 00162 * 00163 IF (LSAME(TRANS,'N')) THEN 00164 NROWA = N 00165 ELSE 00166 NROWA = K 00167 END IF 00168 UPPER = LSAME(UPLO,'U') 00169 * 00170 INFO = 0 00171 IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN 00172 INFO = 1 00173 ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND. 00174 + (.NOT.LSAME(TRANS,'T')) .AND. 00175 + (.NOT.LSAME(TRANS,'C'))) THEN 00176 INFO = 2 00177 ELSE IF (N.LT.0) THEN 00178 INFO = 3 00179 ELSE IF (K.LT.0) THEN 00180 INFO = 4 00181 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN 00182 INFO = 7 00183 ELSE IF (LDB.LT.MAX(1,NROWA)) THEN 00184 INFO = 9 00185 ELSE IF (LDC.LT.MAX(1,N)) THEN 00186 INFO = 12 00187 END IF 00188 IF (INFO.NE.0) THEN 00189 CALL XERBLA('DSYR2K',INFO) 00190 RETURN 00191 END IF 00192 * 00193 * Quick return if possible. 00194 * 00195 IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR. 00196 + (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN 00197 * 00198 * And when alpha.eq.zero. 00199 * 00200 IF (ALPHA.EQ.ZERO) THEN 00201 IF (UPPER) THEN 00202 IF (BETA.EQ.ZERO) THEN 00203 DO 20 J = 1,N 00204 DO 10 I = 1,J 00205 C(I,J) = ZERO 00206 10 CONTINUE 00207 20 CONTINUE 00208 ELSE 00209 DO 40 J = 1,N 00210 DO 30 I = 1,J 00211 C(I,J) = BETA*C(I,J) 00212 30 CONTINUE 00213 40 CONTINUE 00214 END IF 00215 ELSE 00216 IF (BETA.EQ.ZERO) THEN 00217 DO 60 J = 1,N 00218 DO 50 I = J,N 00219 C(I,J) = ZERO 00220 50 CONTINUE 00221 60 CONTINUE 00222 ELSE 00223 DO 80 J = 1,N 00224 DO 70 I = J,N 00225 C(I,J) = BETA*C(I,J) 00226 70 CONTINUE 00227 80 CONTINUE 00228 END IF 00229 END IF 00230 RETURN 00231 END IF 00232 * 00233 * Start the operations. 00234 * 00235 IF (LSAME(TRANS,'N')) THEN 00236 * 00237 * Form C := alpha*A*B' + alpha*B*A' + C. 00238 * 00239 IF (UPPER) THEN 00240 DO 130 J = 1,N 00241 IF (BETA.EQ.ZERO) THEN 00242 DO 90 I = 1,J 00243 C(I,J) = ZERO 00244 90 CONTINUE 00245 ELSE IF (BETA.NE.ONE) THEN 00246 DO 100 I = 1,J 00247 C(I,J) = BETA*C(I,J) 00248 100 CONTINUE 00249 END IF 00250 DO 120 L = 1,K 00251 IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN 00252 TEMP1 = ALPHA*B(J,L) 00253 TEMP2 = ALPHA*A(J,L) 00254 DO 110 I = 1,J 00255 C(I,J) = C(I,J) + A(I,L)*TEMP1 + 00256 + B(I,L)*TEMP2 00257 110 CONTINUE 00258 END IF 00259 120 CONTINUE 00260 130 CONTINUE 00261 ELSE 00262 DO 180 J = 1,N 00263 IF (BETA.EQ.ZERO) THEN 00264 DO 140 I = J,N 00265 C(I,J) = ZERO 00266 140 CONTINUE 00267 ELSE IF (BETA.NE.ONE) THEN 00268 DO 150 I = J,N 00269 C(I,J) = BETA*C(I,J) 00270 150 CONTINUE 00271 END IF 00272 DO 170 L = 1,K 00273 IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN 00274 TEMP1 = ALPHA*B(J,L) 00275 TEMP2 = ALPHA*A(J,L) 00276 DO 160 I = J,N 00277 C(I,J) = C(I,J) + A(I,L)*TEMP1 + 00278 + B(I,L)*TEMP2 00279 160 CONTINUE 00280 END IF 00281 170 CONTINUE 00282 180 CONTINUE 00283 END IF 00284 ELSE 00285 * 00286 * Form C := alpha*A'*B + alpha*B'*A + C. 00287 * 00288 IF (UPPER) THEN 00289 DO 210 J = 1,N 00290 DO 200 I = 1,J 00291 TEMP1 = ZERO 00292 TEMP2 = ZERO 00293 DO 190 L = 1,K 00294 TEMP1 = TEMP1 + A(L,I)*B(L,J) 00295 TEMP2 = TEMP2 + B(L,I)*A(L,J) 00296 190 CONTINUE 00297 IF (BETA.EQ.ZERO) THEN 00298 C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2 00299 ELSE 00300 C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 + 00301 + ALPHA*TEMP2 00302 END IF 00303 200 CONTINUE 00304 210 CONTINUE 00305 ELSE 00306 DO 240 J = 1,N 00307 DO 230 I = J,N 00308 TEMP1 = ZERO 00309 TEMP2 = ZERO 00310 DO 220 L = 1,K 00311 TEMP1 = TEMP1 + A(L,I)*B(L,J) 00312 TEMP2 = TEMP2 + B(L,I)*A(L,J) 00313 220 CONTINUE 00314 IF (BETA.EQ.ZERO) THEN 00315 C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2 00316 ELSE 00317 C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 + 00318 + ALPHA*TEMP2 00319 END IF 00320 230 CONTINUE 00321 240 CONTINUE 00322 END IF 00323 END IF 00324 * 00325 RETURN 00326 * 00327 * End of DSYR2K. 00328 * 00329 END
1.7.3 
