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