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