|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) 00002 * .. Scalar Arguments .. 00003 DOUBLE COMPLEX ALPHA 00004 DOUBLE PRECISION BETA 00005 INTEGER K,LDA,LDB,LDC,N 00006 CHARACTER TRANS,UPLO 00007 * .. 00008 * .. Array Arguments .. 00009 DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) 00010 * .. 00011 * 00012 * Purpose 00013 * ======= 00014 * 00015 * ZHER2K performs one of the hermitian rank 2k operations 00016 * 00017 * C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + beta*C, 00018 * 00019 * or 00020 * 00021 * C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + beta*C, 00022 * 00023 * where alpha and beta are scalars with beta real, C is an n by n 00024 * hermitian matrix and A and B are n by k matrices in the first case 00025 * and k by n matrices in the second case. 00026 * 00027 * Arguments 00028 * ========== 00029 * 00030 * UPLO - CHARACTER*1. 00031 * On entry, UPLO specifies whether the upper or lower 00032 * triangular part of the array C is to be referenced as 00033 * follows: 00034 * 00035 * UPLO = 'U' or 'u' Only the upper triangular part of C 00036 * is to be referenced. 00037 * 00038 * UPLO = 'L' or 'l' Only the lower triangular part of C 00039 * is to be referenced. 00040 * 00041 * Unchanged on exit. 00042 * 00043 * TRANS - CHARACTER*1. 00044 * On entry, TRANS specifies the operation to be performed as 00045 * follows: 00046 * 00047 * TRANS = 'N' or 'n' C := alpha*A*conjg( B' ) + 00048 * conjg( alpha )*B*conjg( A' ) + 00049 * beta*C. 00050 * 00051 * TRANS = 'C' or 'c' C := alpha*conjg( A' )*B + 00052 * conjg( alpha )*conjg( 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 = 'C' or 'c', K specifies the number of rows of the 00066 * matrices A and B. K must be at least zero. 00067 * Unchanged on exit. 00068 * 00069 * ALPHA - COMPLEX*16 . 00070 * On entry, ALPHA specifies the scalar alpha. 00071 * Unchanged on exit. 00072 * 00073 * A - COMPLEX*16 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 - COMPLEX*16 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 - COMPLEX*16 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 hermitian 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 hermitian 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 * Note that the imaginary parts of the diagonal elements need 00121 * not be set, they are assumed to be zero, and on exit they 00122 * are set to zero. 00123 * 00124 * LDC - INTEGER. 00125 * On entry, LDC specifies the first dimension of C as declared 00126 * in the calling (sub) program. LDC must be at least 00127 * max( 1, n ). 00128 * Unchanged on exit. 00129 * 00130 * Further Details 00131 * =============== 00132 * 00133 * Level 3 Blas routine. 00134 * 00135 * -- Written on 8-February-1989. 00136 * Jack Dongarra, Argonne National Laboratory. 00137 * Iain Duff, AERE Harwell. 00138 * Jeremy Du Croz, Numerical Algorithms Group Ltd. 00139 * Sven Hammarling, Numerical Algorithms Group Ltd. 00140 * 00141 * -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1. 00142 * Ed Anderson, Cray Research Inc. 00143 * 00144 * ===================================================================== 00145 * 00146 * .. External Functions .. 00147 LOGICAL LSAME 00148 EXTERNAL LSAME 00149 * .. 00150 * .. External Subroutines .. 00151 EXTERNAL XERBLA 00152 * .. 00153 * .. Intrinsic Functions .. 00154 INTRINSIC DBLE,DCONJG,MAX 00155 * .. 00156 * .. Local Scalars .. 00157 DOUBLE COMPLEX TEMP1,TEMP2 00158 INTEGER I,INFO,J,L,NROWA 00159 LOGICAL UPPER 00160 * .. 00161 * .. Parameters .. 00162 DOUBLE PRECISION ONE 00163 PARAMETER (ONE=1.0D+0) 00164 DOUBLE COMPLEX ZERO 00165 PARAMETER (ZERO= (0.0D+0,0.0D+0)) 00166 * .. 00167 * 00168 * Test the input parameters. 00169 * 00170 IF (LSAME(TRANS,'N')) THEN 00171 NROWA = N 00172 ELSE 00173 NROWA = K 00174 END IF 00175 UPPER = LSAME(UPLO,'U') 00176 * 00177 INFO = 0 00178 IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN 00179 INFO = 1 00180 ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND. 00181 + (.NOT.LSAME(TRANS,'C'))) THEN 00182 INFO = 2 00183 ELSE IF (N.LT.0) THEN 00184 INFO = 3 00185 ELSE IF (K.LT.0) THEN 00186 INFO = 4 00187 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN 00188 INFO = 7 00189 ELSE IF (LDB.LT.MAX(1,NROWA)) THEN 00190 INFO = 9 00191 ELSE IF (LDC.LT.MAX(1,N)) THEN 00192 INFO = 12 00193 END IF 00194 IF (INFO.NE.0) THEN 00195 CALL XERBLA('ZHER2K',INFO) 00196 RETURN 00197 END IF 00198 * 00199 * Quick return if possible. 00200 * 00201 IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR. 00202 + (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN 00203 * 00204 * And when alpha.eq.zero. 00205 * 00206 IF (ALPHA.EQ.ZERO) THEN 00207 IF (UPPER) THEN 00208 IF (BETA.EQ.DBLE(ZERO)) THEN 00209 DO 20 J = 1,N 00210 DO 10 I = 1,J 00211 C(I,J) = ZERO 00212 10 CONTINUE 00213 20 CONTINUE 00214 ELSE 00215 DO 40 J = 1,N 00216 DO 30 I = 1,J - 1 00217 C(I,J) = BETA*C(I,J) 00218 30 CONTINUE 00219 C(J,J) = BETA*DBLE(C(J,J)) 00220 40 CONTINUE 00221 END IF 00222 ELSE 00223 IF (BETA.EQ.DBLE(ZERO)) THEN 00224 DO 60 J = 1,N 00225 DO 50 I = J,N 00226 C(I,J) = ZERO 00227 50 CONTINUE 00228 60 CONTINUE 00229 ELSE 00230 DO 80 J = 1,N 00231 C(J,J) = BETA*DBLE(C(J,J)) 00232 DO 70 I = J + 1,N 00233 C(I,J) = BETA*C(I,J) 00234 70 CONTINUE 00235 80 CONTINUE 00236 END IF 00237 END IF 00238 RETURN 00239 END IF 00240 * 00241 * Start the operations. 00242 * 00243 IF (LSAME(TRANS,'N')) THEN 00244 * 00245 * Form C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + 00246 * C. 00247 * 00248 IF (UPPER) THEN 00249 DO 130 J = 1,N 00250 IF (BETA.EQ.DBLE(ZERO)) THEN 00251 DO 90 I = 1,J 00252 C(I,J) = ZERO 00253 90 CONTINUE 00254 ELSE IF (BETA.NE.ONE) THEN 00255 DO 100 I = 1,J - 1 00256 C(I,J) = BETA*C(I,J) 00257 100 CONTINUE 00258 C(J,J) = BETA*DBLE(C(J,J)) 00259 ELSE 00260 C(J,J) = DBLE(C(J,J)) 00261 END IF 00262 DO 120 L = 1,K 00263 IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN 00264 TEMP1 = ALPHA*DCONJG(B(J,L)) 00265 TEMP2 = DCONJG(ALPHA*A(J,L)) 00266 DO 110 I = 1,J - 1 00267 C(I,J) = C(I,J) + A(I,L)*TEMP1 + 00268 + B(I,L)*TEMP2 00269 110 CONTINUE 00270 C(J,J) = DBLE(C(J,J)) + 00271 + DBLE(A(J,L)*TEMP1+B(J,L)*TEMP2) 00272 END IF 00273 120 CONTINUE 00274 130 CONTINUE 00275 ELSE 00276 DO 180 J = 1,N 00277 IF (BETA.EQ.DBLE(ZERO)) THEN 00278 DO 140 I = J,N 00279 C(I,J) = ZERO 00280 140 CONTINUE 00281 ELSE IF (BETA.NE.ONE) THEN 00282 DO 150 I = J + 1,N 00283 C(I,J) = BETA*C(I,J) 00284 150 CONTINUE 00285 C(J,J) = BETA*DBLE(C(J,J)) 00286 ELSE 00287 C(J,J) = DBLE(C(J,J)) 00288 END IF 00289 DO 170 L = 1,K 00290 IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN 00291 TEMP1 = ALPHA*DCONJG(B(J,L)) 00292 TEMP2 = DCONJG(ALPHA*A(J,L)) 00293 DO 160 I = J + 1,N 00294 C(I,J) = C(I,J) + A(I,L)*TEMP1 + 00295 + B(I,L)*TEMP2 00296 160 CONTINUE 00297 C(J,J) = DBLE(C(J,J)) + 00298 + DBLE(A(J,L)*TEMP1+B(J,L)*TEMP2) 00299 END IF 00300 170 CONTINUE 00301 180 CONTINUE 00302 END IF 00303 ELSE 00304 * 00305 * Form C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + 00306 * C. 00307 * 00308 IF (UPPER) THEN 00309 DO 210 J = 1,N 00310 DO 200 I = 1,J 00311 TEMP1 = ZERO 00312 TEMP2 = ZERO 00313 DO 190 L = 1,K 00314 TEMP1 = TEMP1 + DCONJG(A(L,I))*B(L,J) 00315 TEMP2 = TEMP2 + DCONJG(B(L,I))*A(L,J) 00316 190 CONTINUE 00317 IF (I.EQ.J) THEN 00318 IF (BETA.EQ.DBLE(ZERO)) THEN 00319 C(J,J) = DBLE(ALPHA*TEMP1+ 00320 + DCONJG(ALPHA)*TEMP2) 00321 ELSE 00322 C(J,J) = BETA*DBLE(C(J,J)) + 00323 + DBLE(ALPHA*TEMP1+ 00324 + DCONJG(ALPHA)*TEMP2) 00325 END IF 00326 ELSE 00327 IF (BETA.EQ.DBLE(ZERO)) THEN 00328 C(I,J) = ALPHA*TEMP1 + DCONJG(ALPHA)*TEMP2 00329 ELSE 00330 C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 + 00331 + DCONJG(ALPHA)*TEMP2 00332 END IF 00333 END IF 00334 200 CONTINUE 00335 210 CONTINUE 00336 ELSE 00337 DO 240 J = 1,N 00338 DO 230 I = J,N 00339 TEMP1 = ZERO 00340 TEMP2 = ZERO 00341 DO 220 L = 1,K 00342 TEMP1 = TEMP1 + DCONJG(A(L,I))*B(L,J) 00343 TEMP2 = TEMP2 + DCONJG(B(L,I))*A(L,J) 00344 220 CONTINUE 00345 IF (I.EQ.J) THEN 00346 IF (BETA.EQ.DBLE(ZERO)) THEN 00347 C(J,J) = DBLE(ALPHA*TEMP1+ 00348 + DCONJG(ALPHA)*TEMP2) 00349 ELSE 00350 C(J,J) = BETA*DBLE(C(J,J)) + 00351 + DBLE(ALPHA*TEMP1+ 00352 + DCONJG(ALPHA)*TEMP2) 00353 END IF 00354 ELSE 00355 IF (BETA.EQ.DBLE(ZERO)) THEN 00356 C(I,J) = ALPHA*TEMP1 + DCONJG(ALPHA)*TEMP2 00357 ELSE 00358 C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 + 00359 + DCONJG(ALPHA)*TEMP2 00360 END IF 00361 END IF 00362 230 CONTINUE 00363 240 CONTINUE 00364 END IF 00365 END IF 00366 * 00367 RETURN 00368 * 00369 * End of ZHER2K. 00370 * 00371 END
1.7.3 
