![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsyr2k.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 dsyr2k.f CVS log | Up to [local] / rpl / lapack / blas |
1.1 bertrand 1: SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
2: * .. Scalar Arguments ..
3: DOUBLE PRECISION ALPHA,BETA
4: INTEGER K,LDA,LDB,LDC,N
5: CHARACTER TRANS,UPLO
6: * ..
7: * .. Array Arguments ..
8: DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
9: * ..
10: *
11: * Purpose
12: * =======
13: *
14: * DSYR2K performs one of the symmetric rank 2k operations
15: *
16: * C := alpha*A*B' + alpha*B*A' + beta*C,
17: *
18: * or
19: *
20: * C := alpha*A'*B + alpha*B'*A + beta*C,
21: *
22: * where alpha and beta are scalars, C is an n by n symmetric matrix
23: * and A and B are n by k matrices in the first case and k by n
24: * matrices 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*B' + alpha*B*A' +
47: * beta*C.
48: *
49: * TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A +
50: * beta*C.
51: *
52: * TRANS = 'C' or 'c' C := alpha*A'*B + alpha*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 = 'T' or 't' or 'C' or 'c', K specifies the number
66: * of rows of the matrices A and B. K must be at least zero.
67: * Unchanged on exit.
68: *
69: * ALPHA - DOUBLE PRECISION.
70: * On entry, ALPHA specifies the scalar alpha.
71: * Unchanged on exit.
72: *
73: * A - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION 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 symmetric 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 symmetric 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: *
121: * LDC - INTEGER.
122: * On entry, LDC specifies the first dimension of C as declared
123: * in the calling (sub) program. LDC must be at least
124: * max( 1, n ).
125: * Unchanged on exit.
126: *
127: * Further Details
128: * ===============
129: *
130: * Level 3 Blas routine.
131: *
132: *
133: * -- Written on 8-February-1989.
134: * Jack Dongarra, Argonne National Laboratory.
135: * Iain Duff, AERE Harwell.
136: * Jeremy Du Croz, Numerical Algorithms Group Ltd.
137: * Sven Hammarling, Numerical Algorithms Group Ltd.
138: *
139: * =====================================================================
140: *
141: * .. External Functions ..
142: LOGICAL LSAME
143: EXTERNAL LSAME
144: * ..
145: * .. External Subroutines ..
146: EXTERNAL XERBLA
147: * ..
148: * .. Intrinsic Functions ..
149: INTRINSIC MAX
150: * ..
151: * .. Local Scalars ..
152: DOUBLE PRECISION TEMP1,TEMP2
153: INTEGER I,INFO,J,L,NROWA
154: LOGICAL UPPER
155: * ..
156: * .. Parameters ..
157: DOUBLE PRECISION ONE,ZERO
158: PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
159: * ..
160: *
161: * Test the input parameters.
162: *
163: IF (LSAME(TRANS,'N')) THEN
164: NROWA = N
165: ELSE
166: NROWA = K
167: END IF
168: UPPER = LSAME(UPLO,'U')
169: *
170: INFO = 0
171: IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
172: INFO = 1
173: ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
174: + (.NOT.LSAME(TRANS,'T')) .AND.
175: + (.NOT.LSAME(TRANS,'C'))) THEN
176: INFO = 2
177: ELSE IF (N.LT.0) THEN
178: INFO = 3
179: ELSE IF (K.LT.0) THEN
180: INFO = 4
181: ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
182: INFO = 7
183: ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
184: INFO = 9
185: ELSE IF (LDC.LT.MAX(1,N)) THEN
186: INFO = 12
187: END IF
188: IF (INFO.NE.0) THEN
189: CALL XERBLA('DSYR2K',INFO)
190: RETURN
191: END IF
192: *
193: * Quick return if possible.
194: *
195: IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
196: + (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
197: *
198: * And when alpha.eq.zero.
199: *
200: IF (ALPHA.EQ.ZERO) THEN
201: IF (UPPER) THEN
202: IF (BETA.EQ.ZERO) THEN
203: DO 20 J = 1,N
204: DO 10 I = 1,J
205: C(I,J) = ZERO
206: 10 CONTINUE
207: 20 CONTINUE
208: ELSE
209: DO 40 J = 1,N
210: DO 30 I = 1,J
211: C(I,J) = BETA*C(I,J)
212: 30 CONTINUE
213: 40 CONTINUE
214: END IF
215: ELSE
216: IF (BETA.EQ.ZERO) THEN
217: DO 60 J = 1,N
218: DO 50 I = J,N
219: C(I,J) = ZERO
220: 50 CONTINUE
221: 60 CONTINUE
222: ELSE
223: DO 80 J = 1,N
224: DO 70 I = J,N
225: C(I,J) = BETA*C(I,J)
226: 70 CONTINUE
227: 80 CONTINUE
228: END IF
229: END IF
230: RETURN
231: END IF
232: *
233: * Start the operations.
234: *
235: IF (LSAME(TRANS,'N')) THEN
236: *
237: * Form C := alpha*A*B' + alpha*B*A' + C.
238: *
239: IF (UPPER) THEN
240: DO 130 J = 1,N
241: IF (BETA.EQ.ZERO) THEN
242: DO 90 I = 1,J
243: C(I,J) = ZERO
244: 90 CONTINUE
245: ELSE IF (BETA.NE.ONE) THEN
246: DO 100 I = 1,J
247: C(I,J) = BETA*C(I,J)
248: 100 CONTINUE
249: END IF
250: DO 120 L = 1,K
251: IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
252: TEMP1 = ALPHA*B(J,L)
253: TEMP2 = ALPHA*A(J,L)
254: DO 110 I = 1,J
255: C(I,J) = C(I,J) + A(I,L)*TEMP1 +
256: + B(I,L)*TEMP2
257: 110 CONTINUE
258: END IF
259: 120 CONTINUE
260: 130 CONTINUE
261: ELSE
262: DO 180 J = 1,N
263: IF (BETA.EQ.ZERO) THEN
264: DO 140 I = J,N
265: C(I,J) = ZERO
266: 140 CONTINUE
267: ELSE IF (BETA.NE.ONE) THEN
268: DO 150 I = J,N
269: C(I,J) = BETA*C(I,J)
270: 150 CONTINUE
271: END IF
272: DO 170 L = 1,K
273: IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
274: TEMP1 = ALPHA*B(J,L)
275: TEMP2 = ALPHA*A(J,L)
276: DO 160 I = J,N
277: C(I,J) = C(I,J) + A(I,L)*TEMP1 +
278: + B(I,L)*TEMP2
279: 160 CONTINUE
280: END IF
281: 170 CONTINUE
282: 180 CONTINUE
283: END IF
284: ELSE
285: *
286: * Form C := alpha*A'*B + alpha*B'*A + C.
287: *
288: IF (UPPER) THEN
289: DO 210 J = 1,N
290: DO 200 I = 1,J
291: TEMP1 = ZERO
292: TEMP2 = ZERO
293: DO 190 L = 1,K
294: TEMP1 = TEMP1 + A(L,I)*B(L,J)
295: TEMP2 = TEMP2 + B(L,I)*A(L,J)
296: 190 CONTINUE
297: IF (BETA.EQ.ZERO) THEN
298: C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
299: ELSE
300: C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
301: + ALPHA*TEMP2
302: END IF
303: 200 CONTINUE
304: 210 CONTINUE
305: ELSE
306: DO 240 J = 1,N
307: DO 230 I = J,N
308: TEMP1 = ZERO
309: TEMP2 = ZERO
310: DO 220 L = 1,K
311: TEMP1 = TEMP1 + A(L,I)*B(L,J)
312: TEMP2 = TEMP2 + B(L,I)*A(L,J)
313: 220 CONTINUE
314: IF (BETA.EQ.ZERO) THEN
315: C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
316: ELSE
317: C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
318: + ALPHA*TEMP2
319: END IF
320: 230 CONTINUE
321: 240 CONTINUE
322: END IF
323: END IF
324: *
325: RETURN
326: *
327: * End of DSYR2K.
328: *
329: END