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