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