1: SUBROUTINE ZHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
2: * .. Scalar Arguments ..
3: DOUBLE COMPLEX ALPHA
4: DOUBLE PRECISION BETA
5: INTEGER K,LDA,LDB,LDC,N
6: CHARACTER TRANS,UPLO
7: * ..
8: * .. Array Arguments ..
9: DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
10: * ..
11: *
12: * Purpose
13: * =======
14: *
15: * ZHER2K performs one of the hermitian rank 2k operations
16: *
17: * C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + beta*C,
18: *
19: * or
20: *
21: * C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + beta*C,
22: *
23: * where alpha and beta are scalars with beta real, C is an n by n
24: * hermitian matrix and A and B are n by k matrices in the first case
25: * and k by n matrices in the second case.
26: *
27: * Arguments
28: * ==========
29: *
30: * UPLO - CHARACTER*1.
31: * On entry, UPLO specifies whether the upper or lower
32: * triangular part of the array C is to be referenced as
33: * follows:
34: *
35: * UPLO = 'U' or 'u' Only the upper triangular part of C
36: * is to be referenced.
37: *
38: * UPLO = 'L' or 'l' Only the lower triangular part of C
39: * is to be referenced.
40: *
41: * Unchanged on exit.
42: *
43: * TRANS - CHARACTER*1.
44: * On entry, TRANS specifies the operation to be performed as
45: * follows:
46: *
47: * TRANS = 'N' or 'n' C := alpha*A*conjg( B' ) +
48: * conjg( alpha )*B*conjg( A' ) +
49: * beta*C.
50: *
51: * TRANS = 'C' or 'c' C := alpha*conjg( A' )*B +
52: * conjg( alpha )*conjg( B' )*A +
53: * beta*C.
54: *
55: * Unchanged on exit.
56: *
57: * N - INTEGER.
58: * On entry, N specifies the order of the matrix C. N must be
59: * at least zero.
60: * Unchanged on exit.
61: *
62: * K - INTEGER.
63: * On entry with TRANS = 'N' or 'n', K specifies the number
64: * of columns of the matrices A and B, and on entry with
65: * TRANS = 'C' or 'c', K specifies the number of rows of the
66: * matrices A and B. K must be at least zero.
67: * Unchanged on exit.
68: *
69: * ALPHA - COMPLEX*16 .
70: * On entry, ALPHA specifies the scalar alpha.
71: * Unchanged on exit.
72: *
73: * A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is
74: * k when TRANS = 'N' or 'n', and is n otherwise.
75: * Before entry with TRANS = 'N' or 'n', the leading n by k
76: * part of the array A must contain the matrix A, otherwise
77: * the leading k by n part of the array A must contain the
78: * matrix A.
79: * Unchanged on exit.
80: *
81: * LDA - INTEGER.
82: * On entry, LDA specifies the first dimension of A as declared
83: * in the calling (sub) program. When TRANS = 'N' or 'n'
84: * then LDA must be at least max( 1, n ), otherwise LDA must
85: * be at least max( 1, k ).
86: * Unchanged on exit.
87: *
88: * B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is
89: * k when TRANS = 'N' or 'n', and is n otherwise.
90: * Before entry with TRANS = 'N' or 'n', the leading n by k
91: * part of the array B must contain the matrix B, otherwise
92: * the leading k by n part of the array B must contain the
93: * matrix B.
94: * Unchanged on exit.
95: *
96: * LDB - INTEGER.
97: * On entry, LDB specifies the first dimension of B as declared
98: * in the calling (sub) program. When TRANS = 'N' or 'n'
99: * then LDB must be at least max( 1, n ), otherwise LDB must
100: * be at least max( 1, k ).
101: * Unchanged on exit.
102: *
103: * BETA - DOUBLE PRECISION .
104: * On entry, BETA specifies the scalar beta.
105: * Unchanged on exit.
106: *
107: * C - COMPLEX*16 array of DIMENSION ( LDC, n ).
108: * Before entry with UPLO = 'U' or 'u', the leading n by n
109: * upper triangular part of the array C must contain the upper
110: * triangular part of the hermitian matrix and the strictly
111: * lower triangular part of C is not referenced. On exit, the
112: * upper triangular part of the array C is overwritten by the
113: * upper triangular part of the updated matrix.
114: * Before entry with UPLO = 'L' or 'l', the leading n by n
115: * lower triangular part of the array C must contain the lower
116: * triangular part of the hermitian matrix and the strictly
117: * upper triangular part of C is not referenced. On exit, the
118: * lower triangular part of the array C is overwritten by the
119: * lower triangular part of the updated matrix.
120: * Note that the imaginary parts of the diagonal elements need
121: * not be set, they are assumed to be zero, and on exit they
122: * are set to zero.
123: *
124: * LDC - INTEGER.
125: * On entry, LDC specifies the first dimension of C as declared
126: * in the calling (sub) program. LDC must be at least
127: * max( 1, n ).
128: * Unchanged on exit.
129: *
130: * Further Details
131: * ===============
132: *
133: * Level 3 Blas routine.
134: *
135: * -- Written on 8-February-1989.
136: * Jack Dongarra, Argonne National Laboratory.
137: * Iain Duff, AERE Harwell.
138: * Jeremy Du Croz, Numerical Algorithms Group Ltd.
139: * Sven Hammarling, Numerical Algorithms Group Ltd.
140: *
141: * -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1.
142: * Ed Anderson, Cray Research Inc.
143: *
144: * =====================================================================
145: *
146: * .. External Functions ..
147: LOGICAL LSAME
148: EXTERNAL LSAME
149: * ..
150: * .. External Subroutines ..
151: EXTERNAL XERBLA
152: * ..
153: * .. Intrinsic Functions ..
154: INTRINSIC DBLE,DCONJG,MAX
155: * ..
156: * .. Local Scalars ..
157: DOUBLE COMPLEX TEMP1,TEMP2
158: INTEGER I,INFO,J,L,NROWA
159: LOGICAL UPPER
160: * ..
161: * .. Parameters ..
162: DOUBLE PRECISION ONE
163: PARAMETER (ONE=1.0D+0)
164: DOUBLE COMPLEX ZERO
165: PARAMETER (ZERO= (0.0D+0,0.0D+0))
166: * ..
167: *
168: * Test the input parameters.
169: *
170: IF (LSAME(TRANS,'N')) THEN
171: NROWA = N
172: ELSE
173: NROWA = K
174: END IF
175: UPPER = LSAME(UPLO,'U')
176: *
177: INFO = 0
178: IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
179: INFO = 1
180: ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
181: + (.NOT.LSAME(TRANS,'C'))) THEN
182: INFO = 2
183: ELSE IF (N.LT.0) THEN
184: INFO = 3
185: ELSE IF (K.LT.0) THEN
186: INFO = 4
187: ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
188: INFO = 7
189: ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
190: INFO = 9
191: ELSE IF (LDC.LT.MAX(1,N)) THEN
192: INFO = 12
193: END IF
194: IF (INFO.NE.0) THEN
195: CALL XERBLA('ZHER2K',INFO)
196: RETURN
197: END IF
198: *
199: * Quick return if possible.
200: *
201: IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
202: + (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
203: *
204: * And when alpha.eq.zero.
205: *
206: IF (ALPHA.EQ.ZERO) THEN
207: IF (UPPER) THEN
208: IF (BETA.EQ.DBLE(ZERO)) THEN
209: DO 20 J = 1,N
210: DO 10 I = 1,J
211: C(I,J) = ZERO
212: 10 CONTINUE
213: 20 CONTINUE
214: ELSE
215: DO 40 J = 1,N
216: DO 30 I = 1,J - 1
217: C(I,J) = BETA*C(I,J)
218: 30 CONTINUE
219: C(J,J) = BETA*DBLE(C(J,J))
220: 40 CONTINUE
221: END IF
222: ELSE
223: IF (BETA.EQ.DBLE(ZERO)) THEN
224: DO 60 J = 1,N
225: DO 50 I = J,N
226: C(I,J) = ZERO
227: 50 CONTINUE
228: 60 CONTINUE
229: ELSE
230: DO 80 J = 1,N
231: C(J,J) = BETA*DBLE(C(J,J))
232: DO 70 I = J + 1,N
233: C(I,J) = BETA*C(I,J)
234: 70 CONTINUE
235: 80 CONTINUE
236: END IF
237: END IF
238: RETURN
239: END IF
240: *
241: * Start the operations.
242: *
243: IF (LSAME(TRANS,'N')) THEN
244: *
245: * Form C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) +
246: * C.
247: *
248: IF (UPPER) THEN
249: DO 130 J = 1,N
250: IF (BETA.EQ.DBLE(ZERO)) THEN
251: DO 90 I = 1,J
252: C(I,J) = ZERO
253: 90 CONTINUE
254: ELSE IF (BETA.NE.ONE) THEN
255: DO 100 I = 1,J - 1
256: C(I,J) = BETA*C(I,J)
257: 100 CONTINUE
258: C(J,J) = BETA*DBLE(C(J,J))
259: ELSE
260: C(J,J) = DBLE(C(J,J))
261: END IF
262: DO 120 L = 1,K
263: IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
264: TEMP1 = ALPHA*DCONJG(B(J,L))
265: TEMP2 = DCONJG(ALPHA*A(J,L))
266: DO 110 I = 1,J - 1
267: C(I,J) = C(I,J) + A(I,L)*TEMP1 +
268: + B(I,L)*TEMP2
269: 110 CONTINUE
270: C(J,J) = DBLE(C(J,J)) +
271: + DBLE(A(J,L)*TEMP1+B(J,L)*TEMP2)
272: END IF
273: 120 CONTINUE
274: 130 CONTINUE
275: ELSE
276: DO 180 J = 1,N
277: IF (BETA.EQ.DBLE(ZERO)) THEN
278: DO 140 I = J,N
279: C(I,J) = ZERO
280: 140 CONTINUE
281: ELSE IF (BETA.NE.ONE) THEN
282: DO 150 I = J + 1,N
283: C(I,J) = BETA*C(I,J)
284: 150 CONTINUE
285: C(J,J) = BETA*DBLE(C(J,J))
286: ELSE
287: C(J,J) = DBLE(C(J,J))
288: END IF
289: DO 170 L = 1,K
290: IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
291: TEMP1 = ALPHA*DCONJG(B(J,L))
292: TEMP2 = DCONJG(ALPHA*A(J,L))
293: DO 160 I = J + 1,N
294: C(I,J) = C(I,J) + A(I,L)*TEMP1 +
295: + B(I,L)*TEMP2
296: 160 CONTINUE
297: C(J,J) = DBLE(C(J,J)) +
298: + DBLE(A(J,L)*TEMP1+B(J,L)*TEMP2)
299: END IF
300: 170 CONTINUE
301: 180 CONTINUE
302: END IF
303: ELSE
304: *
305: * Form C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A +
306: * C.
307: *
308: IF (UPPER) THEN
309: DO 210 J = 1,N
310: DO 200 I = 1,J
311: TEMP1 = ZERO
312: TEMP2 = ZERO
313: DO 190 L = 1,K
314: TEMP1 = TEMP1 + DCONJG(A(L,I))*B(L,J)
315: TEMP2 = TEMP2 + DCONJG(B(L,I))*A(L,J)
316: 190 CONTINUE
317: IF (I.EQ.J) THEN
318: IF (BETA.EQ.DBLE(ZERO)) THEN
319: C(J,J) = DBLE(ALPHA*TEMP1+
320: + DCONJG(ALPHA)*TEMP2)
321: ELSE
322: C(J,J) = BETA*DBLE(C(J,J)) +
323: + DBLE(ALPHA*TEMP1+
324: + DCONJG(ALPHA)*TEMP2)
325: END IF
326: ELSE
327: IF (BETA.EQ.DBLE(ZERO)) THEN
328: C(I,J) = ALPHA*TEMP1 + DCONJG(ALPHA)*TEMP2
329: ELSE
330: C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
331: + DCONJG(ALPHA)*TEMP2
332: END IF
333: END IF
334: 200 CONTINUE
335: 210 CONTINUE
336: ELSE
337: DO 240 J = 1,N
338: DO 230 I = J,N
339: TEMP1 = ZERO
340: TEMP2 = ZERO
341: DO 220 L = 1,K
342: TEMP1 = TEMP1 + DCONJG(A(L,I))*B(L,J)
343: TEMP2 = TEMP2 + DCONJG(B(L,I))*A(L,J)
344: 220 CONTINUE
345: IF (I.EQ.J) THEN
346: IF (BETA.EQ.DBLE(ZERO)) THEN
347: C(J,J) = DBLE(ALPHA*TEMP1+
348: + DCONJG(ALPHA)*TEMP2)
349: ELSE
350: C(J,J) = BETA*DBLE(C(J,J)) +
351: + DBLE(ALPHA*TEMP1+
352: + DCONJG(ALPHA)*TEMP2)
353: END IF
354: ELSE
355: IF (BETA.EQ.DBLE(ZERO)) THEN
356: C(I,J) = ALPHA*TEMP1 + DCONJG(ALPHA)*TEMP2
357: ELSE
358: C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
359: + DCONJG(ALPHA)*TEMP2
360: END IF
361: END IF
362: 230 CONTINUE
363: 240 CONTINUE
364: END IF
365: END IF
366: *
367: RETURN
368: *
369: * End of ZHER2K.
370: *
371: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zher2k.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas